1. Introduction

The understanding of the physiological and pathophysiological behavior of the heart can be improved by computer aided simulations of the electrical and mechanical processes in the myocardium. These simulations are of importance e.g. for studying cardiac arrhythmias and the influence of medicaments as well as for the computer aided planning of surgical interventions. Different approaches allow the simulation of electromechanical processes in the myocardium. Commonly, the approaches are restricted to a reconstruction of either electrophysiological or mechanical processes. Hereby, the influence of feedback mechanisms is neglected. The approach presented in this work consists of combining electrophysiological models of single myocardial cells with models of the electrical current flow through the myocardium. The electrophysiological models are coupled with models of the cellular force development. The calculated forces are used by an elastomechanical model of the myocardium. The resulting deformation is delivered to all models e.g. to consider feedback mechanisms. The approach is used to generate a three dimensional model of an area of the ventricular myocardium. The model allows the simulation of the electrical excitation propagation, force development and deformation. The model is derived from a detailed electrophysiological model of ventricular cells [Noble et al., 1998]. The electrical coupling of cells is achieved using of an extension of the bidomain model, which consists of calculation of the intracellular and extracellular current flow [Sachse et al., 2000]. Therefore, Poisson's equation for electrical current fields is applied on grids, which are geometrically deformed by the stretch. The force model described in [Rice et al., 1999] uses parameters delivered by the electrophysiological model. The elastomechanical model is generated from a strain energy function proposed by [Hunter et al., 1998] and numerical methods of continuum mechanics [Bathe, 1996]. The calculation takes the orientation and lamination of myocytes into account, which leads to an anisotropy of the electrical conductivity and the elastomechanical behavior. 

2. Methods 

2.1 Microscopic Anatomy of the Human Heart 

2.1.1 Overview 

The myocard consists primarily of discrete myocytes, which are arranged in an oriented and laminated structure [Streeter and Bassett, 1966,Streeter, 1979,LeGrice et al., 1995]. The myocytes are surrounded by the extracellular space, which contains e.g. collagen fibers, water, and electrolytes. Additionally, the myocard is pervaded by nerves, capillaries, blood vessels, and lymphatic vessels. It can be distinguished between myocytes assigned to the working myocard and myocytes with the task of initiating and transmitting electrical excitation. This work is focussed upon the myocytes of the working myocard. 

2.1.2 Myocytes of the Working Myocard 

Myocytes are enclosed by the sarcolemma, a phospholipid bilayer, which delimits the extracellular from the intracellular space. In the intracellular space resides the nucleus, the mitochondria, and the myofibrils. Myofibrils are tube shaped contractile elements taking a high part of the cellular volume. They have a thickness of approximately 1 μm and are divided every approximately 2.5 μm by the Z disk into the sarcomeres. The sarcomere contains the myofilaments, which are of importance for the mechanical contraction. The myofilaments are composed of actin (thin) and myosin (thick) filaments [Bers, 1991]. The thin filaments lead from the Z disks approximately 1 μm towards the center of the sarcomere. The filaments consist of long, flexible, coiled coil molecules, i.e. tropomyosin, at which the actin and troponin is attached. At every seventh actin a troponin is bound to the tropomyosin. Troponin consists of three components: troponin T, troponin I and troponin C. Troponin T is connected to the tropomyosin and the two further components. Troponin I inhibits interactions between actin and myosin. Troponin C is of importance for the initiation of the mechanical contraction by binding of calcium. The thick filaments have a length of approximately 1.6 μm and are centered in the sarcomere. They are parallel to and arranged between the thin filaments. The thick filament is composed primarily of the myosin molecule, which has a length of approximately 160 nm. Myosin in myocytes of cardiac, skeletal and smooth muscles is subdivided into two heads and tails. The head shows an ATP and an actin binding site. The tails are coiled coil molecules and form the main axis of the thick filaments. The sarcolemma intrudes at the adjacencies of the Z disks to form the transversal tubuli and to infold the myofibrils. The system of longitudinal tubuli plays in contrast to the skeletal myocytes in cardiomyocytes only a subsidiary role [Bers, 1991]. Myocytes are of irregular shape, but a dominant principal axis can be assigned. This leads to a macroscopic anisotropic electrical and elastomechanical behavior. The ratios of volume and area of the cell components are differing for atrial and ventricular myocytes [Bers, 1991]. 

2.1.3 Gap Junctions 

The intracellular space of myocytes is coupled by gap junctions (nexus), which are located by bundle at the intercalated disks [Dhein, 1998,Dellmar et al., 1999]. Intercalated disks are disk shaped segments, which mechanically couple cells [Forbes and Sperelakis, 1985]. Primarily, the disks can be found at or near to the ends of myocytes [Saffitz and Yamada, 1999]. A gap junction is cylinder or barrel shaped with a diameter of 1.5 - 2.0 nm and a length of approximately 2 - 12 nm. A molecule of atomic weight up to 1 kD can pass through the gap junction [Jongsma and Rook, 1999,Yeager, 1999], e.g. nutrients, metabolites, and ions. A gap junction is built up by two hemi-channels, so-called connexons, piercing the sarcolemma of the involved cells. The connexons are formed by six integral membrane proteins, so-called connexins. More than one dozen of connexins have been cloned, which are named by their atomic weight ranging from 25 to 50 kD. The most abundant connexin in the mammalian myocard is connexin43, which is also expressed in ovary, uterus, kidney, and lens epithelium. The connexin32 is found in liver, stomach, kidney, and brain, but not in any part of the heart. In the myocard each myocyte is coupled by gap junctions non-uniformly with other myocytes, e.g. in canine a myocytes is coupled with 9.1 ± 2.2 myocytes [Hoyt et al., 1989]. A distinction can be made between longitudinal and transversal gap junctions. A longitudinal gap junction is oriented in approximately the same direction as the first principal axis of adjacent myocytes, a transversal gap junction is oriented perpendicular thereto. The density and distribution of orientations of gap junctions differs depending on the tissue, e.g. the density in the sinus and atrioventricular node is smaller than in the ventricular myocard. The average density of gap junctions in longitudinal orientation is larger than the transversal orientation [Hoyt et al., 1989]. The average length of longitudinal gap junctions is smaller than the length of transversal gap junctions. Both circumstances lead to a macroscopic anisotropic electrical conductivity. 

2.1.4 Cardiac Collagen Network 

The myocard is pervaded and surrounded by a mesh of extracellular collagen fibers, which are composed of a multitude of collagen fibrils synthesized by the cardiac fibrocytes. Collagen takes 2 - 5 % of the weight of the heart. Furthermore, fibers of elastin draw through the myocard. The content and the structure of the connective tissue are dependent on age, pathologies and species. Both, collagen and elastin are fibrous proteins of the extracellular matrix, which is a determinant for the viscoelastic behavior of the myocard. The network serves for the mechanical coupling of the myocytes, capillaries, and lymphatic vessels. Collagen fibrils have a thickness from 30 - 70 nm [Caulfield and Borg, 1979]. The fiber thickness is in physiologic cases between 120 and 150 nm. An increase up to 250 - 300 nm is possible in pathophysiologic cases, e.g. hypertrophy, hypertension and myocardial infarction [Abrahams et al., 1987,Ju and Dixon, 1996]. The density of collagen fibers is depending on the tissue. E.g. a small density can be found in papillary muscle and trabeculae, a high density in the left ventricular myocard, the endocardium and epicardium [Caulfield and Borg, 1979,Weber et al., 1994].

2.2 Modeling of Cellular Electrophysiology 

Table 1: Electrophysiological models of cardiac cells.
 

Cell Type	Species	Reference
Purkinje fiber	-	[Noble, 1962]
Purkinje fiber	-	[McAllister et al., 1975]
Ventricular myocard	mammal	[Beeler and Reuter, 1977]
Purkinje fiber	mammal	[DiFrancesco and Noble, 1985]
Atrial myocard	rabbit	[Hilgemann and Noble, 1987]
Atrial myocard	rabbit	[Earm and Noble, 1990]
Ventricular myocard	mammal	[Luo and Rudy, 1991]
Ventricular myocard	guinea-pig	[Luo and Rudy, 1994a,Luo and Rudy, 1994b]
Sinus node	mammal	[Demir et al., 1994]
Ventricular myocard	canine	[Demir et al., 1996]
Atrial myocard	human	[Courtemanche et al., 1998]
Ventricular myocard	guinea-pig	[Jafri et al., 1998]
Ventricular myocard	guinea-pig	[Noble et al., 1998]
Atrial myocard	human	[Nygren et al., 1998]
Ventricular myocard	human	[Priebe and Beuckelmann, 1998]
Ventricular myocard	canine	[Winslow et al., 1999,O'Rourke et al., 1999]

2.2.1 Overview 

The electrophysiological state of cells can be described by a spatial distribution of ion concentrations, which are changed by passive and active transport mechanisms. The transport of ions is time dependent as well as influenced by gradients of concentrations and the electrical field. Most electrophysiological models of cells are constructed from the classical work of Hodgkin and Huxley, who described quantitatively the active and passive behavior of a neuron membrane [Hodgkin and Huxley, 1952]. Hodgkin and Huxley used voltage clamp techniques to measure currents through the membrane of giant axons of squids. From this data an equivalent circuit consisting of resistors, a capacity, and voltage sources was parameterized. Partly, the resistors were nonlinear time and voltage dependent. In the last years a large number of models of cardiac cells were constructed (table 1), with increasing abilities to describe the different electrophysiological mechanisms. Primarily, the models are produced from animal experiments. Modern models include detailed descriptions of the behavior of intracellular structures as well as of the influence of drugs and deformation. 

2.2.2 Noble-Varghese-Kohl-Noble Model 

A foundation of this work is the Noble-Varghese-Kohl-Noble model of a ventricular cell [Noble et al., 1998]. Thereby, the time derivative of the transmembrane voltage V_m is described by: 

(1)

with the membrane capacity C_m and the summary current through the membrane I_sum

with the following parameters: 

background K current
time-independent K current
time-dependent, delayed K current	, ,
ACh-dependent K current
background Na current
fast Na current
voltage dependent Na current
background Ca current
L-type Ca current	,  ,  ,
	,  ,
Na-K exchange pump current
Na-Ca exchange pump current	,
stretch activated current	,  ,  ,
	,

Figure 1: Simulations with the Noble-Varghese-Kohl-Noble model. The transmembrane voltage V_m is dependent on the stimulus frequency. For each frequency a single action potential is visualized [Sachse, 2001].

Figure 2: Simulations with the Noble-Varghese-Kohl-Noble model. The intracellular calcium concentration [Ca²⁺]_i is dependent on the stimulus frequency. For each frequency a single course of the calcium concentration is visualized [Sachse, 2001].

The figures 1 and 2 shows the influence of stimulus frequency to the course of the transmembrane voltage V_m and intracellular calcium concentration [Ca²⁺]_i. Thereby, stretch activated currents are neglected. With higher stimulus frequency the resting voltage is increased and the duration of the action voltage is decreased.

2.2.3 Intracellular Mechano-Electrical Feedback 

The Noble-Varghese-Kohl-Noble model includes dependencies of electrophysiological parameters on the length or tension of the sarcomere. The mechano-electrical feedback is realized by introducing selective and non selective stretch-activated ion conductances, a length/tension dependent modulation of calcium binding to troponin and a length/tension dependent modulated sarcoplasmatic leak current. A modification of the model is performed, whereby an adaption according to measurements published in [White et al., 1993] is obtained [Sachse et al., 2000]. This modification concerns the stretch dependent action potential duration. 

Figure 3: (a) Transmembrane voltage, (b) calcium concentration in the cytoplasm, (c) concentration of calcium bound to troponin C, and (d) concentration of calcium in the release part of the sarcoplasmic reticulum dependent on length of sarcomere calculated with Noble-Varghese-Kohl-Noble model (from [Sachse, 2001]). The cell is excited by applying a stimulus current at t = 25 ms with a length of 3 ms. The sarcomere length ranges from 2.0 to 2.2 μm. The default length of the sarcomere is 2.2 μm. The stimulus frequency was set to 1 Hz.

Figure 4: Initiation of action impulse by stretch of sarcomere. The cell is excited by applying a stimulus current at t = 1 s. At t = 2 s a mechanical stretch of 5 ms was performed delivering a sarcomere length SL from 2 to 2.9 μm. The default length of the sarcomere is 2 μm (from [Sachse et al., 2000]).

The influence of stretch on the run of the transmembrane potential is illustrated in figure 3. Thereby, the stretch is specified by the length of the sarcomere with a default of 2 μm. The resting potential as well as the progression of the action potential are dependent on the length of the sarcomere. The initiation of an excitation by mechanical stretch is depicted in figure 4. The stretch is applied for a duration of 5 ms with varying strength. Depending on the strength of stretch an effect ranging from a small change of the resting potential to an excitation of the cell can be achieved.

2.3 Modeling of Excitation Propagation 

2.3.1 Approaches 

Figure 5: Modeling of electrical intercellular coupling. The myocytes are coupled via gap junctions and through the extracellular space.

A propagation of electrical excitation from one cell to neighboring cells is primarily achieved by intercellular transport of ions via the gap junctions. Also extracellular potentials resulting from the electrical activity of cells or from an external current flow can modulate the propagation and initiate an excitation (figure 5). Two different classes of approximations of the excitation propagation in the myocard can be distinguished: microscopic and macroscopic approaches. The macroscopic approach allows the combining of groups of cells and their common treatment. In contrast, microscopic models at a cellular level split cells in components, which are separately treated. In the last years different approaches for the macroscopic excitation propagation were developed: 

• Cellular automatons. Rules are included defining the time delay and the neighborhood for the propagation [Eifler and Plonsey, 1975,Killmann et al., 1991,Saxberg and Cohen, 1991,Wei et al., 1995,Siregar et al., 1996,Werner et al., 1998,Siregar et al., 1998]. 
• Excitable dynamics equations or reaction diffusion systems 

[FitzHugh, 1961,Rogers and McCulloch, 1994,Panfilov, 1999]. 
• Resistor networks/monodomain models. These models incorporate the effect of coupling the intracellular space with gap junctions [Rudy and Quan, 1989,Virag et al., 1999]. 
• Bidomain models. Bidomain models are an extension of monodomain models including the effects of the extracellular space [Henriquez and Plonsey, 1989,Sepulveda and Wikswo, 1994,Henriquez et al., 1996]. 

All these models allow the inclusion of anisotropic effects resulting from the orientation of myocytes, e.g. by using conductivity tensors. Microscopic models deliver information concerning the stochastic behavior of the myocard [Spach and Heidlage, 1995,Spach et al., 1999]. Anisotropic effects are implicitly included by the cell geometry as well as by the distribution and orientation of gap junctions. 

2.3.2 Bidomain Model 

The bidomain model treats the electrical behavior of tissue in two domains, in the intracellular and extracellular space, which are separated by the cell membrane. In each domain Poisson's equation for fields of stationary electrical current is fulfilled: 

			(2)
			(3)

with the intracellular potential Φ_i , the extracellular potential Φ_e, the intracellular conductivity tensor σ_i, the extracellular conductivity tensor σ_e, the intracellular current source density I_si, the extracellular current source density I_se, and the surface to volume ratio β of cells. The intracellular conductivity σ_i consists of conductivities for the intracellular components and for the gap junctions. The domains are coupled by the current density I_m through the cell membrane. 

Figure 6: Bidomain modeling of cardiac electrophysiology.

The following method can be chosen to couple the bidomain equations with the electrophysiological cell models (figure 6) [Hooke et al., 1992]. The method bases on the iterative solving of Poisson's equation with numerical techniques: 

• Potentials are determined from current source densities. 
• Current sources are calculated from potentials. 

Therefore, commonly the finite-difference or finite-element method is applied [Press et al., 1992,Bathe, 1982]. In a first step the current source density I_im delivered by the transmembrane potential V_m = Φ_i - Φ_e is determined: 

(4)

In a second step the extracellular potential Φ_e is calculated from the current source density I_im:

(5)

The calculation of Φ_e is commonly numerically expensive, because the solving of a large system of linear equations is necessary. In a third step the intracellular source density I_si is determined and delivered to the electrophysiological cell model: 

(6)

The left side of this equation describes a current source density delivered by the intracellular potentials: 

(7)

2.3.3 Extension of the Mono- and Bidomain Model 

Figure 7: Coupling of myocytes with gap junctions and through the extracellular space. The deformation of a region changes the intra- and extracellular conductivity. The resistor yielded by the gap junction is not changed. 

In a previous paper an extension of the mono- and bidomain model was introduced, which allows to take the deformation of tissue into account (figure 7) [Sachse et al., 2000]. This extension delivers conductivity tensors σ_i and σ_e for the intra- and extracellular space respectively, which follow the rules of model assumptions. In principal the method consists of extracting the stretch of regions resulting from an arbitrary deformation. The extracted stretch is used to construct a conductivity tensor in a local coordinate system. Different weights allow to choose a specific behavior of the conductivity. The local conductivity tensor is transformed into the global coordinate system.

2.4 Modeling of Cellular Force Development 

2.4.1 Overview 

The development of force in the contractile elements of myocytes is provoked by an increase of the concentration of intracellular calcium [Ca]_i. The progression of the force is modulated by the progression of the concentration [Ca]_i. Commonly, the increase of the concentration [Ca]_i is a result of an electrical excitation. The progression of the electrical excitation influences the progression of the force development (electro-mechanical feedback). Therefore, many models of cellular force development use the concentration [Ca]_i to define rate coefficients, which depict the interaction between states [Landesberg and Sideman, 1994a,Landesberg and Sideman, 1994b,Rice et al., 1999,Rice et al., 2000]. The states describe e.g. the binding of intracellular Ca²⁺ to the troponin complex and the cross-bridge cycling. Further parameters influencing the rate coefficients are the sarcomere length and the state variables. 

2.4.2 Rice-Winslow-Hunter Model 

Table 2: Tropomyosin and cross bridge states of
Rice-Winslow-Hunter Model 3 of cardiac cells.

state	Tropomyosin	No. of cross bridges
N0	non permissive	0
N1	non permissive	1
P0	permissive	0
P1	permissive	1

 

Table 3: Ca²⁺ binding states of Rice-Winslow-Hunter
Model 3 of cardiac cells.

state	Ca2+binding to Troponin C
T	no
TCa	yes

A foundation of this work are the Rice-Winslow-Hunter models of cardiac muscle [Rice et al., 1999]. As an example of the modeling a short description of the 3rd model is given. This model consists of 6 states, N0, N1, P0, P1, T, and TCa (tables 2 and 3) with: 

			(8)
			(9)

The interaction between the states of the model is described by a system of 1st order differential equations: 

(10)

with the 6 × 6 matrix R consisting of rate coefficients. Partly, the rate coefficients are dependent on the sarcomere length SL and the concentration [Ca]_i. The normalized force F is determined by 

(11)

with the sarcomere overlap function  α = α(SL) and the maximal force F_max. The states P1 and N1 are the force generating states.

2.5 Modeling of Elastomechanical Behavior

2.5.1 Principle of Virtual Displacements 

The equilibrium of a body is achieved if the internal and external forces are balanced [Bathe, 1996]. The equilibrium at time  t + Δt can be expressed using the principle of virtual displacements: 

(12)

with the volume   ^{t + Δt}V , the components of the Cauchy stress tensor   ^{t + Δt}τ_ij , the strain tensor  δ_t + Δt e_ij , and the external virtual work R. The formula uses the summation convention of Einstein, where repeated subscripts become the designation for summation. The strain tensor is defined as 

(13)

with the components of the virtual displacement vector δu_i. The external virtual work R is sub-dived in applied force densities ^{t + Δt} f_i^B   and surface tensions ^{t + Δt} f_i^S  

(14)

with the surface  ^{t + Δt} S _f .

2.5.2 Strain Energy Density Function 

The strain energy density function W proposed by [Hunter et al., 1997,Hunter et al., 1998] takes the anisotropic and inhomogeneous behavior of the myocard into account. Three micro-structural, orthogonal axes are distinguished: the fiber, sheet and sheet normal axis. For each axis i a set of parameters, k_i, a_i, and β_i, describes its contribution to the strain energy density, called pole-zero law: 

(15)

with the diagonal components of the Green-Lagrange strain tensor e_ii. The parameter k_i is set to zero, if e_ii is negative. The strain energy density function W is defined for |e_ii| < a_i . The function shows large values for e_ii approaching a_i , reflecting the steep rise in tension coming upon a strain limit. The strain energy was extended by terms representing the incompressibility of the myocard. The energy does not comprise shear and viscoelastic effects. The parameterization of the function W was performed by uniaxial measurements of canine ventricle in the different directions of the axes. Hereby, the parameters of different regions in the myocard were collected.

3. Results 

The developed numerical model has the purpose to achieve knowledge concerning the cardiac deformation and its influence to the initiation and propagation of electrical excitation and to the force development. The model combines and extends the presented cellular and macroscopic models. It consists of 

• a single cell electrophysiological model with stretch dependent behavior 
• an extended bidomain model taking stretch into account 
• a single cell model of the force development with inclusion of stretch effects 
• an elastomechanical model 

3.1 Modeling of the Cardiac Electromechanics 

Figure 8: Modeling of cardiac electromechanics.

The interdependencies of the different data are depicted in figure 8. As an electrophysiological model the modified Noble-Varghese-Kohl-Noble model was used [Noble et al., 1998], whereby stretch dependent ion channels were included [Sachse et al., 2001]. The intercellular electrical coupling through the gap junctions and extracellular space was performed with the extended bidomain model [Sachse et al., 2000]. The engaged force model was the Rice-Winslow-Hunter model (type 3). The elastomechanical behavior was modeled by numerical methods of continuum mechanics [Bathe, 1996] using the strain energy function proposed by [Hunter et al., 1997]. 

Figure 9: Model of myocardial area. The electrical excitation is initiated by applying current at plane z=0. The area is mechanically fixed at point P. The orientation of myocytes is indicated by f, the sheet orientation by s, and the sheet normal by n.

3.2 Simulations 

Simulations with the electromechanical model of a myocardial area were performed. The results were visualized with surface based techniques (figures 10 and 11). The applied anatomical model consisted of 20 × 20 × 20 cubic voxels, each with a size of 0.1 mm × 0.1 mm × 0.1 mm (figure 9). The activation starts at time 0 ms by application of a sufficient electrical current at plane z = 0. The principal axis of myocytes f was chosen parallel to the z-axis. The lamination of the myocytes is determined by the sheet orientation s and the sheet normal n. Transversal isotropy of the electrical conductivities was set. Anisotropy of the elastomechanical parameters and incompressibility was assumed. The central position of the plane z=0 was fixed, i. e. the displacements were set to zero. The electrophysiological modeling (equation 1 and 2) and the force development modeling (equation 11) was performed with the Euler method using a time step of 20 μs [Press et al., 1992]. The bidomain model  used a Gauss-Seidel iteration every 20 μs [Press et al., 1992]. The deformation was calculated with a time step of 1 ms. The system of linear equations resulting from equation 13 was solved by the conjugate gradient method. 

Figure 10: Transmembrane voltage at time (a) 0 ms, (b) 3 ms, (c) 5 ms, and (d) 8 ms in an anisotropic model of myocardial area. The model consists of 20 × 20 × 20 cubic voxels with a size of 0.1 mm × 0.1 mm .

The simulations with the combined model show processes of different time scale. The process of excitation propagation is rapidly spreading over the myocardium (figure 10). The force development and the resulting deformation follows with a significant delay (figure 11).

Figure 11: Color-coded normalized force and deformation at time (a) 0 ms , (b) 50 ms, (c) 100 ms, (d) 150 ms, (e) 200 ms, and (f) 250 ms in an anisotropic model of myocardial area. The model consists of 20 × 20 × 20 cubic voxels with a size of 0.1 mm × 0.1 mm × 0.1 mm. The central position of the plane z=0 was fixed, i. e. the displacements were set to zero. The wire frame shows the reference configuration.

4. Discussion and Conclusions 

The presented model describes aspects of the electromechanical behavior of a myocardial region. The model combines an electrophysiological, an excitation propagation, a force development and an elastomechanical model. The performed simulations illustrate effects of myocardial electromechanical behavior. These effects are of great significance for the development of realistic models of the whole heart. The presented combination of electrophysiological and mechanical modeling allows the simulation of therapeutical interactions with pharmaceutical, electrical and mechanical methods. This functionality is of importance e.g. for studying of cardiac arrhythmias and for the computer aided planning of surgical interventions. Further work will be focussed on the inclusion of inertia and surface forces by blood pressure as well as modeling of viscoelastic tissue properties. Of interest is also the integration of models of the metabolism. 

Acknowledgments

The authors want to thank Dr. P. Kohl, University Laboratory of Physiology, Oxford, UK, for his help to parameterize the electrophysiologic model. R. Mayer, Rechenzentrum, Univerität Karlsruhe (TH), supported our work by providing the necessary visualization and computing resources.