IJBEM, Vol. 2, No. 2, pp. 159-182, 2000.

Institute of Biomedical Engineering, University of Karlsruhe, 76128 Karlsruhe, Germany

Abstract. A model of the electromechanical behavior of a myocardial region is presented. The model combines an electrophysiological, a force development and an excitation propagation model. All of these models incorporate the effects of deformation of the myocardium. An extension of the traditional bidomain model for excitation propagation is proposed. The extension describes the stretch dependency of the conductivity tensor of the intra- and extracellular space and is constructed outgoing from physically motivated assumptions, which simplify the behavior of the conductivity tensor. The extension makes usage of the deformation gradient tensor, which is a foundation in the theory of continuums mechanics. The performed simulations illustrate some effects of myocardial electromechanical behavior.

Keywords: Cardiac Modeling, Mechano-Electrical Feedback, Electro-Mechanical Feedback, Cellular Models, Electrophysiology, Excitation-Propagation, Bidomain Model

Introduction

Knowledge concerning the development and propagation of electrical excitation in the myocardium is of importance for the understanding of the physiological and pathophysiological behavior of the heart. Different approaches allow the simulation of an excitation development and propagation in the myocardium. An approach consists of combining detailed electrophysiological models of single myocardial cells with models of the electrical conductivity of the gap junctions and extracellular space.

Commonly, the electrophysiological models describe the concentration and flow of ions as well as the resistor of cellular structures and the transmembrane potential by a set of coupled differential equations. Some of these models are capable of reproducing effects of mechanical stimuli, eg stretching the cell can lead to an excitation of a cell and influences the propagation of excitation. Stretch of the myocardium can occur eg by an enhanced blood inflow leading to an increase of the end-diastolic volume as well as a mechanical procedure from outside of the body and inside of the cavities. Under the influence of stretch arrhythmias can be initiated and sustained.

In this work a three dimensional model of a myocardial area is presented. The model allows the simulation of an electrical excitation propagation and of force development under the influence of stretch. The model is derived from a modified Noble-Kohl-Varghese-Noble model of a ventricular cell. The electrical coupling of cells is achieved by usage of an extension of the bidomain model, which consists of calculation of the intracellular and extracellular current flow. Therefore, Poisson's equation for electrical current fields is applied on grids, which are geometrically deformed by the stretch. A new approach is introduced to model the change of the electrical properties of the intra- and extracellular space depending on the deformation. The calculation takes into account the orientation of myocytes, which leads to an anisotropy of the electrical conductivity.

The method is applied to determine the stretch induced changes of the conduction velocity and to simulate the excitation initiation by stretch. Finally, a realistic model of region of the ventricular heart wall is presented implying the anisotropy introduced by the orientation of the myocytes. Demonstrated and visualized is the propagation of excitation outgoing from endocardial Purkinje fibers, the concentration of intracellular calcium and the generated force.

Methods

Microscopic Anatomy of the Human Heart

Overview

The myocardium consists primarily of discrete myocytes, which are arranged in an oriented and laminated structure [1][2][3]. The myocytes are surrounded by the extracellular space, which contains eg collagen fibers, water, and electrolytes. Additionally, the myocardium is pervaded by nerves, capillaries, blood vessels, and lymphatic vessels. It can be distinguished between myocytes assigned to the working myocardium and myocytes with the task of initiating and transmitting an electrical excitation. This work is focussed upon the myocytes of the working myocardium.

Myocytes of the Working Myocardium

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 ≈ 1 μm and are divided every ≈ 2.5 μm by the Z disk into the sarcomeres. The sarcomere contains the actin and myosin filaments, which are of importance for the mechanical contraction. 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 [4].

Myocytes are of irregular shape, but a dominant principal axis can be assigned. This leads to a macroscopic anisotropic electrical and elastomechanic behavior.

The ratios of volume and area of the cell components are differing for atrial and ventricular myocytes [4].

Gap Junctions

The intracellular space of myocytes is coupled by gap junctions, which are located by bundle at the intercalated disks. Intercalated disks are disk shaped segments, which mechanically couple cells. Primarily, the disks can be found at or near to the ends of myocytes [5].

It is distinguished 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. Eg the density in the sinus and atrioventricular node is smaller than in the ventricular myocardium. The average density of gap junctions in longitudinal orientation is larger than in the transversal orientation [6]. 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.

Cardiac Collagen Network

The myocardium is pervaded and surrounded by a mesh of extracellular collagen fibers, which are composed of a multitude of collagen fibrils. The collagen network is a determinant for the viscoelastic behavior of the myocardium. The network serves for the mechanical coupling of the myocytes, capillaries, and lymphatic vessels.

Collagen fibrils have a thickness from 30-70 nm [7]. The fiber thickness is in physiologic cases between 120 and 150 nm . An increase up to 250-300 is possible in pathophysiologic cases (hypertrophy) [8].

The density of collagen fibers is depending on the tissue. Eg a small density can be found in papillary muscle and trabeculae, a high density in the subendocardial and subepicardial left ventricular myocardium [7][9].

Modeling of Myocytes

Modeling of Electrophysiology

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 outgoing from the classical work of Hodgkin and Huxley, who described quantitatively the active and passive behavior of a neuron membrane [10]. 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 (see Table 1), with increasing abilities to describe the different electrophysiological mechanisms. Primarily, the models are produced outgoing from animal experiments. Modern models include detailed descriptions of the behavior of intracellular structures as well as of the influence of pharmaceutical and deformation.

TABLE 1. Electrophysiological models of cardiac cells.

Date	Publisher	Cell Type	Animal	Reference
1977	Beeler, Reuter	ventricular myocardium	mammal	[11]
1985	DiFrancesco, Noble	Purkinje fiber	mammal	[12]
1990	Earm, Noble	atrial myocardium	rabbit	[13]
1991	Luo, Rudy	ventricular myocardium	guinea pig	[14]
1994	Luo, Rudy	ventricular myocardium	guinea pig	[15][16]
1994	Demir, Clark, Murphey, Giles	sinus node	mammal	[17]
1996	Demir, O'Rourke, Tomaselli, Marban, Winslow	ventricular myocardium	canine	[18]
1998	Nygren, Fiset, Firek, Clark, Lindblad et al.	atrial myocardium	human	[19]
1998	Noble, Varghese, Kohl, Noble	ventricular myocardium	guinea pig	[20]
1998	Jafri, Rice, Winslow	ventricular myocardium	guinea pig	[21]
1999	Winslow, Rice, Jafri. Marban, O'Rourke	ventricular myocardium	canine	[22][23]

Noble-Varghese-Kohl-Noble Model

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

δV_m
δt

–

1
C_m

I_sum

I_b,K + I_K₁ + I_Kr₁ + I_Kr₂ + I_Ks + I_K,ACh +
I_b,Na + I_Na + I_p,Na +
I_b,Ca + I_Ca,L,K + I_Ca,L,Na + I_Ca,L,Ca +
I_Ca,L,K,ds + I_Ca,L,Na,ds + I_Ca,L,Ca,ds +
I_Na,K + I_Na,Ca + I_Na,Ca,ds +
I_Na-stretch + I_K-stretch + I_Ca-stretch +
I_Ns-stretch + I_An-stretch

with the following parameters:

membrane capacity C_m

background K current I_b,K

time-independent
K current I_K1

time-dependent,
delayed K current I_K_r1, I_K_r2, I_K_rs

ACh-dependent
K current I_K,ACh

background Na current I_b,Na

fast Na current I_Na

voltage dependent
Na current I_p,Na

background Ca current I_b,Ca

L-type Ca current I_Ca,L,Ca , I_Ca,L,Ca,ds , I_Ca,L,K ,
I_Ca,L,K,ds , I_Ca,L,Na , I_Ca,L,Na,ds

Na-K exchange
pump current I_NaK

Na-Ca exchange
pump current I_NaCa, I_NaCa,ds

stretch activated current I_Na-stretch , I_K-stretch , I_Ca-stretch , I_Ns-stretch , I_An-stretch

Intracellular Mechano-Electric Feedback

The Noble-Varghese-Kohl-Noble model includes dependencies of electrophysiological parameters on the length or tension of the sarcomere. The mechano-electric feedback is realized by introducing

selective and non selective stretch-activated ion conductances
a length/tension dependent modulation of calcium binding to troponin
a length/tension dependent modulated sarcoplasmatic leak current

A modification of the model is performed, whereby an adaption according to measurements [24] is obtained. This modification concerns the stretch dependent action potential duration. The influence of stretch on the run of the transmembrane potential is illustrated in figure 1. Hereby, 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.

$\begin{figure}\center\epsfig{file=ActionPotentialSarcomereLength.eps,width=.7\textwidth}\end{figure}$

Figure 1. Transmembrane potential dependent on length of sarcomere. The cell is excited by applying a stimulus current at t = 1 s. The sarcomere length SL ranges from 1.5 to 2.4 μm. The default length of a sarcomere is 2 μm

The initiation of an excitation by mechanical stretch is depicted in figure 2. 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.

$\begin{figure}\center\epsfig{file=ImpulseSarcomereLength.eps,width=.9\textwidth}\end{figure}$

Figure 2. 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,

Also the model includes electro-mechanic feedback, which is discussed in the following section, and the influence of neurotransmitters (ACh).

Modeling of Force Development

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-mechanic feedback).

Therefore, many models of cellular force development use the concentration [Ca]_i to define rate coefficients, which depict the interaction between states [25] [26] [27][28]. The states describe eg the bounding 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.

TABLE 2. States of force models of cardiac cells.

state	Tropomyosin	No. of cross bridges
N0	non permissive	0
N1	non permissive	1
P0	permissive	0
P1	permissive	1
state	Ca²⁺ bounding to Troponin
T	no
TC_a	yes

Rice-Winslow-Hunter Model

A foundation of this work are the Rice-Winslow-Hunter models of cardiac muscle [27]. 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 with N0 + N1 + P0 + P1 = 0 and T + TCa = 0 (see Table 2). The interaction between the states of the model is described by a system of 1st order differential equations:

δ
δt

N0
N1
P0
P1
T
TC_a

= R

N0
N1
P0
P1
T
TC_a

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

F =

α(P1 + N1)

F_max

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

Modeling of Excitation Propagation

Approaches

$\begin{figure}\center\epsfig{file=MyoGap.eps,width=.8\textwidth}\end{figure}$

Figure 3. 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 3).

Two different classes of approximations of the excitation propagation in the myocardium 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 [29] [30] [31] [32] [33] [34] [35].
Excitable dynamics equations or reaction diffusion systems [36][37][38].
Resistor networks/monodomain models. These models incorporate the effect of coupling the intracellular space with gap junctions [39][40].
Bidomain models. Bidomain models are an extension of monodomain models including the effects of the extracellular space [41][42][43].

All these models allow the inclusion of anisotropic effects resulting from the orientation of myocytes, eg by using conductivity tensors (see Appendix Electrical Conductivity Tensor).

Microscopic models deliver information concerning the stochastic behavior of the myocardium [44][45]. Anisotropic effects are implicitly included by the cell geometry as well as by the distribution and orientation of gap junctions.

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:

∇(σ_i∇ϕ_i) = βI_m - I_si (1)

∇(σ_e∇ϕ_e) = βI_m - I_se (2)

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.

$\begin{figure}\center\epsfig{file=bidomain.eps,width=1\textwidth}\end{figure}$

Figure 4. Bidomain modeling of cardiac electrophysiology.

The following method can be chosen to couple the bidomain equations with the electrophysiological cell models (see Figure 4) [46]. The method bases on the iterative solving of Poisson's equation with numerical techniques:

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

Therefore, commonly the finite-difference or finite-element method is applied [47][48]. In a first step the current source density I_im delivered by the transmembrane potential V_m = ϕ_i - ϕ_e is determined:

∇(σ_i∇V_m) = I_im

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

∇((σ_i + σ_e) ∇ϕ_e) = I_im

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:

∇(σ_i∇V_m) + ∇(σ_i∇ϕ_e) = β (C_m

δV_m
δt

– I_ion) – I_si

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

∇(σ_i∇V_m) + ∇(σ_i∇ϕ_e) = ∇(σ_i∇(V_m + ϕ_e)) = ∇(σ_i∇ϕ_i)

Modeling of Mechano-Electrical Feedback

Figure 5. Modeling of Mechano-Electrical Feedback.

Overview

Knowledge concerning the influence of stretch and tension to the initiation and propagation of electrical excitation, and to the force development in the heart can be achieved by development of a model, which combines and extends the presented cellular electrophysiological and propagation model. The model consists of

a single cell electrophysiological model with stretch dependent behavior
a single cell model of the force development with inclusion of stretch effects
an extended bidomain model taking stretch into account (see next section)
an elastomechanical model

The interdependencies of the different data are depicted in figure 5. As an electrophysiological model the modified Noble-Kohl-Varghese-Noble model is used. Stretch dependent ion channels are included. The engaged force model is the Rice-Winslow-Hunter model (type 3). The elastomechanical behavior is modeled by affine transformations and the equations of Navier [49].

Extension of the Mono- and Bidomain Model

Figure 6. 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 this paper an extension of the mono- and bidomain model is introduced, which allows to take the deformation of tissue into account (see figure 6). 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 in the global coordinate system.

Model Assumptions

In the following simulations a different behavior due to stretch of the extra- and intracellular conductivity tensor is assumed:

The intracellular conductivity tensor σ_i is primarily influenced by the resistors R_gj of the gap junctions. These resistors R_gj are stretch independent.
The extracellular space behaves like an incompressible fluid. The extracellular conductivity tensor σ_e is not influenced by stretch.

These physically motivated assumptions simplify the behavior of the conductivity tensors and allow their efficient calculation.

Extraction of Stretch

The stretch of the myofibers in a voxel is determined by applying the deformation gradient tensor F (see Appendix Deformation Gradient Tensor) to the base of the material coordinate system M :

M' = FM

The base vectors of the material coordinate system M are the fiber orientation M_x , the sheet orientation M_y , and the sheet normal M_z :

M' = ( M_x | M_y | M_z )

The stretch s_x, s_y, and s_z in direction of the fiber, the sheet and the sheet normal respectively are the lengths of the transformed corresponding base vectors M'_x , M'_y , and M'_z :

s_x = √(M'_x²)
s_y = √(M'_y²)
s_z = √(M'_z²)

Construction and Scaling of Conductivity Tensor

The conductivity tensor influenced by stretch in a local coordinate system is determined by:

σ_s,local =
	(((s_x/s_ys_z) - 1)a_x + 1)σ_x	0	0
	0	(((s_y/s_xs_z) - 1)a_y + 1)σ_y	0
	0	0	(((s_z/s_xs_y) - 1)a_z + 1)σ_z

with the conductivity σ_x in direction of the fiber, σ_y in direction of the sheet, and σ_z in direction of the sheet normal. The weighting parameters a_x, a_y, and a_z allow to select a specific material behavior. For isovolumetric deformations, s_xs_ys_z = 1 the conductivity tensor reduces to:

σ_s,local =
	((s²_x - 1)a_x + 1)σ_x	0	0
	0	((s²_y - 1)a_y + 1)σ_y	0
	0	0	((s²_z - 1)a_z + 1)σ_z


$\begin{figure}\center\epsfig{file=cubedeform.eps,width=.6\textwidth}\end{figure}$

Figure 7. Cube and deformed cube. The undeformed cube is parameterized by a conductivity σ, width w, height h, and depth d. The deformed cube's parameters are a conductivity σ', width w', height h', and depth d'.

Examples for Weighting Parameters

To illustrate the influence of the weighting parameters, two cases are exemplary examined. The examples base on a cube with a conductivity σ (in x-direction), a height h, a width w = h, and a depth d = h. Hereby, the resistor R (in x-direction) of the cube is determined by

R =

1
hσ

The cube is isovolumetric, transversal isotropic deformed with a stretch s, which delivers a quader of width w' = hs, height h' = h/√s, and depth d' = d/√s. The conductivity of the quader σ' is scaled with the weighting factor a to

σ' = ((s² - 1)a + 1)

The resistor of the deformed cube R' is calculated by:

R' =

s²
hσ'

s²
h((s2 - 1)a + 1)σ

Case 1: a = 0. The influence of stretch is inhibited:

R'_{a = 0}

s²
hσ

resulting in a quadratic dependency of a corresponding resistor with respect to the stretch.
Case 2: a = 1. The corresponding resistor is independent on stretch:

R'_{a =1}

1
hσ

Coordinate System Transformation

A polar decomposition of the deformation gradient tensor F delivers the rotation matrix R_F, which is needed in conjunction with the rotation of the material coordinate system R_M to construct a summary rotation matrix R:

R = R_M R_F

The summary rotation matrix R is applied to transform the conductivity tensor σ_s,local from the local coordinate into the global coordinate system σ_s,global:

σ_s,global = Rσ_s,local R^T

Results

Deformation Induced Changes of Propagation Velocity

$\begin{figure}\center\epsfig{file=stretchLattice.eps,width=.5\textwidth}\end{figure}$

Figure 8. Lattice for determination of the deformation induced changes of the propagation velocity. A stimulus is applied at plane z =

. The velocity is determined outgoing from measurement of the transmembran voltage at the reading points P

and P

Numerical experiments are performed to determine the deformation induced changes of the propagation velocity. Hereby, a lattice (see Figure 8) of 13 x 13 x 25 cubic voxels was isovolumetrically deformed by a scaling matrix U:

	s_f	0	0
U =	0	s_t	0
	0	0	s_t

The scaling factor in fiber direction s_f ranges between 0.8 and 1.2. The transversal scaling factor s_t is chosen to achieve isovolumetry. Depending on the calculation of the longitudinal and transversal velocity the fibers have z- and y-direction, respectively. The initial size of a voxel is 0.1 mm x 0.1 mm x 0.1 mm. An excitation was initiated by applying a sufficiently large intracellular stimulus current at slice z = 1, which leads to a plane front of excitation wandering in z-direction.

The velocity of the excitation propagation is calculated by usage of two reading points. The excitation of a voxel is decided by defining a threshold to the transmembran voltage V_m. The threshold value is set to 0 mV. The time delay between the first exceed at the reading points is divided by their spatial distance:

v =	t_{V_m(P}_2)>0 – t_{V_m(P}_{1)>0 —————————— \|P1 - P2\|}

The reading points P1 and P2 are located at the positions (6, 6, 15) and P2 = (6, 6, 19), respectively.

Different environments are examined:

traditional bidomain model and stretch independent electrophysiological model (Stretch=0, SigmaStretch=0)
extended bidomain model and stretch independent electrophysiological model (Stretch=0, SigmaStretch=1)
traditional bidomain model and stretch depended electrophysiological model with constant sarcomere length (Stretch=1, SL=2, SigmaStretch=0)
extended bidomain model and stretch dependent electrophysiological model with constant sarcomere length (Stretch=1, SL=2, SigmaStretch=1)
traditional bidomain model and stretch dependent electrophysiological model (Stretch=1, SigmaStretch=0)
extended bidomain model and stretch dependent electrophysiological model (Stretch=1, SigmaStretch=1)

In the environments with the extended bidomain model (SigmaStretch=1) the weighting factors a_x , a_y , and a_z , for the intracellular conductivity are set to 1. The weighting factors for the extracellular conductivity are set to 0. The extra- and intracellular conductivity tensors are chosen according to [43].

The stretch independent electrophysiological model of the first and second environment (Stretch=0) is derived from the Noble-Varghese-Kohl-Noble model by setting the stretch-activated ion conductances to 0. In the environments with the stretch dependent electrophysiological model and non-constant sarcomere length (Stretch=1) the sarcomere length SL is scaled with s_f.

$\begin{figure}\center\epsfig{file=longms.eps,width=.9\textwidth}\end{figure}$

Figure 9. Progression of the velocity of longitudinal propagation by variation of stretch in fiber direction s_x. Different environments with isovolumetric deformation are examined.


$\begin{figure}\center\epsfig{file=longvs.eps,width=.9\textwidth}\end{figure}$

Figure 10. Progression of the voxel velocity of longitudinal propagation by variation of stretch in fiber direction s_x. Different environments with isovolumetric deformation are examined.


$\begin{figure}\center\epsfig{file=transms.eps,width=.9\textwidth}\end{figure}$

Figure 11. Progression of the velocity of transversal propagation by variation of stretch in fiber direction s_x. Different environments with isovolumetric deformation are examined.


$\begin{figure}\center\epsfig{file=transvs.eps,width=.9\textwidth}\end{figure}$

Figure 12. Progression of the voxel velocity of transversal propagation by variation of stretch in fiber direction s_x. Different environments with isovolumetric deformation are examined.

The results of the numerical experiments, the calculated velocities of excitation propagation influenced by stretch, are depicted in figure 9, 10, 11, and 12. Longitudinal and transversal velocities are determined.
The results of the simulations with the extended bidomain model and with the stretch dependent electrophysiological model are significant different from the results with the traditional bidomain model. The tendency of decreasing the propagation velocity by increasing the stretch in fiber direction observable with the traditional bidomain model is inverted with the extended model. The influence of the stretch dependent electrophysiological model is starting at a stretch larger than 1 and leads to an approximately exponential progression of the propagation velocity in the selected range of stretch. This effect is resulting from the raise of the resting potential by stretch with the consequence that a smaller and therefore faster deliverable amount of electrical charge is sufficient to initiate an excitation.

Deformation Induced Initiation of Excitation

The initiation of an excitation by stretch is simulated with the extended bidomain model and the modified Noble-Varghese-Kohl-Noble stretch dependent model. The lattice includes 16 x 16 x 16 voxels. The initial size of a voxel is 0.1 mm x 0.1 mm x 0.1 mm. The conductivity tensors are constructed to model an alignment of the fibers in x-direction.

The lattice is deformed by usage of Navier's equation [49]. The Lame coefficients are chosen to assure the approximation of an isovolumetric deformation. The eight corners of the lattice are fixed to their position. Forces are applied to the voxels at position (6,6,15), (6,7,15), (7,6,15), and (7,7,15) in negative x-direction.

Figure 13. Initiation of an electrical excitation in the myocardium by stretch. The surface with transmembran potential equal to 0 mV is visualized using the Marching Cube Algorithm [50].

The forces act on the lattice until an excitation front is initiated (see Figure 13 and movie). At first the excitation picks up only some cells near to the force sources. These cell receive the largest stretch. This small group of cells cannot deliver enough charge to excite their neighbors. After further acting of the forces a sufficiently large area is excited and an excitation front is initiated. At the beginning the excitation front moves primarily in direction of the fibers with only a small radial component. After reaching the border at plane

the front wanders in radial direction. The excitation includes the complete lattice.

Electromechanics in the Ventricular Heart Wall

The electromechanics in a region of the ventricular heart wall is simulated with the classic bidomain, the Noble-Varghese-Kohl-Noble stretch independent and the Rice-Winslow-Hunter force development model. The lattice includes 64 x 64 x 64 voxels. The size of a voxel is 0.2 mm x 0.2 mm x 0.2 mm.

The orientation of the fibers varies from the subepicardial to the subendocardial myocardium. An angle of -70 $^\circ$ is assigned for the orientation at the ventricular subepicardial myocardium, an angle of 70 $^\circ$ at the subendocardial myocardium [1][2]. The orientation in the space lying in between is interpolated outgoing from these boundary conditions by iterative averaging [51].

A modification of the K⁺ depolarizing current I_to of the electrophysiological model is integrated. The current is contingent on the distance to the epicardial and endocardial border. Three different regions are distinguished: subepicardial, midmyocardial, and subendocardial [52].

The excitation of the heart wall is initiated at the subendocardial myocardium by applying pointwise a sufficiently large intracellular current. The applying of current starts at apical points modeling myocytes with connections to Purkinje fibers and wanders basal.

Figure 14: Excitation propagation of the heart wall. The progression of the transmembran potential V_m in a three dimensional model of the ventricular free wall is simulated. The excitation is initiated at endocardial points modeling myocyte-Purkinje fiber connections.

Figure 15: Electrophysiology of the heart wall. The progression of the concentration of intracellular calcium

C_a

V_m

Figure 16: Electromechanics in the heart wall. The progression of the force development in a model of the heart wall is simulated and visualized in combination with the transmembran potential V_m and intracellular calcium concentration

C_a

The excitation propagation in a ventricular free wall is described in Figure 14 and movie by visualizing the transmembran potential V_m. The distribution of the concentration of intracellular calcium [C_a]_i is depicted in Figure 15 and movie in conjunction with the transmembran potential V_m. The developed force is shown in Figure 16 and movie in conjunction with the transmembran potential V_m and the concentration of intracellular calcium [C_a]_i .

Conclusions

The presented model describes some aspects of the electromechanical behavior of a myocardial region. The model combines an electrophysiological, a force development and an excitation propagation model. All of these models incorporate the effects of deformation of the myocardium.

An extension of the traditional bidomain model for excitation propagation is proposed. The extension describes the stretch dependency of the conductivity tensor of the intra- and extracellular space and is constructed outgoing from physically motivated assumptions, which simplify the behavior of the conductivity tensor. The extension makes usage of the deformation gradient tensor, which is a foundation in the theory of continuums mechanics.

The performed simulations illustrate some effects of myocardial electromechanical behavior. These effects are of great significance for the development of realistic models of the hole heart. The extension of models of cardiac excitation propagation allows mechanical interaction. This functionality is of importance eg for studying of cardiac arrhythmias and for the computer aided planning of surgical interventions.

Acknowledgments

The authors want to thank Miss P. Noble, University Laboratory of Physiology, Oxford, UK, for her aid to parameterize the electrophysiologic model. Dr. J. Rice, Johns Hopkins University School of Medicine, Baltimore, US, helped us to implement the force model. R. Mayer, Rechenzentrum, Univerität Karlsruhe (TH), supported our work by providing the necessary visualization and computing resources.

This work is supported in parts by SFB 414 "Information Technology in Medicine: Computer- and Sensor-based Surgery", a co-operation between the University of Karlsruhe, the University of Heidelberg and the German Cancer Research Center (DKFZ).

Appendix

Deformation Gradient Tensor

The finite continuous medium Ω is deformed outgoing from the reference configuration R₀ at time t = 0 [53]. The following configurations R_t are defined at time t . The coordinates of a point P in Ω are described at the reference configuration R₀ by the vector X . At time t the point's coordinates are:

x = x(X,t)

The deformation gradient tensor F is defined by differentiating x(X,t) with respect to the initial coordinates X :

F =

δx_i
δX_j

δx₀ δX₀	δx₀ δX₁	δx₀ δX₂
δx₁ δX₀	δx₁ δX₁	δx₁ δX₂
δx₂ δX₀	δx₂ δX₁	δx₂ δX₂

The tensor F is of second order. It converts an elementary segment dX of the reference configuration R₀ into a segment dx in the configuration R_t . The deformation gradient tensor can be divided in a rotation tensor R and a right stretch tensor U :

F = RU

The right stretch tensor U can be represented by a diagonal matrix.

Electrical Conductivity Tensor

In a stationary electrical current field the current density J is coupled with the gradient of the potential ϕ by the conductivity tensor σ:

J = – σ∇ϕ

Commonly in bioelectrical tasks an orthotrophy of the material properties is assumed. The material properties are described by a symmetric conductivity tensor $\sigma$ of second order. It has the following form in a local coordinate system:

	σ_x	0	0
σ_local =	0	σ_y	0
	0	0	σ_z

with a scalar conductivity σ_x in direction of x-axis, σ_y in direction of the y-axis, and σ_z in direction of the z-axis. In the transversal isotropic case the diagonal element σ_y is equal to σ_z . In the isotropic case all diagonal elements are equal: σ_x = σ_y = σ_z . Hereby, the tensor σ can be substituted by a scalar conductivity.

Implementation

$\begin{figure}\center\epsfig{file=overview.eps,width=1\textwidth}\end{figure}$

Figure 17. Implementation.

The model is implemented in the C++ programming language using interprocess communication mechanisms. The following data are kept in shared memory:

states of cellular electrophysiologic models
states of cellular force models
extracellular potentials
intra- and extracellular conductivity tensors
deformations
elastomechanic material parameters

The program can take advantage of multiprocessor computers with shared memory architecture using OpenMP. The modules communicate via message passing. Semaphores are used for synchronization of processes.

References

[1] D. D. Streeter, Jr. and D. L. Bassett, ``An engineering analysis of myocardial fiber orientation in pig's left ventricle in systole,'' Anatomical Record, vol. 155, pp. 503-512, 1966.

[2] D. D. Streeter, ``Gross morphology and fiber geometry of the heart,'' in Handbook of Physiology: The Cardiovascular System (B. Bethesda, ed.), vol. I, pp. 61-112, American Physiology Society, 1979.

[3] I. J. LeGrice, B. H. Smaill, L. Z. Chai, S. G. Edgar, J. B. Gavin, and P. J. Hunter, ``Laminar structure of the heart: Ventricular myocyte arrangement and connective tissue architecture in the dog,'' Am J Physiol, vol. 269, pp. H571-H582, 1995.

[4] D. M. Bers, Excitation-Contraction Coupling and Cardiac Contractile Force.
Dordrecht, Netherlands: Kluwer Academic Publishers, 1991.

[5] J. E. Saffitz and K. A. Yamada, ``Gap junction distribution in the heart,'' in Cardiac Electrophysiology. From Cell to Bedside (D. P. Zipes and J. Jalife, eds.), ch. 21, pp. 271-277, Philadelphia: W. B. Saunders Company, 1999.

[6] R. H. Hoyt, M. L. Cohen, and J. E. Saffitz, ``Distribution and three-dimensional structure of intercellular junctions in canine myocardium,'' Circ Res., vol. 64, pp. 563-574, 1989.

[7] J. B. Caulfield and T. Borg, ``The collagen network of the heart,'' Laboratory Investigation, vol. 40, no. 3, pp. 364-372, 1979.

[8] C. Abrahams, J. S. Janicki, and K. T. Weber, ``Myocardial hypertrophy in macaca fascicularis: Structural remodeling of the collagen matrix,'' Laboratory Investigation, vol. 56, pp. 676-683, 1987.

[9] K. T. Weber, Y. Sun, S. C. Tyagi, and J. P. M. Cleutjens, ``Collagen network of the myocardium: Function, structural remodeling and regulatory mechanisms,'' J Mol Cell Cardiol, vol. 26, pp. 279-292, 1994.

[10] A. L. Hodgkin and A. F. Huxley, ``A quantitative description of membrane current and its application to conduction and excitation in nerve,'' J. Physiol, vol. 177, pp. 500-544, 1952.

[11] G. W. Beeler and H. Reuter, ``Reconstruction of the action potential of ventricular myocardial fibres,'' J. Physiol., vol. 268, pp. 177-210, 1977.

[12] D. DiFrancesco and D. Noble, ``A model of cardiac electrical activity incorporating ionic pumps and concentration changes,'' Phil. Trans. R. Soc. Lond., vol. 307, pp. 353-398, 1985.

[13] Y. E. Earm and D. Noble, ``A model of single atrial cell: Relation between calcium current and calcium release,'' Proc. R. Soc. Lond., vol. 240, pp. 83-96, 1990.

[14] C.-H. Luo and Y. Rudy, ``A model of the ventricular cardiac action potential,'' Circ Res., vol. 68, no. 6, pp. 1501-1526, 1991.

[15] C.-H. Luo and Y. Rudy, ``A dynamic model of the ventricular cardiac action potential: I. simulations of ionic currents and concentration changes,'' Circ Res., vol. 74, no. 6, pp. 1071-1096, 1994.

[16] C.-H. Luo and Y. Rudy, ``A dynamic model of the ventricular cardiac action potential: II. afterdepolarizations, triggered activity, and potentiation,'' Circ Res., vol. 74, no. 6, pp. 1097-1113, 1994.

[17] S. S. Demir, J. W. Clark, C. R. Murphey, and W. R. Giles, ``A mathematical model of a rabbit sinoatrial node cell,'' Am. J. Physiol., vol. 35, pp. 832-852, 1994.

[18] S. S. Demir, B. O'Rourke, G. F. Tomaselli, E. Marban, and R. L. Winslow, ``Action potential variation in canine ventricle: A modeling study,'' in Proc. Computers in Cardiology, pp. 221-224, 1996.

[19] A. Nygren, C. Fiset, L. Firek, J. W. Clark, D. S. Lindblad, R. B. Clark, and W. R. Giles, ``Mathematical model of an adult human atrial cell,'' Circ Res., vol. 82, pp. 63-81, 1998.

[20] D. Noble, A. Varghese, P. Kohl, and P. Noble, ``Improved guinea-pig ventricular cell model incorporating a diadic space, ${I}_{Kr}$ and ${I}_{Ks}$ , and length- and tension-dependend processes,'' Canadian Journal of Cardiology, vol. 14, pp. 123-134, Jan. 1998.

[21] M. S. Jafri, J. J. Rice, and R. L. Winslow, ``Cardiac ${C}a^{2+}$ dynamics: The roles of ryanodine receptor adapation and sarcoplasmic reticulum load,'' Biophysical J, vol. 74, pp. 1149-1168, Mar. 1998.

[22] R. L. Winslow, J. J. Rice, S. Jafri, E. Marbán, and B. O'Rourke, ``Mechanisms of altered excitation-contraction coupling in canine tachycardia-induced heart failure, II model studies,'' Circ. Res., vol. 84, pp. 571-586, 1999.

[23] B. O'Rourke, D. A. Kass, G. F. Tomaselli, S. Kaab, R. Tunin, and E. Marban, ``Mechanisms of altered excitation-contraction coupling in canine tachycardia-induced heart failure, I experimental studies,'' Circ Res, vol. 84(5), pp. 562-570, 1999.

[24] E. White, J.-Y. L. Guennec, J. M. Nigretto, F. Gannier, J. A. Argibay, and D. Garnier, ``The effects of increasing cell length on auxotonic contractions; membrane potential and intracellular calcium transients in single guninea-pig ventricular myocytes,'' Experimental Physiology, pp. 65-78, 1993.

[25] A. Landesberg and S. Sideman, ``Coupling calcium binding to troponin C and cross-bridge cycling in skinned cardiac cells,'' American Journal of Physiology, vol. 266, pp. H1260-H1271, 1994.
Heart Circ. Physiol. 35.

[26] A. Landesberg and S. Sideman, ``Mechanical regulation of cardiac muscle by coupling calcium kinetics with cross-bridge cycling: A dynamic model,'' American Journal of Physiology, vol. 267, pp. H779-H795, 1994.
Heart Circ. Physiol. 36.

[27] J. J. Rice, R. L. Winslow, and W. C. Hunter, ``Comparison of putative cooperative mechanisms in cardiac muscle: length dependence and dynamic responses,'' Am. J. Physiol. Circ. Heart., vol. 276, pp. H1734-H1754, 1999.

[28] J. J. Rice, M. S. Jafri, and R. L. Winslow, ``Modeling short-term interval-force relations in cardiac muscle,'' Am. J. Physiol. Circ. Heart., vol. 278, pp. H913-H931, 2000.

[29] W. J. Eifler and R. Plonsey, ``A cellular model for the simulation of activation in the ventricular myocardium,'' J. Electrocardiology, vol. 8, no. 2, pp. 117-128, 1975.

[30] R. Killmann, P. Wach, and F. Dienstl, ``Three-dimensional computer model of the entire human heart for simulation of reentry and tachycardia: Gap phenomenon and Wolff-Parksinson-White syndrome,'' Basic Research in Cardiology, vol. 86, no. 5, 1991.

[31] B. E. H. Saxberg and R. J. Cohen, ``Cellular automata models of cardiac conduction,'' in Theory of Heart (L. Glass, P. Hunter, and A. McCulloch, eds.), pp. 437-476, Berlin, Heidelberg, New York: Springer, 1991.

[32] D. Wei, O. Okazaki, K. Harumi, E. Harasawa, and H. Hosaka, ``Comparative simulation of excitation and body surface electrocardiogram with isotropic and anisotropic computer heart models,'' IEEE Transactions on Biomedical Engineering, vol. 42, pp. 343-357, Apr. 1995.

[33] P. Siregar, J. P. Sinteff, M. Chahine, and P. L. Beux, ``A cellular automata model of the heart and its coupling with a qualitative model,'' Computers and Biomedical Research, vol. 29, pp. 222-246, 1996.

[34] C. D. Werner, F. B. Sachse, and O. Dössel, ``Applications of the visible man dataset in electrocardiology: Simulation of the electrical excitation propagation,'' in Proc. Second Users Conference of the National Library of Medicine's Visible Human Project, 1998.

[35] P. Siregar, J. P. Sinteff, N. Julen, and P. L. Beux, ``An interactive 3D anisotropic cellular automata model of the heart,'' Computers and Biomedical Research, vol. 31, pp. 323-347, 1998.

[36] R. A. FitzHugh, ``Impulses and physiological states in theoretical models of nerve membran,'' Biophys J, vol. 1, pp. 445-466, 1961.

[37] J. M. Rogers and A. D. McCulloch, ``A collocation-Galerkin finite element model of cardiac action potential propagation,'' IEEE Transactions on Biomedical Engineering, vol. 41, pp. 743-757, Aug. 1994.

[38] A. V. Panfilov, ``Three-dimensional wave propagation in mathematical models of ventricular fibrillation,'' in Cardiac Electrophysiology. From Cell to Bedside (D. P. Zipes and J. Jalife, eds.), ch. 31, pp. 271-277, Philadelphia: W. B. Saunders Company, 1999.

[39] Y. Rudy and W. Quan, ``Mathematical model of reentry of cardiac excitation,'' in Computers in Cardiology, pp. 135-136, 1989.

[40] N. Virag, O. Blanc, J. M. Vesin, J. Koerfer, and L. Kappenberger, ``Study of the mechanisms of arrhythmias in an anatomical computer model of human atria,'' in Proc. Computers in Cardiology, pp. 113-116, 1999.

[41] C. S. Henriquez and R. Plonsey, ``A bidomain model for simulating propagation in multicellular cardiac tissue,'' in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 4, p. 1266, 1989.

[42] N. G. Sepulveda and J. P. Wikswo, ``Bipolar stimulation of cardiac tissue using an anisotropic bidomain model,'' Cardiovascular Electrophysiology, vol. 5, pp. 258-267, May 1994.

[43] C. S. Henriquez, A. L. Muzikant, and C. K. Smoak, ``Anisotropy, fiber curvature and bath loading effects on activation in thin and thick cardiac tissue preparations: Simulations in a three-dimensional bidomain model,'' J Cardiovascular Electrophysiology, vol. 7, pp. 424-444, May 1996.

[44] M. S. Spach and J. F. Heidlage, ``The stochastic nature of cardiac propagation at a microscopic level,'' Circ. Res., vol. 76, no. 3, pp. 118-130, 1995.

[45] M. S. Spach, J. F. Heidlage, and P. C. Dolber, ``The dual nature of anisotropic discontinous conduction in the heart,'' in Cardiac Electrophysiology. From Cell to Bedside (D. P. Zipes and J. Jalife, eds.), ch. 25, pp. 213-222, Philadelphia: W. B. Saunders Company, 1999.

[46] N. Hooke, C. S. Henriquez, P. Lanzkron, and D. Rose, ``Linear algebraic transformations of the bidomain equations: Implications to numerical methods,'' Crit Rev Biomed Eng, vol. 120, pp. 127-145, 1992.

[47] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C.
Cambridge, New York, Melbourne: Cambridge University Press, 2 ed., 1992.

[48] K.-J. Bathe, Finite Element Procedures in Engineering Analysis.
Englewood Cliffs/NJ: Prentice Hall, 1982.

[49] Y. C. Fung, Biomechanics: Motion, Flow, Stress and Growth.
New York, Berlin, Heidelberg: Springer, 1990.

[50] W. E. Lorensen and H. E. Cline, ``Marching cubes: A high resolution 3D surface construction algorithm,'' Computer Graphics, vol. 21, no. 4, pp. 163-169, 1987.

[51] F. B. Sachse, M. Wolf, C. Werner, and K. Meyer-Waarden, ``Extension of anatomical models of the human body: Three dimensional interpolation of muscle fiber orientation based on restrictions,'' Journal of Computing and Information Technology, vol. 6, no. 1, pp. 95-101, 1998.

[52] D.-W. Liu, G. A. Gintant, and C. Antzelevitch, ``Ionic bases for electrophysiological distinctions among epicardial, midmyocardial, and endocardial myocytes from the free wall of the canine left ventricle,'' Circ Res., vol. 72, no. 3, pp. 671-687, 1993.

[53] W. Maurel, Y. Wu, N. Magnenat Thalmann, and D. Thalmann, Biomechanical Models for Soft Tissue Simulation.
Berlin: Springer, 1998.