Finite element analysis of left ventricle during cardiac cycles in viscoelasticity

Computers in Biology and Medicine 75 (2016) 63–73

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Finite element analysis of left ventricle during cardiac cycles in viscoelasticity Jing Jin Shen a,n, Feng Yu Xu a, Wen An Yang b a b

School of Automation, Nanjing University of Posts and Telecommunications Nanjing, Jiangsu 210023, China College of Mechanical and Electrical Engineering, Nanjing University of Aeronautics and Astronautics Nanjing, Jiangsu 210016, China

art ic l e i nf o

a b s t r a c t

Article history: Received 3 December 2015 Received in revised form 18 May 2016 Accepted 20 May 2016

To investigate the effect of myocardial viscoeslasticity on heart function, this paper presents a finite element model based on a hyper-viscoelastic model for the passive myocardium and Hill's three-element model for the active contraction. The hyper-viscoelastic model considers the myocardium microstructure, while the active model is phenomenologically based on the combination of Hill's equation for the steady tetanized contraction and the specific time–length–force property of the myocardial muscle. To validate the finite element model, the end-diastole strains and the end-systole strain predicted by the model are compared with the experimental values in the literature. It is found that the proposed model not only can estimate well the pumping function of the heart, but also predicts the transverse shear strains. The finite element model is also applied to analyze the influence of viscoelasticity on the residual stresses in the myocardium. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Viscoelasticity Active contraction Finite element analysis Hill's equation

1. Introduction Cardiovascular diseases, being the leading cause of death and morbidity, have been studied from different aspects, such as diagnostic tests, surgical procedures, and medications. Both from clinical practice and laboratory research it has been realized that the mechanical performance of the ventricular muscle is one of the most important factors that can influence the pumping function of the heart. For instance, the weak contractive ability can cause a reduced ventricular ejection fraction, and the high stiffness of the myocardium muscle can cause impaired filling [43]. Thus, to better treat cardiovascular diseases, it is essential to investigate the relation between the heart's function and its mechanical properties. In general, stress and strain are two typical kinds of parameters in the mechanical characterization of local myocardial properties, and have been used as the vital determinants of many cardiac physiological and pathophysiological functions [15]. To avoid the experimental limitations encountered in the direct measurement of wall stresses in an intact ventricle [16], several mathematical models for the passive left ventricle (LV) have been presented to predict the stresses. Earlier mathematical models based on the shell analysis theory paid attention to the overall pressure–volume behavior of the LV. This kind of global model has been explored extensively and used in clinical diagnosis [20]. With n

Corresponding author. E-mail address: [email protected] (J.J. Shen).

http://dx.doi.org/10.1016/j.compbiomed.2016.05.012 0010-4825/& 2016 Elsevier Ltd. All rights reserved.

the development of continuum mechanics, especially the finite element method, mathematical models quantified in terms of local deformation and stress have been proven to be able to reflect some characteristics of the ventricular mechanics. Up to now, a few local-property models, including the pole–zero model [24], Fung-type models [30], and strain–energy function [4], have been proposed in the finite elasticity framework to describe the heterogeneous, anisotropic properties of the LV. Recently, a structurally based model of the passive myocardium, which takes into account the morphology and structure of the myocardium, was proposed by Holzapfel and Ogden [14]. Due to its simple invariantbased formulation and its small set of material parameters, the Holzapfel–Ogden model is particularly attractive in present practice [28]. In fact, as Fung mentioned, the myocardium is viscoelastic. Although the viscoelastic property of the myocardium is closely related to its physiological functions, especially in situations involving a high heart rate, the relevant study is still immature. To the author's knowledge, there are two viscoelastic myocardial models, both constructed in the quasi-linear viscoelasticity theory, that can be found in the literature [17,23]. To understand the active ventricular mechanics during systole, models characterizing the contractile properties of the cardiac muscle fibers are necessary. At present, several current models of cardiac muscle contraction under different assumptions, e.g., the time-varying elastance model, the modified Hill's model, and the fully history-dependent model, can be found in the literature. The time-varying elastance model describing the relation of the LV pressure to the LV volume and a slack volume reflects a load-

64

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

independent intrinsic contractility property of the LV [27]. The modified Hill's model predicts the fiber stresses by modifying shortening or lengthening of the sarcomere based on the force– velocity relation [4]. The fully history-dependent model proposed by [10] is a general deactivation model of cardiac contraction in rapid length perturbation. Based on active and passive models, various numerical techniques, such as the finite element (FE) method, have been used to estimate the deformation patterns of the LV under different circumstances. For instance, Vetter and McCulloch [40] developed an anatomically detailed FE model to analyze the stress state during passive filling of the left ventricle. Bovendeerd et al. [2] developed an FE model in which the LV is approximated by a truncated confocal ellipsoid, and studied the influence of the fiber orientation on the LV mechanics. Dorri et al. [5] constructed an FE model for a representative human ventricle to investigate the role played by the LV myocytes in cardiac function. Jiang et al. [19] used the smoothed finite element method to guarantee the accuracy of the analysis of the passive LV to be volumetric locking-free. As an important functional indicator of various biological organs, residual stresses/strains have been extensively investigated from the aspects of experimental measurement and theoretical analysis. Rachev and Greenwald [29] reviewed the experimental methods for measuring the residual stresses in the artery wall, and found that residual stresses can increase arterial compliance, thereby making it more suitable to be a blood conduit. The residual strains in the myocardium were well illustrated by Omens and Fung [25], who radially cut an equatorial slice of rat LV, and defined the opening angle as a residual strain index. Summerour et al. [37] followed Omens and Fung's line to study how the remodelings of collagen fibers and myocytes change the residual strains in the ischemic ventricular myocardium. Costa et al. [4] developed an experimental procedure to measure the transmural distribution of the residual strains by a biplane radiographic scanning of the coordinates of beads implanted in the heart wall. To incorporate residual stresses into the constitutive modelling, Nash and Hunter [24] introduced a growth tensor that can map the load-free configuration to the stress-free configuration in the analysis of the beating heart. Alternatively, Wang et al. [42] applied the virtual configuration approach proposed by Hoger [12] to an FE analysis based on the Holzapfel–Ogden model. In this paper, efforts are made to combine an orthotropic finiteviscoelastic model for the passive myocardium developed by the authors and a modified Hill's model for the active contraction to analyze the deformations and stresses experienced by the LV over the whole of the cardiac cycles. Special attention is paid to the residual stresses and the end-diastolic strains, whose shear components are hard to correctly predict by an elastic constitutive model. In addition, to improve the stability and convergence of the FE analysis, the consistent tangential moduli of the modified Hill model are derived and implemented.

2. A hyper-viscoelastic model of the passive myocardium The viscoelastic model of the passive myocardium used to analyze the LV is briefly described in this section. For its detailed derivation and application to a cylindrical LV model, cf. [31].

Fig. 1. Schematic view of the fiber distributions across the LV.

myocardial microstructural architecture [22] and the slippage of adjacent muscle layers along the cleavage planes [21]. Therefore, there are three mutually orthogonal preferred material directions: the fibre axis f0 , the sheet axis s0 , and the sheet normal axis n 0 , as shown in Fig. 1. It is well known that the isotropic constitutive relation merely depends on the deformation quantities, but the anisotropic relation is simultaneously affected by both the deformation quantities and the intrinsical material structural quantities. In general, the structural quantities can be represented by the structural tensors defined by the tensor products of the preferred direction vectors. In orthotropicity, three irreducible structural tensors, L 0 , M0 , and N0 , can be expressed as

L 0 = f0 ⊗ f0,

(1)

N0 = n 0 ⊗ n 0

The symbol ⊗ denotes the tensor product of two vectors. To establish the thermodynamically viscoelastic constitutive relation, the Helmholtz free energy function Ψ, which is the thermodynamical counterpart of the elastic strain energy, has to be suitably postulated for the specific material. In this paper, if we only consider isothermal mechanical phenomena, then orthotropicity leads to Ψ = Ψ (C, L 0, M0, N0, Q1, … , Q q), which is a scalarvalued function with 4 + q tensor arguments. Among these, C is the right Cauchy–Green strain tensor, and the internal variables Q1, …, Q q , defined in the reference configuration, represent the dissipative mechanism associated with the viscoelasticity. Based on the experimental evidence that the bulk deformation of soft tissue is elastic, and the viscoelasticity mainly occurs in the deviatoric deformation, the free energy function can be decoupled in the physical sense as [13] 1 ^ Ψ (C, L 0, M0 , N0 , Q 1, … , Q q) = Ψ 0 (J ) + Ψ 0 (C, L 0, M0 , N0 , Q 1, … , Q q) − vol iso 2

q

∑ i =1

⎛ q ⎞ ^ C: Q i + Ξ ⎜⎜ ∑ Q i⎟⎟ ⎝ i =1 ⎠

(2)

^ where J is the Jacobian determinant of the deformation gradient, C 0 is the initial freeis the isochoric right Cauchy–Green tensor, Ψvol 0 describes energy function describing the bulk elastic response, Ψiso the deviatoric elastic response, and Ξ is a given function of the internal variables, which can be determined via the Clausius–Plank inequality. The reason for the introduction of the specific viscoen ^ 1 lastic form, − 2 ∑i = 1 C: Q i , in the free-energy function is twofold. First, it is a simple matter to find a function Ξ that makes Ψ satisfy the second law of thermodynamics. Second, this kind of free energy function can be viewed as the finite deformation extension of the linear Maxwell viscoelastic model. Through the Coleman–Noll procedure, the second Piola–Kirchhoff stress and the free-energy function have the relation

2.1. The decoupled free-energy function From a morphological and microstructural standpoint, Holzapfel and Ogden [14] have pointed out that the passive myocardium tissue can be locally characterized as an orthotropic material in the reference configuration. The orthotropy of the myocardium is also illustrated by the experiments about the

M 0 = s0 ⊗ s0,

S=J

⎛ ∂Ψ 0 0 ∂Ψvol C−1 + J−2/3 : ⎜⎜ 2 iso − ^ ∂J ⎝ ∂C

with 1

 =  − 3 C −1 ⊗ C

q

⎞

i=1

⎠

∑ Q i⎟⎟

(3)

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

Table 1 Values of parameters in the hyper-viscoelastic model of passive myocardium.

Q ij (t ) = 2

a1

b1

a4s

a4f

b4s

b4f

a8fs

b8fs

γf

γs

τf

τs

0.93

4.34

16.13

136.3

11.9

13.31

3.91

5.24

0.71

0.61

0.13

0.39

where  is the fourth-rank identity tensor and  is a projection tensor which maps a material quantity to its deviatoric component. In agreement with the elastic model proposed for the passive myocardium in [14], the exponential function is adopted as the free-energy function 0 Ψvol = c (J − 1)2

0 Ψiso =

(4)

a exp [b (I1 − 3)] + 2b

∑ j=f ,s

aj {exp [bj (I4j − 1)2] − 1} 2bj

afs + [exp (bfs I82fs ) − 1] 2bfs

(5)

where c, a, b, af, bf, as, bs, afn, bfn, afs and bfs are positive material constants whose values are shown in Table 1.

2.2. Evolution rule of the internal variables Around the myocardial myocytes, there is a complex weave of collagen fibers. Moreover, a couple of the weaves are loosely connected to form long collagen fibers. The weave structure is thought to be the origin of the myocardium viscoelastic properties [2]. Accordingly, we surmise that the viscoelastic property of the myocardium is caused by the deformation of the collagens and muscle fibres. That is, the evolution of the internal variables are only linked to the invariants I4f and I4s . This assumption can reduce the number of material parameters in the stress–strain relation, making the parameter identification procedure more reliable. Further, we assume that the effects imposed on each internal variable, for example Q i , by the collagens and muscle fibres are uncoupled. Therefore, Q i has the expression

Q i = Q if + Q is,

i = 1, …, q

(6)

where the evolution of the component Q if is governed by I4f , while that of Q if is governed by I4s . A variety of evolution rules for the internal variables can be established based on the specific relations connecting Q if with I4f and Q is with I4s . In the present paper, the linear differential relation proposed by Simo and Hughes [36] is adopted, namely, 0 ∂I γij ∂Ψiso 1 4j Q̇ ij (t ) + Q ij (t ) = 2 : , ^ ∂I4j ∂C τij τij

j = f, s

0 ∂I ∂Ψiso 4j ^ ∂I4j ∂C

τij

t

⎛

⎞

∫−∞ exp ⎜⎜ wτ− t ⎟⎟ : Φj

⎡ dw = 2γij ⎢ : Φj − ⎢⎣

⎝ t

∫−∞

ij

(8)

Upon introducing the integrating factor exp (t /τij ), the first-order differential linear evolution rule in Eq. (7) can be integrated:

⎠

⎛w − exp ⎜ ⎝ τij

⎤ t⎞ ⎟ d (: Φj ) ⎥ ⎠ ⎦⎥

(9)

The last term inside the brackets in the above expression is evaluated via integration by parts. After substituting Eq. (9) into Eq. (3), the second Piola–Kirchhoff stress can be further expressed as 0 ∂Ψvol C−1 + 2J−2/3  ∂J ⎧ q ⎡ ⎪ ∂Ψ 0 : ⎨ iso − ∑ ∑ γij ⎢ Φj − ^ ⎢⎣ ⎪ i=1 j=f ,s ⎩ ∂C

S=J

t

⎛

⎞

⎤⎫ ⎪

∫−∞ exp ⎜⎝ wτ− t ⎟⎠ d (: Φj ) ⎥⎥ ⎬⎪ ⎦⎭

ij

(10)

If the material under investigation is strictly incompressible, the internal constraint J ¼1 must be taken into account. As a result, the rate of change of the Cauchy–Green deformation tensor Ċ is no longer arbitrary, which makes the classical Coleman–Noll procedure fail. However, if the free-energy function Φ (C) is linearly homogeneous, i.e., if Φ (α C) = αΦ (C), where α is an arbitrary constant, the second Piola–Kirchhoff stress for incompressible materials can be expressed as

S = − pJ C−1 + 2J−2/3  ⎧ q ⎡ ⎪ ∂Ψ 0 : ⎨ iso − ∑ ∑ γij ⎢ Φj − ^ ⎢⎣ ⎪ i=1 j=f ,s ⎩ ∂C

t

⎛

⎞

⎤⎫ ⎪

∫−∞ exp ⎜⎝ wτ− t ⎟⎠ d (: Φj ) ⎥⎥ ⎬⎪ ij

⎦⎭

(11)

in which p denotes the hydrostatic pressure, which can only be determined from the equilibrium equations and the boundary conditions. In the spatial configuration, the Cauchy stress σ can be obtained by applying the push-forward operation, that is, σ = FSFT /J . In practice, we can find that the constitutive relation with one internal variable suffices to approximate well the experimental data. And the corresponding relaxation times and the ratios of the relaxation modulus to the instantaneous modulus are summarized in Table 1.

3. Constitutive model for active contraction 3.1. Hill's three-element model Hill's equation was originally proposed to characterize the contractile dynamics of a steady tetanized skeleton muscle during a quick-release experiment. Complying with the thermomechanical laws, Hill's equation relating the tension in the muscle, Ph, with the velocity of contraction, v, can be written as

(v + b)(Ph + a) = b (Ph0 + a),

(7)

with the initial condition Q ij ( − ∞) = 0 , τij ≥ 0 denotes the relaxation time, and γij ranging from 0 to 1 represents the ratio of the relaxation modulus to the instantaneous modulus. This evolution rule is motivated by the Maxwell model in linear viscoelasticity; it has been successfully applied in rubber and tissue modelling. For convenience, we define henceforth

Φj =

γij

65

b = av0/Ph0

(12)

where Ph0 represents the isometric maximum tension corresponding to v ¼ 0 and a is the coefficient of shortening heat. In order to extend the application of Hill's equation to the myocardial contraction, two distinct features of the cardiac muscle have to be considered: (1) the cardiac muscle can only be excited by single twitches; (2) the myocardium in the resting state is relatively stiff. Therefore, the myocardial active state is a transient process, and the passive stresses cannot be neglected in the circumstances of large active stresses. To describe the mechanical behavior of a muscle, the most commonly used model is Hill's three-element model shown in Fig. 2, which consists of a contractile element, a series element, and a parallel element. It is worth noting that the arrangement of the three elements is not

66

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

3.2. An expression for the active state function A (ts, ls, lc ) In this paper, the dependence of A on (ts, ls, lc ) is assumed to have the following decoupled form:

A (ts, ls ) = A1 (lc ) A2 (ts, ls )

Fig. 2. Schematic of Hill's three-element model.

unique. In this paper, we associate the parallel viscoelastic element to the mechanical property of the passive myocardium discussed in Section 2, the series elastic element to the interaction between the actin/myosin molecules and cross-bridges, and Hill's element to the sarcomere. Since it is very hard, if not impossible, to experimentally distinguish the influence due to the series elastic element from that due to Hill's element, these two elements often are considered as one active unit. As mentioned above, to describe the length–time-dependent property of the myocardial activation, it is necessary to modify the classic Hill equation. A simple and efficient modification approach is to introduce an active state function An which depends on the stimulation time ts, the length of the sarcomere ls, and the length of the contractile element lc [1]. Then Hill's equation is modified to

⎛ P ⎞⎛ v⎞ ⎜ 1 + h ⎟ ⎜ 1 + ⎟ = A⁎ (ts, ls , lc ) + 1 cp ⎠ ⎝ cv ⎠ ⎝

(13)

where cp and cv are constant parameters. Because the lengths of lc and ls are very close in common cases, hence, it is acceptable to use ls as the length of the active unit. As can be seen in Fig. 2, in equilibrium, the stresses generated by Hill's element and the series elastic element are equal. Hence, the active stress can be characterized by the stiffness of the elastic element and its change in length ls − lc . To make the analysis tractable, we make the following assumptions, whose validation can be found in [2]: (1) the deformation of the cardiac muscle at each time instant is composed of one isosarcomeric contraction followed by one quickrelease contraction; (2) the relaxation of the contractile state due to an abrupt length change is negligible; (3) given the sarcomere length, the change of length of the series elastic element at the isometric maximum tension can be approximated by a constant, Δl . Under these assumptions, the relation of the active stress with the change in length has the following simple form: ⁎

Ph = (ls − lc )

cp A P h0 = (ls − lc ) Δl Δl

(14)

At this stage, if the active stress is assumed to be along the direction of the muscle fiber, then the first Piola–Kirchoff active stress P can be obtained by making use of Eq. (14): ⁎

P = Ea A (ls − lc ) f ⊗ f0 Ea = cp/Δl

with

⁎

A = Pmax A (ts, ls , lc ), (15)

(16)

where two normalized functions, A1 and A2, reflect, respectively, the length-dependence and time-dependence properties of the active contraction. To determine the formulation of A1 (lc ), it is worth understanding the difference between isometric contraction and isosarcometric contraction. In isometric contraction, the length of the active unit, including Hill's element, and the series elastic element is fixed. In the isosarcometric contraction, the length of Hill's element is fixed. Strictly speaking, the function A1 (ls ) reflecting the stress-length dependence of the sarcomere must be established by isosarcometric experiments. However, few experiments on myocardial muscle contraction are isosarcometric. In an isosarcometric experiment conducted in [38], it is found that the maximum contractile stress linearly increases along with the sarcomere length ls ranging from 1.5 μm to 2.1 μm, the maximum stress levels off for ls > 2.1 μm , and no decreasing phenomenon presents in the physiological length range of the sarcomere. Thus, a relation approximating the experimental data can be formulated as

A1 (lc ) = 1 + ax (lc − lx ) −

(ax (lc − lx ))2 + 0.01

(17)

where ax is a constant coefficient. To determine A2 (ts, ls ), the time courses of the isosarcomeric contractions with different sarcomere lengths need to be considered. In a study of the heat production of isolated trabecula, Han et al. [2012] found that the length-dependence primarily influences the decaying session of the active state, with a minor contribution for the rising counterpart. Therefore, we can express A2 (ts, ls ) in a multiplicative form as

A2 (ts, ls ) = Br (ts ) Bd (ts, ls )

(18)

where Br represents the rising session approximated by

Br (ts ) = 1 −

1 1 + (ts/tr )4

(19)

and tr is the time interval of the rising session, which is assumed to be constant in the model. On the other hand, Bd is approximated by

Bd (ts, ls ) = 1 −

1 , 1 + ((te − ts )/td )4

ts ≤ te

(20)

where td is the time interval of the decaying session, and te denotes the total duration of the activation. In [11,39], it is verified that the twitch duration is linear with the sarcomere length. Therefore, the following linear relation between te and ls holds:

te = b (ls − ld )

(21)

where b is a constant coefficient and ld is the extrapolated sarcomere length at which te ¼0. To determine lc, we can find the time derivative of lc by transforming Eq. (13) as ⁎

where f0 and f are the unit vectors in the direction of the fiber ⁎ before and after deformation, respectively, and Ea A represents the stiffness of the series elastic element. In the stiffness term, An is expressed as the multiplication of the maximum isometric stress Pmax by the normalized active function A. In Section 3.2, we will determine a specific expression for A (ts, ls, lc ) by decoupling it into two parts, the length-dependent part and the time-dependent part.

cp A − Ph dl c cv = dt cp + Ph

(22)

The initial value of lc at ts ¼0 is the sarcomere length ls at ts ¼ 0. 3.3. The consistent tangent moduli of the active constitutive relation As illustrated by Eq. (16), ls and lc at any instant are essential parameters in the numerical simulation of the heartbeat. Let γ be

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

the initial sarcomere length. Then, the sarcomere length in a deformed configuration with deformation gradient F is

ls =

γ Ff0·γ Ff0 = γ Fij f0 j Fik f0 k

(23)

where the Einstein summation convention is used for the repeated indices i, j, k which can range within the set {1, 2, 3}. To determine the instantaneous length of the contractile element lc, we need to integrate the velocity equation (22). As illustrated in Eq. (15), the developed active force is irrelevant for the velocity of the contraction. Therefore, the constitutive relation governing the active contraction is of the time-dependent, rate-independent type. In this case, the so-called consistent or algorithmic tangent moduli, defined as the differentiation of the discretized constitutive equation with respect to the strain increment, are advantageous in numerical analysis, which can make the iterative update solution scheme achieve an asymptotic quadric rate of convergence [8]. In this section, the derivation of the consistent tangent moduli of the active constitutive relation is based on the integration of Eq. (22) through the second-order Adams–Bashforth–Moulton method. Given the values of all state variables in the constitutive relation, lc is calculated by the following procedure. By Eq. (22), the approximation of lc by the Adams–Bashforth predictor can be explicitly expressed as

lc n + 1 = lc n −

⁎ ⁎ Ph n − 1 − cp A n − 1 ⎞ Δt ⎛ Ph n − cp A n ⎜3 ⎟ cv − 2 ⎝ Ph n + cp Ph n − 1 + cp ⎠

(24)

With the Adams–Moulton corrector, lc expressed in Eq. (24) can be corrected in the following implicit form as

lc n + 1 = lc n −

Δt ⎛ Ph n − ⎜ 2 ⎝ Ph n + cp

⁎ cp A n

+

⁎ cp A n + 1 ⎞

Ph n + 1 + Ph n + 1 + cp

⎟ cv ⎠

(25)

where the indices n − 1, n and n + 1 of a physical quantity refer, respectively, to its values at the time points tn  1, tn and tn + 1, and Δt = tn + 1 − tn is the time step. To simplify the notation, we define C1, C2 and C3 as ⁎

Δt Ph n − cp A n C1 = lc n − cv 2 Ph n + cp cp Δt = − cv 2 Ph n + 1 + cp

Δt Ph n + 1 C2 = − cv 2 Ph n + 1 + cp

C3

lc n + 1 = C1 + C2 +

(26)

where C1 is a constant parameter that can be calculated from the known variables at tn, while C2 and C3 are functions of Ph n + 1. To ⁎ determine the discrete formulation of An + 1 at tn + 1, substituting Eqs. (17) and (20) into Eq. (15) leads to ⁎ A n + 1 = C4 n + 1 ⎡⎣ 1 + ax (lc n + 1 − lx ) −

(ax (lc n + 1 − lx ))2 + 0.01 ⎤⎦

(27)

with

C4 n + 1 = Pmax Br (ts n + 1) Bd (ts n + 1, ls n + 1)

(28)

In the framework of continuum mechanics, the tangential moduli of the constitutive Eq. (15) based on the first Piola– Kirchhoff stress can be obtained by

=

∂P ∂P ∂f ⊗ f0 = h f ⊗ f0 + Ph ∂F ∂F ∂F

⁎ ⎛ ∂l ⎞ ∂A n + 1 ∂Ph n + 1 ∂l ⁎ = Ea (ls n + 1 − lc n + 1) + Ea A n + 1 ⎜ s n + 1 − c n + 1 ⎟ ∂Fij n + 1 ∂Fij n + 1 F F ∂ ∂ ⎝ ij n + 1 ij n + 1 ⎠

(29)

where  is the tangential moduli. The term ∂(f0 ⊗ f) /∂F in Eq. (29) can be widely encountered in soft tissue modeling and is well elaborated in [13]. In the following, attention is paid to deriving the derivative ∂Ph/∂F . From Eq. (15), the temporal discretization of ∂Ph/∂F can be expressed as

(30)

All partial derivative terms in the right-hand side of the above equation can be calculated as follows: First, making use of Eq. (23), the partial derivative of ls with respect to F has the following analytical expression:

ls f Fmk f0 k ∂ls l ∂(Fij f0 j Fik f0 k ) = s = 0n ∂Fmn ∂Fmn 2γ γ

(31)

Second, the partial derivative of lc with respect to F can be calculated from Eq. (26): ⁎

⁎ ∂A n + 1 ∂lc n + 1 Δtcv cp (1 + A n + 1) ∂Ph n + 1 = − + C3 2 2 (Ph n + 1 + cp ) ∂Fij n + 1 ∂Fij n + 1 ∂Fij n + 1 ⁎

= C3

∂A n + 1 ∂P + C4 h n + 1 ∂Fij n + 1 ∂Fij n + 1

(32)

with ⁎

C4 = −

Δtcv cp (1 + A n + 1) 2 (Ph n + 1 + cp )2

Third, taking into account Eqs. (15), (16) and (18), the partial derivative of An with respect to F can be explicitly expressed as ⁎ ⎡ ∂A ∂A (l ) ∂lc ∂B (t , l ) ∂ls ⎤ ⎥ = Pmax Br (ts ) ⎢ Bd (ts, ls ) 1 c + A1 (lc ) d s s ∂Fij ∂lc ∂Fij ∂ls ∂Fij ⎦ ⎣

(33)

Furthermore, from Eqs. (17) and (20), two terms in Eq. (33) can be obtained as

∂A1 (lc ) = ax − ∂lc

ax (lc − lx ) 0.01 + ax (lc − lx

)2

4btd4 (te − ts )3 ∂Bd (ts, ls ) = 2 ∂ls td4 + (te − ts )4

(

)

Then, by substituting Eq. (32) into Eq. (33), we can get

⁎

∂A n + 1 ∂Fij n + 1

⎡ ⎤ ∂B ∂l ∂A ∂P Pmax Br n + 1 ⎢ A1 n + 1 d n + 1 s n + 1 + C4 Bd n + 1 1 n + 1 h n + 1 ⎥ l F l F ∂ ∂ ∂ ∂ ⎣ s n+ 1 ij n + 1 c n+ 1 ij n + 1 ⎦ = ∂A1 n + 1 (1 − C3 Pmax Br n + 1Bd n + 1) ∂lc n + 1

Then, Eq. (25) can be written as ⁎ C3 A n + 1

67

(34)

At this stage, the value of ∂Ph/∂F can be easily evaluated by substituting Eqs. (31), (32) and (34) into Eq. (30). Now the consistent tangential modulus  can be calculated by a routine procedure [13]. For different finite element formulations, the above consistent tangent moduli need to be transformed to the appropriate form. For example, the consistent tangent moduli based on the second Piola–Kirchhoff stress tensor and the Cauchy stress tensor are necessary for the total Lagrangian formulation and updated Lagrangian formulation, respectively. In essence, the transformation between different forms of the consistent tangent modulus is a kind of push-forward or pull-back operations in differential geometry [36].

4. Finite element implementation 4.1. A geometrical model of the LV In order to improve the computational efficiency and stability, this paper adopts a simplified axisymmetric geometry instead of a more complex realistic geometry. As in [9], the three-dimensional

68

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

Fig. 3. LV model with the Cartesian coordinates (x1, x2, x3 ) and the prolate spheroidal coordinates (η1, η2, η3 ) , fiber angle α defined in (η2, η3 ) , sheet angle β defined in the (η1, η2 ) plane.

anatomical measurements of the anterior free wall of the diastolic canine heart are used to establish the reference geometric model of the LV. In [9], a finite element model of the LV with an elastance model for the active contraction and an elastically transversely isotropic model for the passive myocardium was established in prolate spheroidal coordinates (η1, η2, η3 ). Although prolate spheroidal coordinates have the advantage of approximating the LV geometry with a small set of approximation parameters, it may be not convenient for some practical situations. For example, the variational form of the contact between the tactile sensor and the LV is hard to develop in prolate spheroidal coordinates, which is an important part of our research [31,33–35]. As a result, the Cartesian coordinates (x1, x2, x3 ) shown in Fig. 3 are chosen in this paper to describe the geometry of the LV. Mathematically, the following transformation equations hold between prolate spheroidal coordinates and Cartesian coordinates: x1 = d sinh η1 sin η2 cos η3, x2 = d sinh η1 sin η2 sin η3, x3 = d cosh η1 cos η2. Based on the specific values of (η1, η2, η3 ) describing the endocardial and epicardial borders in [9], we establish a geometrical model of LV with the focal lengths d set to be 29 mm. In addition, there is a small hole at the apex of the LV (η3 = [0°, 3°]) to avoid the numerical singularity due to the element degeneration. 4.2. Material property and fiber orientation With regard to the passive mechanical deformations, the local material property of any point of the LV is characterized by the orthotropic hyper-viscoelastic model, and its parameters are determined by fitting to the experimental data of simple shear deformations and shown in Table 1 [31,32]. On the other hand, the parameters involved in the active contraction model are determined by fitting to the experimental data reported in [18], which include the force–time diagram under sarcomere length control, duration of systolic phase in sarcomere-isometric contractions, time course and duration of final relaxation phase in sarcomere-isometric contraction. In the fitting process, a genetic algorithm is used to solve the multi-parameter optimization problem, and the initial value of each parameter is chosen to be that given in [2]. The obtained values are summarized in Table 2. Furthermore, the initial length of the sarcomere is set to be 2.0 μm, and is assumed to be uniformly distributed from endocardium to epicardium.

To determine the fiber/sheet orientation in the myocardium, we use a rule-based approach which can generate a homogeneous fiber orientation [6]. In the rule-based approach, the fiber angle, denoted by α, cubically varies along the transmural direction, from α = − 60° at the endocardium to circumferential at the midwall to α = 60° at the epicardium. The sheet angle β, defined in the (η1, η2 ) plane, linearly varies from −85° at the endocardial border to 85° at the epicardial border.

4.3. Boundary and loading conditions Boundary conditions are necessary to eliminate the rigid-body displacement of the LV during a finite element analysis. For a normal heart, the outer border of the LV is covered by the pericardial sac, and its basal part is constrained by the aorta. Furthermore, there are interactions between the LV and the right ventricle, mitral valve, etc. These constraint factors are too complex to be quantitatively described. Therefore, reasonable approximate conditions are implemented in the finite element analysis. The displacements of all nodes at the base between the epicardium and the midwall are fixedly constrained, and the other nodes can freely move. Besides the above displacement boundary conditions, the pressure interaction between the intraventricular blood and the endocardium must be represented as loading boundary conditions. In general, the complex fluid dynamics of the blood needs to be considered to determine the distribution of the blood pressure, which is dependent on location and time. Taking into consideration that the pressure gradient in the LV is much less than its absolute value during ejection, we assume that the blood pressure is a function of time by neglecting its spatial dependence. This assumption can greatly reduce the computational workload in the finite element analysis. Thus, a typical pressure–time diagram in [7,31] is used in the present paper. In addition, the initial conditions must be appropriately defined for a viscoelastic analysis. Without loss of generality, we assume that the myocardium starts beating from a natural state, i.e., no residual stresses and strains. Moreover, the values of all the state variables in the passive viscoelastic constitutive equation are set to be zero.

Table 2 Values of parameters in active contraction model. T1 (kPa)

Ea (μm  1)

ax (μm  1)

lx (μm)

ld (μm)

tr (ms)

td (ms)

cv (μm s  1)

b (s mm−1)

122.68

18.7

1.46

2.57

−0.37

71.74

78.41

7.53

143.39

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

69

Fig. 4. Finite element model of the LV, the geometrical surfaces of the endocardium, midwall and epicardium.

4.4. Finite element discretization After specifying the LV geometry, its material properties, and the initial/boundary conditions, the discrete model of the LV, consisting of multiple unidirectional fiber lamina, was developed in the nonlinear finite element software, ABAQUS FEA (Simulia. Providence, RI, USA). The material properties of the myocardium developed in Sections 2 and 3 were implemented in the finite element model via a UMAT subroutine. The element for discretizing the LV is chosen to be C3D8R which is an eight-node brick with reduced integration and hourglass control. Inside each element, ten equal-thickness composite lamina whose muscle fiber/sheet directions were defined in Section 4.2 are created. In total, the LV is discretized to be a 234-element mesh as shown in Fig. 4. Comparing with standard solid-element model, the composite continuum-element model has more than 50% savings in degrees of freedom. When implementing the model in a laboratory computer (3.60 GHz CPU, 16 GB RAM), the computation time for one cardiac cycle gradually decreases as the simulation proceeds, with less than 9 h for the first cardiac cycle. To achieve a balance between computational convergence and efficiency, the model was solved by the ABAQUS/standard iterative solution strategy, which uses the consistent tangent moduli. However, it should be noted that the consistent tangent moduli of the active contraction model are nonsymmetric, hence, the solver for the non-symmetric matrix equations is suitable.

5. Results In this section, the transmural distribution of the end-diastolic and end-systolic strains are computed and compared with the excremental measurements to validate the proposed finite element model. Then, based on the model, we investigate the influence of the viscoelasticity on the pressure–volume relation of the heart and the residual stresses in the myocardium. 5.1. End-diastole strains By using biplane radiography of three transmural columns of radiopaque beads, Omens et al. [26] measured the end-diastolic strains with reference to the zero-pressure state in the equatorial area of seven isolated dog hearts. When the left ventricular enddiastolic pressure is 1.1 ± 0.5 kPa , the normal strain components in cardiac coordinates increase in magnitude from the epicardium to the endocardium. Furthermore, the circumferential strain Ecc and the longitudinal strain Ell remain positive through the LV wall thickness, whereas the radial Err is negative. This kind of strain distribution is the underlying microstructural manifestation of wall thinning. In contrast to the normal strain components, the in-plane shear strains have small magnitudes and transmural gradients. The estimation of the strain components from the FE model in this paper and the experimental measurements performed in [26] are shown in Fig. 5.

From Fig. 5, it can be seen that all strain components are within the measurement range. Although the change trend of the shear strain component Ecr is different from the measurement data, it is better than that predicted in [24], where Ecr is more than ten times less than the measurement value. At the same time, the distribution of another shear strain component Ecl has a negative magnitude. This is consistent with the torsional deformation along the apex–base direction of the LV during diastolic filling, that is, the apex rotates clockwise with respect to a fixed base, as viewed from the apex. In addition, the magnitude decrease of Ecl from the endocardium to the epicardium reflects how the torsion decreases with wall depth, which is also consistent with the experimental observations. 5.2. End-systole strains With cineradiography, Waldman et al. [41] examined the threedimensional deformations for seven open-chest dogs during contraction, and computed the end-systole strains with respect to the end-diastole configuration. At the equator, the transmural distributions of the end-systole strains from [41] and those predicted by the proposed finite element model are illustrated in Fig. 6. From this figure, we can see that the normal components match the experimental values better than the shear components in an overall sense. Especially, the in-plane shear strain Ecl apparently lies outside the range of the corresponding experimental values. The same discrepancy also appears in the analysis performed by Guccione et al. [9]. One of the advantages of the proposed finite element model is that the prediction of the circumferential radial strain Ecr has the same variation tendency as the experimental data. Hence, the model can be used to estimate the apex–base rotation. 5.3. Effect of viscoelasticity on myocardium contraction Pressure–Volume (P–V) loops reflects the global mechanical functioning of the heart. To study the effect of viscoelasticity, the P–V loops have been computed for the 11th, 7th, and 3rd cardiac cycles, and shown in Fig. 7(a). It can be seen that both the endsystole and end-diastole volume increases from the initial cardiac cycle and achieves a steady state after 11 cycles. This is a manifestation of myocardial viscoelastic creep. Furthermore, the volume increase at the end diastole is greater than that at the end systole, which can be supported by the change of the strain, Ecc, in Fig. 7(b). 5.4. Residual stresses Due to the viscoelastic property of the myocardium, the deformation/stress responses of the LV to periodic variation of blood pressure should gradually tend to be steady. Through simulations, it is found that the residual stresses corresponding to zero blood pressure achieve a maximum value, then level off after 11 beating

70

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

Fig. 5. Transmural distributions of the strain components around the equatorial area at end-diastole.

cycles. The effect imposed by the viscoelastic property of the myocardium on LV passive filling has been investigated in [31]. In that paper, the residual stresses at the equator area predicted by the viscoelastic model and Guccione's elastic model are compared and it was found that its shear components are much less than the normal components. As a further development of Guccione's model, Nash's model proposed in [24] can predict the residual stresses in the fiber direction over all regions of LV. In the present paper, the residual stresses in the fiber, sheet, and normal to sheet directions are calculated by the developed finite element model. Fig. 8 shows the calculated results and the fiber residual stresses predicted by Nash's model at the base, the equator and the apex of the LV. As can be seen in Fig. 8, there are apparent discrepancies in the fiber residual stresses predicted by the FEA and those by Nash's model. Hence, we can state that other factors besides the viscoelasticity of the LV, such as tissue growth, should play a more important role in the generation of the residual stresses. In addition, the magnitudes of all normal residual stresses predicted by the FEA are of the same order. This means that the residual stresses caused by viscoelasticity are nondirectionally dominate. Fig. 8 also shows the dependence of the residual stresses on location, e.g., its

magnitude around the apex is much greater than that around the equator and the base. This difference can be explained by the fact that the apex of the LV is subjected to less constraints and so experiences a larger deformation than other areas. As another important point, the residual stresses at the equator and base predicted by the FEA are very small. This is not consistent with the results for passive filling of the LV computed by [31], in which large values were obtained. This inconsistency can be attributed to the active contraction. Without the active fiber stress, the largest deformation of the LV occurs around the peak value of the blood pressure, while if the active stresses are considered, the largest deformation occurs at end-diastole. Furthermore, the largest deformation without the active stresses is ten times greater than its counterpart with the active stresses. Therefore, taking account of the memory effect of viscoelasticity, the residual stresses estimated by the FEA should be much smaller than those in [31]. 6. Summary and discussion In the present paper, a finite element model was built on the viscoelastic model of the passive myocardium and a model based

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

Fig. 6. Transmural distributions of the strain components around the equatorial area at end-systole.

Fig. 7. Influence of viscoelasticity on the LV P–V and the circumferential strain Ecc at the endocardium of the equator.

71

72

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

geometrical and boundary conditions. Since the viscoelastic model of the myocardium in this paper is determined by fitting to the experimental data of simple shear deformations, the shear strains predicted by the proposed model are apparently better than those predicted by other models in the literature. However, it should be noted that the trend of Ecr from the endocardium to the epicardium does not match the experimental curve. Therefore, to further improve the prediction of Ecr, both the geometrical and boundary conditions must be refined. Using this model, the three-dimensional residual stresses/ strains were investigated in the ventricular wall. As mentioned above, the residual stresses/strains can be produced by nonuniform growth, turnover of LV constituents, or other means. By comparing the residual stresses estimated by this paper with those by other models in the literature, we have found that the residual stresses due to the myocardial viscoelasticity are very small in the basal and equatorial areas. In contrast, around the apical area, there is a big discrepancy between the viscoelastic residual stresses and their elastic counterparts from Nash's model. Therefore, this suggests that the traditional method for measuring the residual strain, such as the open angle method, may not be suitable for this area. Finally, the proposed finite element model is subject to several other limitations. For instance, the interaction between the left ventricle, right ventricle and papillary muscle are not considered.

Conflict of interest statement None declared.

Acknowledgements The work was supported by the Scientific Research Foundation of NJUPT (NY214020), Science and Technology Support Program of Jiangsu Province (CN)(BK20150844) and Natural Science Foundation of the Jiangsu Higher Education Institutions of China (15KJB310010). These research funds sponsor the author in the development of computational model and the following simulations. Fig. 8. The transverse distribution of residual stresses; sff, sss and snn, respectively, denote the residual stresses in the fiber, sheet, and normal to sheet directions.

References on Hill's equation of the active contraction. To validate the proposed finite element model, the end-diastole and end-systolic strains in the equatorial area were calculated and compared with the experimental measurements. It was found that the circumferential component of the strain, Ecc, shows a good agreement with the experimental data. Since Ecc is the main volumetric index of the LV cavity, it is reasonable to judge that the finite element model in this paper can well represent the pump functioning of the LV. On the other hand, the shear component of the end-diastolic strain, Ecr, is relevant to the rotation of the endocardial and epicardial surfaces about the longitudinal axis of LV. Since the deformation process of Ecr is subjected to few kinematical constraints, a large variability can be expected. As a consequence, the results obtained from different experiments in measuring Ecr differ from each other [3]. Up to now, there have been few computational models that have the ability to well estimate the shear strain, especially Ecr. The discrepancies between the calculated and the measured strains may be caused by several factors, such as the inadequacy of the myocardial constitutive model, neglecting the crossover of the myofibers in the myocardium, and simplifying the

[1] T. Arts, A model of the mechanics of the left ventricle, Ann. Biomed. Eng. 7 (1979) 299–318. [2] P.H.M. Bovendeerd, T. Arts, J.M. Huyghe, D.H. Campen, R.S. Reneman, Dependence of local left ventricular wall mechanics of myocardial fiber orientation: a model study, J. Biomech. 25 (1992) 1129–1140. [3] P.H.M. Bovendeerd, W. Kroon, T. Delhaas, Determinants of left ventricular shear strain, Am. J. Physiol. – Heart Circ. Physiol. 297 (2009) H1058–H1068. [4] K.D. Costa, J.W. Holmes, A.D. McCulloch, Modeling cardiac mechanical properties in three dimensions, Philos. Trans. R. Soc. A 359 (2001) 1233–1250. [5] F. Dorri, P.F. Niederer, P.P. Lunkenheimer, R.H. Anderso, The architecture of the left ventricular myocytes relative to left ventricular systolic function, Eur. J. Cardio-thorac. Surg. 37 (2010) 384–392. [6] T.S.E. Eriksson, A.J. Prassl, G. Plank, G.A. Holzapfel, Influence of myocardial fiber/sheet orientations on left ventricular mechanical contraction, Math. Mech. Solids 18 (2013) 592–606. [7] S. Göktepe, S.N.S. Acharya, J. Wong, E. Kuhl, Computational modeling of passive myocardium, Int. J. Numer. Methods Biomed. Eng. 27 (2011) 1–12. [8] Q. Gu, J.P. Conte, Z. Yang, A. Elgamal, Consistent tangent moduli for multiyield-surface J2 plasticity model, Comput. Mech. 48 (2011) 97–120. [9] J.M. Guccione, K.D. Costa, A.D. McCulloch, Finite element stress analysis of left ventricular mechanics in the beating dog heart, J. Biomech. 28 (1995) 1167–1177. [10] J.M. Guccione, L.K. Waldman, A.D. McCulloch, Mechanics of active contraction in cardiac muscle: Part II—Cylindrical models of the systolic left ventricle, J. Biomech. Eng. 115 (1993) 82–90. [11] J.C. Han, K. Tran, A.J. Taberner, D.P. Nickerson, R.S. Kirton, P.M.F. Nielsen, M. L. Ward, M.P. Nash, E.J. Crampin, D.S. Loiselle, Myocardial twitch duration and

J.J. Shen et al. / Computers in Biology and Medicine 75 (2016) 63–73

[12] [13] [14]

[15] [16] [17]

[18] [19]

[20] [21]

[22]

[23]

[24] [25] [26]

[27] [28]

[29] [30]

[31]

[32]

[33]

the dependence of oxygen consumption on pressure-volume area: experiments and modelling, J. Physiol. 590 (2012) 4603–4622. A. Hoger, Virtual configurations and constitutive equations for residually stressed bodies with material symmetry, J. Elast. 48 (1997) 125–144. G.A. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach For Engineering Science, John Wiley & Sons Ltd, New York, 2000. G.A. Holzapfel, R.W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization, Philos. Trans. R. Soc. A 367 (2009) 3445–3475. Z. Hu, D. Metaxas, L. Axel, In vivo strain and stress estimation of the heart left and right ventricles from MRI images, Med. Image Anal. 7 (2003) 435–444. R.M. Huisman, P. Elzinga, N. Westerhof, P. Sipkema, Measurement of left ventricular wall stress, Cardiovasc. Res. 14 (1980) 142–153. J.M. Huyghe, D.H. van Campen, T. Arts, R.M. Heethaar, The constitutive behaviour of passive heart muscle tissue: a quasi-linear viscoelastic formulation, J. Biomech. 24 (1991) 841–849. P.M. Janssen, W.C. Hunter, Force, not sarcomere length, correlates with prolongation of isosarcometric contraction, Am. J. Physiol. 269 (1995) H676–H685. C. Jiang, G. Liu, X. Han, Z. Zhang, W. Zeng, A smoothed finite element method for analysis of anisotropic large deformation of passive rabbit ventricles in diastole, Int. J. Numer. Methods Biomed. Eng. 31 (2015) 1–25. K. Sagawa, L. Maughan, H. Suga, K. Sunagawa, Cardiac Contraction and the Pressure-Volume Relationship, Oxford University Press, Oxford, 1998. I. LeGrice, Y. Takayama, J. Covell, Transverse shear along myocardial cleavage planes provides a mechanism for normal systolic wall thickening, Circ. Res. 77 (1) (1995) 182–193. I.J. LeGrice, B. Smaill, L. Chai, S. Edgar, J. Gavin, P.J. Hunter, Laminar structure of the heart: ventricular myocyte arrangement and connective tissue architecture in the dog, Am. J. Physiol. – Heart Circ. Physiol. 269 (2) (1995) H571–H582. C.E. Miller, M.A. Vanni, B.B. Keller, Characterization of passive embryonic myocardium by quasi-linear viscoelasticity theory, J. Biomech. 30 (1997) 985–988. P.M. Nash, P.J. Hunter, Computational mechanics of the heart, J. Elast. 61 (2000) 113–141. J.H. Omens, Y.C. Fung, Residual strain in rat left ventricle, Circ. Res. 66 (1990) 37–45. J.H. Omens, D.A. MacKenna, A.D. McCulloch, Measurement of strain and analysis of stress in resting rat left ventricular myocardium, J. Biomech. 26 (1993) 665–676. B. Oommen, M. Karamanoglu, S.J. Kovacs, Modeling time varying elastance: the meaning of load-independence, Cardiovasc. Eng.: Int. J. 3 (2003) 123–136. S. Pezzutoa, D. Ambrosia, A. Quarteronib, An orthotropic active-strain model for the myocardium mechanics and its numerical approximation. Eur. J. Mech. A/Solids 48 (2014) 83–96. A. Rachev, S. Greenwald, Residual strains in conduit arteries, J. Biomech. 36 (2003) 661–670. H. Schmid, M.P. Nash, A.A. Young, P.J. Hunter, Myocardial material parameter estimation—a comparative study for simple shear, J. Biomech. Eng. 128 (2006) 742–750. J.J. Shen, Research on Heart Contact Mechanics in the Minimally Invasive Surgery for Mitral Valve (Ph.D. thesis), Nanjing University of Aeronautics and Astronautics, 2014. J.J. Shen, A structurally based viscoelastic model for passive myocardium in finite deformation, Comput. Mech. (2016), http://dx.doi.org/10.1007/ s00466-016-1303-1, accepted for publication. J.J. Shen, M. Kalantari, J. Kövecses, J. Angeles, J. Dargahi, Viscoelastic modeling of the contact interaction between a tactile sensor and an atrial tissue, IEEE Trans. Biomed. Eng. 59 (2012) 1727–1738.

73

[34] J.J. Shen, C.G. Li, H.T. Wu, M. Kalantari, The frictional contact analysis between a tactile sensor and atrial tissue in viscoelasticity, Int. J. Appl. Mech. 5 (2013) 13500281–135002828. [35] J.J. Shen, F.Y. Xu, G.P. Jiang, A general three dimensional numerical technique for determining the contact area of an arbitrary punch on an elastic half-space, Int. J. Appl. Mech. 8 (2016) 16500051–165000518. [36] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer, New York, 2000. [37] S.R. Summerour, J.L. Emery, B. Fazeli, J.H. Omens, A.D. McCulloch, Residual strain in ischemic ventricular myocardium, J. Biomech. Eng. 120 (1998) 710–714. [38] H.E.D.J. ter Keurs, W.H. Rijnsburger, R. van Heuningen, M.J. Nagelsmit, Tension development and sarcomere length in rat cardiac trabeculae: evidence of length-dependent activation, Circ. Res. 46 (1980) 703–714. [39] M. Vendelin, P.H.M. Bovendeerd, T. Arts, J. Engelbrecht, D.H. van Campen, Cardiac mechanoenergetics replicated by cross-bridge model, Ann. Biomed. Eng. 28 (2000) 629–640. [40] F.J. Vetter, A.D. McCulloch, Three-dimensional stress and strain in passive rabbit left ventricle: a model study, Ann. Biomed. Eng. 28 (2000) 781–792. [41] L.K. Waldman, D. Nosan, F. Villarreal, J.W. Covell, Relation between transmural deformation and local myofiber direction in canine left ventricle, Circ. Res. 63 (3) (1988) 550–562. [42] H. Wang, X. Luo, H. Gao, R. Ogden, B. Griffith, C. Berry, T. Wang, A modified Holzapfel–Ogden law for a residually stressed finite strain model of the human left ventricle in diastole, Biomech. Model. Mechanobiol. 13 (2014) 99–113. [43] D. Westermann, M. Kasner, P. Steendijk, F. Spillmann, A. Riad, K. Weitmann, W. Hoffmann, W. Poller, M. Pauschinger, H.-P. Schultheiss, et al., Role of left ventricular stiffness in heart failure with normal ejection fraction, Circulation 117 (16) (2008) 2051–2060.

Jing Jin Shen received the Ph.D. degree in Mechanical Engineering from Nanjing University of Aeronautics and Astronautics, Nanjing, China, in 2014. He is currently an Assistant Professor of Automation at Nanjing University of Posts and Telecommunications, Nanjing, China. His current research interests include tactile sensing, tissue modeling, and contact analysis.

Feng Yu Xu received the Ph.D. degree in Mechanical Engineering from Southeast University, Nanjing, China, in 2009. He is currently an Associate Professor of Automation at Nanjing University of Posts and Telecommunications, Nanjing, China. He has published many papers and has six patents. Dr. Xu has been a Principal Reviewer of several major NSFC proposals in the area of mechanical engineering.

Wen An Yang received the Ph.D. degree in Mechanical Engineering from Nanjing University of Aeronautics and Astronautics, Nanjing, China, in 2012. He is currently an Assistant Professor of automation at Nanjing University of Aeronautics and Astronautics, Nanjing, China. His current research interests include robotics, mechanism optimization, and finite element analysis.