CELLULAR AUTOMATON
MODELS OF THE HEART
by
Anne Spalding
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
1996
This thesis for the Master of Science
degree by
Anne Spalding
has been approved
by
William Briggs
Stephen Billups
Weldon Lodwick
(t March
Date
Spalding, Anne (M.S., Applied Mathematics)
Cellular Automaton Models of the Heart
Thesis directed by Professor William Briggs
ABSTRACT
Understanding the mechanisms of the heart is essential to the diagnosis and cure of
many cardiac disorders. Cellular automaton models represent the cells of the heart
with an integer grid and model the ionic interactions of these cells. Cellular automaton
models provide a tool for understanding the correlation between an ECG
(electrocardiogram) and the actual movement of charge through the heart. These
models are used to gain a better understanding of the mechanisms of activation of the
cells of the heart. They can also be used as a diagnostic tool in modeling and
predicting cardiac disorders. This thesis develops a cellular automaton model of the
ventricles of the heart. Discrete rules are developed to model cellular interaction and
recreate the activation sequence of the heart. The model is then used to produce
ECGs corresponding to several known cardiac disorders in order to show the
effectiveness of the model in representing actual cardiac phenomena.
This abstract accurately represents the content of the candidates thesis. I recommend
its publication.
Signed
William Briggs
ACKNOWLEDGMENTS
I would like to take this opportunity to thank Professor Briggs and Professor Hermes
for their guidance in completing this research. I would also like to thank my friends
and my family for their constant support and encouragement in all of my pursuits.
CONTENTS
1. Introduction.................................................... 1
2. Physiology of the Heart........................................... 3
2.1 Anatomy...................................................... 3
2.2 Electrophysiology .......................................... 4
2.3 Modeling Ionic Interchanges ................................ 8
3. The Cellular Automaton Model.................................... 13
3.1 Model Setup ............................................... 13
3.2 Depolarization Algorithm................................... 16
3.3 Repolarization Algorithm................................... 18
3.4 Output..................................................... 21
4. Experiments on the Twave........................................ 22
4.1 Model A .................................................... 23
4.2 Model B ................................................... 24
4.3 Model C ................................................... 25
4.4 Model D ................................................... 26
5. Experiments...................................................... 29
5.1 SelfSustained Ventricular Tachycardia (SSVT)............... 29
5.2 Ventricular Tachycardia.................................... 32
5.3 WolfFParkinsonWhite Syndrome (WPW)....................... 34
5.4 Right and Left Bundle Branch Block ........................ 35
6 Conclusions...................................................... 37
7 References....................................................... 39
v
FIGURES
Figure
2.1 Anatomy and conduction system of the heart ......................... 4
2.2 Ionic currents and channels in an action potential [3].............. 5
2.3 Standard 12Lead ECG placement[7]................................... 6
2.4 Formation of the ECG during the activation sequence of the heart .... 7
2.5 The normal IILead ECG for a SA pulse.............................. 8
2.6 Representation of the septum and the left ventricle of the heart used in
an ionic model ..................................................... 10
3.1 Cross section of the CA model and the model setup for one cell..... 15
3.2 Elliptical wavefronts on a cross section of a modified MillerGeselowitz
model [9] ......................................................... 18
3.3 Schematic representation of repolarization of a target cell in the CA
model............................................................... 20
4.1 ECG for a normal SA pulse without Pwaves.......................... 23
4.2 ECG for a normal SA pulse from Model A: inverted Twave ........... 24
vi
4.3 ECG for a normal SA pulse from Model C: inverted Twave.......... 26
4.4 Repolarization intervals on an xz cross section of Model D...... 27
4.5 ECG for a normal SA pulse from Model D: Twave is correct ........28
5.1 SSVT for Model A with one ectopic and one SA pulse............... 30
5.2 Nonrhythmic SSVT for Model D.................................... 31
5.3 ECG of SSVT for Model D with a capture beat occurring at 1860 ms 32
5.4 Ventricular tachycardia initiated in the upper left ventricle.... 33
5.5 Ventricular tachycardia initiated in the lower right ventricle .. 33
5.6 Lateral slice of the ventricles showing areas for accessory pathways
inWPW............................................................. 34
5.7 Lead I, II, and III ECGs for WPW in a modified MillerGeselowitz
model [9]......................................................... 35
5.8 ECG for simulated right bundle branch block and a normal SA pulse
in Model D
1. Introduction
Cardiac disorders such as infarction, tachycardia, and fibrillation are the cause
of thousands of critical care situations eveiy year. One of the main methods of
diagnosing and analyzing these disorders is through an electrocardiograph (ECG) scan
which records the electronic changes in the heart tissue. The effectiveness of the ECG
in diagnosis is hindered by the fact that there is still an incomplete understanding of
the relationship between the output of the ECG and the changes that actually occur
on a cellular level.
This thesis presents a mathematical model called a cellular automaton (CA)
that helps to explore this relationship. This cellular automaton model, represents the
ventricles of the heart by 2605 points on a threedimensional integer grid. These
points are given attributes and interaction rules so that they resemble the actual cells
of the heart. The interaction of the grid points in the cellular automaton model is
represented by a simulated lead II ECG scan which can be compared to ECG scans
from the actual heart. The main advantage of this model is the ease of simulation
and the ability to relate ECG output with the point of origin and path of electrical
pulses. The cellular automaton is effective in modeling several disorders of the heart
and reproducing the ECG scans for irregular heart activations. Analysis of the ECG
output of the CA model in comparison with that of patients may prove to be a valuable
tool in diagnosis and treatment of these disorders.
The second chapter of this paper discusses the underlying anatomy necessary
for understanding the functions of the heart in order to develop an accurate cellular
automaton model. This includes a brief description of the action potential and the
propagation of action potential through the conduction system of the heart. In the
next chapter one particular CA model is developed in detail. This model simulates
both the depolarization and repolarization sequences in the ventricles. In the fourth
1
chapter, some of the limitations of various CA models such as the inverted Twave are
discussed. Several possible solutions to the discrepancies between the models and the
actual heart are addressed and evaluated. In the last chapter experiments are presented
that show the use of the CA models in simulating several cardiac conditions. These
experiments include ventricular tachycardia, capture beats, WolfFParitinsonWhite
syndrome and right bundle branch block.
The CA model developed in this thesis was originally created by Professor Henry
Hermes of the University of Colorado at Boulder [13]. The research contribution
of this thesis is to modify the code and algorithms of the original model and use
this modified model in simulations of cardiac disorders. This thesis also uses at new
methods in the repolarization algorithm to analyze the inverted Twave which is a
limitation in the original model.
2
2. Physiology of the Heart
2.1 Anatomy
The human heart has four chambers. The right and left atria are the receiving
chambers; the right atrium receives blood from the body and the left atrium receives
blood from the lungs. From the atria the blood flows into the right and left ventricles
which are separated by a thick muscle wall called the septum. The right ventricle
pumps blood out to the lungs while the left ventricle pumps blood to the entire
circulatory system of the body. Thus the left ventricle walls, which contract to
pump the blood, are three times as thick as in the right ventricle. The wall of each
chamber consists of three layers. The inner layer is the endocardium, the middle
layer is the myocardium and is composed of cardiac muscle, and the outer layer is the
epicardium [3], The electrical conduction system of the heart begins at the sinoatrial
(SA) node which is located above the right atrium; it then follows intemodal pathways
through the muscles of the atria to the atrioventricular (AV) node which is located
on the right side of the interatrial septum. In the ventricles, the conduction system
consists of the bundle of His, which begins at the AV node, and travels down the
septum until it branches into the right and left bundle branches [7], These branches
pass through the right and left ventricles until they branch out into the Purkinje fibers
which spread along the muscle walls at the bottom of the ventricles (Figure 2.1).
3
Figure 2.1 Anatomy and Conduction system of the Heart [3].
2.2 Electrophysiology
The pulsing of the heart is controlled by an electrochemical wave which passes
from cell to cell within the heart. In a normal beating heart, the electrical stimulus
originates in the SA node at the top of the atria and propagates through conductive
pathways to the Purkinje fibers. The electrochemical wave then spreads through the
myocardial muscle of the ventricles. This wave causes the sequential and periodic
contraction of the muscles in the heart which forces blood from the atria to the
ventricles and out into the circulatory system of the body. This electrochemical
wave is formed by the propagation of an action potential.
An action potential is the changing of polarity or charge of a cell due to ionic
exchanges between the extracellular fluid and the cell. The action potential of a cell
is defined by five phases. In phase 0 sodium enters the cell through fast sodium
channels which causes the polarity of the cell to change from its resting potential of
90 mV to a potential of +10 mV. This phase is referred to as the depolarization
of the cell. In phase 1 potassium leaves the cell and leads to a slight decrease in
the potential to about 0 mV. This stage is referred to as early repolarization. In
phase 2, which is the plateau, there is a balance between potassium leaving the
cell and calcium entering it. In phase 3, rapid repolarization, the cell continues
to repolarize as potassium leaves the cell and the potential decreases back towards
4
the resting potential. Finally in phase 4, return to resting potential, the amount
of potassium leaving the cell is only slightly higher than the amount entering it and
the potential levels off at 90 mV. During this phase the ionic balance of the cell is
actively restored by the sodium potassium pumps in the cell [3] (Figure 2.2). From
the end of phase 1 and through most of phase 3 the potential of the cell is too large
to allow the cell to depolarize again; this is called the absolute refractory period.
During the end of phase 3 and the beginning of phase 4, the cell can depolarize but
the height of the action potential will be less than normal; this is called the relative
refractory period.
K* ctartncis (}*, tb)
Figure 2.2 Ionic curents and channels in an action potential [3]. The chemical and electro
static movement of potassium (K+), sodium(Na+), and calcium(Ca++) are shown for each
phase of the action potential.
Any cell can begin and maintain a rhythmic firing of action potentials which is
referred to as automaticity. Certain cells of the heart, such as the SA node, have
an inherent automaticity while other cells develop automaticity causing abnormal
conditions in the heart. In a normal heart, the SA node produces an action potential
at regular intervals. The depolarization of the SA node changes the concentration
of sodium and potassium in the extracellular fluid around the SA node causing the
sequential depolarization of neighboring cells. Action potentials propagate quickly
along the conduction pathways and more slowly through the cardiac muscle cells,
5
so that the electrochemical wave will depolarize the entire heart in about 100 ms.
In addition, propagation through the cardiac muscle is anisotropic so that the action
potentials travel faster along the muscle fiber than across muscle fibers. If a cell other
than the SA node begins firing or producing action potentials, then this cell is called
an ectopic pulse or focus. A depolarization wave will spread through the heart from
an ectopic focus only by cell to cell propagation instead of traveling along the normal
conduction system as an SA pulse does. Depending on the position of the ectopic
pulse, the time to depolarize the entire heart can be longer or shorter than the 80100
ms required when the depolarization wave originates at the SA node.
One of the main methods of observing the electrophysiology of the heart is through
an electrocardiogram (ECG) which is an electronic record of the activity of the heart.
This activity is recorded by placing up to twelve electrodes on the surface of the body
and recording the difference in electrical potential between various combinations of
these electrodes.
Figure 2.3 The 12Lead ECG placement and the standard placement of the three main leads.
The bottom enlargement shows the placement of the six precordial leads [7J.
The specific placement of electrodes on the body (Figure 2.3) gives a standard
twelve lead ECG. The three standard leads (1,11,111) view the heart in the frontal
plane and are generally used for monitoring purposes. The three augmented leads
6
(aVR,aVL,aVF) use different combinations of the three standard electrodes in order
to record small deflections in the potentials. Finally, the six precordial leads (vl...v6)
view the heart in the horizontal plane. The pulsing of the SA node is not detectable in
the ECG, since the atria depolarize almost instantaneously after the firing of the SA
node. The average direction of the depolarization wave in the atria is down towards
the patients left side and is recorded in the Pwave of the ECG [7], The repolarization
of the atria is usually not visible in the ECG because it takes place at the same time
as the depolarization of the ventricles. The ventricular depolarization is represented
by the QRS complex in the ECG. In different leads this complex can contain any
combination of the Qwave, Rwave, and Swave. In all leads'^, the Qwave is
the negative wave preceding the positive Rwave, and the Swave is the negative
wave following the Rwave. Finally the Twave represents the repolarization of the
ventricles (Figure 2.4). The PR interval represents the time for atrial depolarization
and the pause at the AV node before ventricular depolarization. The QRS interval
represents the time for ventricular depolarization. The QT interval represents the
time between depolarization and the beginning of repolarization in the ventricles.
This interval varies inversely with heart rate and its length must be compared with a
normal heart to detect abnormalities.
ll ECG
v ^4
QRS%JF (mV) :i_4
uj V
Figure 2.4 Formation of the ECG during the activation sequence of the heart. The shaded
areas represent depolarized cells. The ULead ECG for each step in the depolarization process
is shown below the figure that shows which cells are depolarized [8J.
1
Each of these waves and intervals has a different normal state for each lead and
the comparison of all twelve leads is used to gain a complete picture of the electrical
activity and abnormalities in the heart. The II lead for a normal SA pulse has positive
R and Twaves with the Twave having about 1/3 the height and 3 times the length
of the Rwave (Figure 2.5).
Figure 2.5 The normal IILead ECG for a SA pulse. The ECG has a positive Rwave and a
positive Twave.
2.3 Modeling Ionic Interchanges
An ECG is one method used to understand the ionic interchanges between cells
based on the changes in polarity of portions of the heart, but this process can also be
understood through models that replicate cellular interactions. The first work towards
modeling the electrochemical interchange in cells was done by Hodgkin and Huxley
in the 1960s [14], By performing voltage clamp experiments on a squid giant axon,
they were able to develop a set of differential equations that model the action potential
within the cell and activation conduction between cells. These equations for a neuron
were then adapted in various ways to fit the action potentials of Purkinje fibers and
cardiac cells. In a study by McCallister, Nobel and Tsien, the adaptions were formed
by using similar voltage clamp techniques on the heart. The Purkinje fibers were
found to have nine different ionic currents [17][18] in contrast to the three currents
..^described in [9j for the neuron. These nine currents are shown in Table 2.1; the
specific time and voltage dependence of each current is described in detail in [17],
Table 2.1: Nine variables describing ionic currents in Purkinje Fibers
Type of Channel Variable Description
POTASIUM outward currents changing exponentially with step changes in membrane potential
ixi
SODIUM iNa fast outward current
^si slow outward current
CHLORIDE lqr variable activation and inactivation
BACKGROUND CURRENT iKl potasium outflow
iNa,b sodium influx during rest (Phase 4)
ici,b chloride flow during plateau
The main equation which relates the interaction of these ionic currents to the overall
membrane current is
dE i
dt Cm
(2.1)
where E is the membrane potential, h is the sum of the nine Purkinje currents and
Cm is the membrane capacitance. In a similar study conducted by Beeler and Reuter
[1], the action potentials of myocardial cells were found by defining equations for
five ionic channels. The ionic channels for myocardial cells are shown in Table 2.2,
the specific equations for voltage and time dependency of the ionic exchanges can be
found in [1],
9
Table 2.2: Five variables describing the ionic currents in myocardial cells.
Type of Channel Variable Description
SODIUM *JVo Fast sodium channel
Slow sodium channel
POTASSIUM iKl Time/voltage independant current
xl Time/voltage independant current
BACKGROUND 1 extern Externally applied current during experiment
The main equation relating the interactions of these ionic currents is similar to that
for the Purkinje fibers and is given as
dE (ijVa b is T" ifCl h ^xl extern)
~dt~ On ' ^ '
In [16] the action potential and excitation sequence in the ventricles are modeled by
combining these two ionic current models. The model in [16] is a Jshaped two
dimensional cut of the septum and left ventricle that has 730 cells (Figure 2.6). The
Purkinje fibers are in the uppermost layer and the remaining nine layers are defined
as myocardial cells. The following set of differential equations is used to describe
the ionic state of both Purkinje and myocardial cells located at coordinates as
well as the propagation of an action potential starting at a stimulus point and moving
through the cells of the model [16],
Figure 2.6 Representation of the septum and the left ventricle of the heart used in an ionic
model. P=Purkinje fiber, V=ventricular muscle [16].
10
dEi
dt
dEjj
dt
A,
M
D
where
f(Eij,Xij) + D Aij (change in membrane potential in Purkinje cells)
G(Eij,Xij) (change in membrane potential in myocardial cells)
Eiij + Ei+ij + Eij1 + Eij+1
(R9 S CJ1,
Eij : Membrane potential of unit cell with coordinates (i,j)
Xij : State variables for Purkinje and myocardial cells from tables 2.1 and 2.2
/, G : Functions describing relationship of Xij from eq. 2.1 and 2.2 respectively
Rg : Electrical resistance of gap junction between cells
S : Surface area of membrane of unit cell (larger in Purkinje)
Cm : Membrane capacitance
D : Coupling coefficient between unit cells.
The coupling coefficients and the initial conditions for the variables are assigned
numerical values estimated to correspond with experimental data [1][16][17], The
coupling coefficient D depends on the membrane capacitance Cm which is 12 times
higher in Purkinje fibers than in myocardial cells. Thus, in order to simulate symmet
rical interaction between cells, the value of D from Purkinje fibers to myocardial cells
is set to 1/12 of the value of D from myocardial cells to Purkinje fibers. Simulation
of one cycle of the entire model is done using numerical integration by Euler and
Runga Kutta methods. The surface potential for the simulation is then calculated and
an ECG of the simulation is formed. The ECG scans produced by the simulation cor
responded closely to actual ECGs indicating that the ionic equations are reasonably
accurate in describing the action potential propagation.
The main limitation of this model is the computational complexity of the numer
ical integration of the equations. Because of this complexity it is only possible to
11
complete one cycle of simulation, which greatly restricts the possibilities for model
simulation of activation phenomena of the heart. One way to avoid this computa
tional complexity is to work with the essential features of the differential equations
instead of the equations themselves. While this type of modeling may lose some of
the accuracy and detail of the ionic model, it is much more efficient in modeling
wholeheart phenomena and in analyzing the activation sequence over several peri
ods. The remainder of this thesis will focus on ways to simulate the characteristics
of the .activation sequence without actually solving the differential equations which
produce these characteristics.
12
3. The Cellular Automaton Model
A cellular automaton (CA) model is a set of discrete points which interact
based on a given set of rules. The aim of the CA model of the heart is to reproduce
accurately the depolarization and repolarization waves of the working heart by using
a set of biologically based rules on an anatomically correct discrete representation of
the cells of the myocardium. This is done by letting the discrete points of an integer
grid represent the cells of the myocardium. The rules governing the interaction of
these discrete points are derived from the end results of the equations used in the ionic
models. Using the differential equations of the ionic model, the polarization state of
each cell can be determined at any time step throughout the activation sequence.
In a CA model, the underlying equations are ignored and only the timedependent
states (repolarized or depolarized) of each cell are considered. Likewise, the ionic
model can determine the cell which will depolarize in each time step in relation to the
cells that have depolarized previously. Thus, since the ionic model determines that a
depolarization wave will travel from one cell to a certain number of surrounding cells
in 10 ms, the corresponding rule for the CA is that this number of cells surrounding
a depolarized cell will depolarize 10 ms or one time step after the initial cell. Thus
the time step used in the model is 10 ms. The precise interactions between the cells
of the CA model are described in the following sections.
3.1 Model Setup
Models of the activation sequence of the heart were first used in the 1920s
[20], Since this time, the introduction of CA models has improved the original
models so that currently, CA models can produce very accurate simulations of the
activation waves in the heart. The majority of these models are ventricular models
only. Although some detail is lost by ignoring the atria, the nature of the activation
13
process and the insulation of the ventricles from the atria allow ventricular models to
retain many of the qualities of a full heart model while reducing the computational
demands. The main differences among ventricular models is the accuracy of the
geometrical description of the ventricles and the rules of interaction both within and
between the cells.
The main model we will be looking at was developed by Hermes [13] in order to
analyze the interaction between ectopic and SA pulse firing. This model deals with a
ventricular heart in which the ventricles are assumed to be electrically insulated from
the atria. The epicardium is formed by the ellipsoid
(x ll)2 (y11)2 (z 16)2
+
100 100
The endocardium is formed by the paraboloid
+
225
= 1.
(3.1)
Or 10)2 + (y ll)2 = 4.9(z 4)2
for
z < 15.
(3.2)
The apex of the paraboloid is offcenter of the apex of the ellipsoid to make the left
ventricular wall about three times thicker than the right. The septum is formed by
the intersection of the inside of this paraboloid and the cylinder
132 < (re 25)2 + (y ll)2 < 172. (3.3)
The entire model lies in the first octant of a Cartesian coordinate system where
1 < x < 21, 1 < y < 21, and 1 < z < 15. A cross section of the model in the
x z plane is shown in Figure 3.1. These curves overlay an integer grid so that 2605
points with integer coordinates lie in the myocardial fiber between the ellipsoid and
the paraboloid or in the septum; these points represent cells of the heart. Each cell is
considered to be a neighbor of another cell if it has a Euclidean distance of \fZ or
less from another cell; thus each cell has at most 26 neighbors. The geometry is such
that if a wave passes from a cell to its neighbors in one time step, the QRS complex
will have a duration of 8 time steps (80 ms) when the heart is beating at 60 beats per
minute [13].
Each cell has a record in a threedimensional array, whose indices correspond
to the integer coordinates of the cell. In this record, cell type (epicardium, Purkinje,
14
conduction system, etc.), the layer position of the cell, the polarization state of the
cell and various other characteristics are recorded. The array and record for a sample
cell is shown in Figure 3.1.
Endocardium
Cell Position Ceil Attributes
Array Record
(I6llf9  p=0: depolarized
t=1 :time after depolarization
puft=l Purklnjc Fitter
u=Q not repolarizing
reuser set + 5: repolarlzatlan
interval
ium
Figure 3.1 Cross section of the CA model and the model setup for one cell. Each cell is
represented by a position in an array of records. The indices of the array are the coordinates
of the cell and the record keeps track of the cell attributes.
The general shape of this heart model is similar to most other CA heart models.
The main differences between this and other ventricular models are the definition of a
neighborhood of a cell and the number of cells in the model. In the MillerGeselowitz
model [9] the general outline of the model is the same, but there are 50,000 cells.
This increase in the number of cells leads to more realistic output and gives a better
resolution of small changes in the depolarization and repolarization of the heart. In
the Japanese model [11], the atria are included and the number of cells is 75,000.
The cells in the Japanese model are also on a tilted plane so that the number of
planar neighbors of each cell is increased from eight to twelve. Each cell in this
model has 16 characteristics associated with it, 10 of which define the shape and
duration of the action potential at any time. This provides a more accurate picture of
the activation sequence because it allows partial depolarization to occur during the
relative refractory period as was described in section 2.2. Also, the addition of the
15
atria in this model provides a better simulation of the depolarization wave from the
SA node and allows for investigation of atrial arrythmias and other disorders.
3.2 Depolarization Algorithm
In a normal heart, the depolarization of the ventricles starts with the firing of the
AV node. The pulse then travels down the bundle of His, separates into the right and
left bundle branches and travels to the Purkinje fibers. In the ventricular model, there
is no AV node so normal depolarization is started by firing six cells in the septum,
placed experimentally to obtain the best agreement with the wave travel in the actual
heart. In order to simulate the conduction system, the Purkinje fibers are set to fire
four time steps after the firing of the initial cells in the septum. This makes the firing
of the Purkinje fibers coincide with their firing if the depolarization wave traveled
through the conduction system, as it does in an actual heart.
The firing or depolarization of a cell with coordinates (i, j, k) is recorded by setting
the polarization variable P(i,j,k) equal to 0 and a time variable t(i,j,k) equal to
1. In this model, the depolarization is isotropic so there is no preferred direction
for the depolarization wave to travel. In the next time step, any of the 26 neighbors
of the original cell that have not yet been depolarized (P(il,jl,kl) = 1) may be
depolarized. Thus in general, a cell C with coordinates (il,jl,kl) can be depolarized
at a time step t if it is not yet depolarized and if its neighbor C1 with coordinates
(i,j, k) was depolarized in the previous time step. Once a cell is depolarized, k)
is incremented by one for each time step after depolarization. When t{i,j,k) reaches
the repolarization interval (r), the polarization variable P(i,j,k) is reset to 1 and
the cell can be depolarized again. During depolarization, current is assumed to flow
instantly from the cell depolarized in the last time step to the cell being depolarized
in the current time step.
The model also allows for automaticity of ectopic cells located anywhere in the
heart. The depolarization algorithm for the firing of these cells uses the celltocell
propagation described above, but the depolarization wave spreads from the ectopic
site instead of from the six targeted cells in the septum. Also, the Purkinje fibers are
only depolarized through cell to cell propagation and are not manually activated four
16
time steps after the firing of the ectopic focus.
The model also allows for a combination of ectopic and SA firing so that waves
can spread from two different starting points. In this case each of the waves will
depolarize part of the heart and meet somewhere between the two starting points.
The waves can overlap only if a cell depolarized by one wave has repolarized before
being hit by the second wave. If a cell is already depolarized when it is hit by another
depolarization wave, the wave will not spread to this cell. This interaction of the two
waves is what causes many of the irregular heart beats, arrythmias and tachycardias
that will be discussed in Chapter 5.
In several CA models, the anisotropic nature of the heart muscles is modeled
[9] [10] [11] [15], This is done by giving cells a preferred direction along the
muscle fiber so that the depolarization wave will spread faster in that direction than
across the muscle. The MillerGeselowitz model [9][10] approximates the anisotropic
properties of the cells in each layer from endocardium to epicardium by using a block
of anisotropic tissue. In this block there is an axis of low conductivity oriented normal
to the epicardium which represents the slow flow of charge across muscle fiber. In
all directions parallel to the epicardial surface there are axes of high conductivity
representing the fast flow of charge along the muscle fiber. Thus a wave starting at a
point within this anisotropic tissue spreads out in an elliptical wave front, depolarizing
more of its neighbors within the same muscle fiber than in adjacent muscle fibers
(Figure 3.2). The magnitudes of the slow (vt) and the fast (vp) velocities for the
elliptical wave are determined using
Vp
Vt
L + L
9it Set
L + L
/ 1 n___
(3.4)
where git,get,9ip, and g^p are constants of intercellular and extracellular conductivity
determined from experiments on the cardiac tissue. Using these values, the isochrones
of the depolarization wave were found to match those of a working heart. The values
of vt and vv are constant throughout the heart.
17
Figure 3.2 Elliptical wavefronts on a cross section of a modified MillerGeselowitz model. A
sample ellipse of activation is shown with the cods of slow and fast propagation flOJ.
In the atrial and ventricular Japanese heart model [11], anisotropy is also used in
the depolarization algorithm. The method of determining the anisotropy in this model
involves layering the heart like an onion and then twisting each layer a specific angle
from apex to base. The elliptical wave fronts of the depolarization wave found with
this model are similar to those of the MillerGeselowitz model. When comparing the
global excitation sequences of isotropic and anisotropic models, the authors of the
Japanese model found that there were no significant differences. This is attributed to
both the high conduction velocities of the isotropic Purkinje fibers and to an averaging
effect of the rotation of the fibers over the entire heart. The anisotropic heart did
produce significant improvement over the isotropic heart in various experiments such
as initiation of ventricular fibrillation.
The algorithm used for incorporating anisotropy in either the MillerGeselowitz
or the Japanese model could be added into our model in a similar manner. The main
advantage of incorporating anisotropy into our model would be to obtain a closer
correlation between actual and simulated ECGs in the cardiac experiments which are
discussed in Chapter 5.
3.2 Repolarization Algorithm
As explained in Section 2.2, the action potential of a myocardial cell has five
18
main phases. The third phase, rapid repolarization, is primarily responsible for the
Twave on the ECG because most of the energy transfer between the interior and
exterior of the cell occurs in this phase. However, the repolarization of a cell is not
instantaneous and in fact occurs over a time period which is about two to three times
as long as the duration of the Rwave. The Twave also has an absolute height of
one third of the height of the Rwave.
In our model, the repolarization wave is based on a parameter D which is the
number of time steps required to depolarize all the cells in the heart; usually D is
8 to 10 time steps which corresponds to 80 to 100 ms. Since the repolarization
of a cell is gradual, in the model it is necessary to make each cell repolarize over
several time steps instead of instantly. In order for the duration and height of the
Twave to correspond correctly to that of the Rwave, the time needed for a cell
to completely repolarize should be the same as the time required to depolarize the
entire heart. Thus in each time step, a cell will become jj repolarized and it will
take D time steps for complete repolarization. If a cell depolarizes at time step t then
it will begin repolarizing after the repolarization interval r or at time t + r. The
repolarization interval corresponds to phase 1 through the beginning of phase 3 of
the action potential and is usually 80 to 250 ms long. This interval represents the
absolute refractory period and is set by the user. In the next section which deals with
the Twave, the repolarization interval of each cell will depend both on the user set
interval and on the location of the cell in the heart.
Once a cell begins repolarizing, a repolarization parameter u is used to keep
track of how long the cell has been repolarizing and when u reaches 1 + D, the
cell is totally repolarized. The value of u corresponds to the potential of the cell as
was discussed in Chapter 2; so that u = 1 corresponds to a potential of 0 mV and
u = 10 corresponds to a potential of 90 mV. The height of the repolarization wave
is recorded in a repolarization array hr(t) at each time step t. Thus, if a target cell
repolarizes in the current time step, the contribution of this cell to the repolarization
array is found by analyzing the potential of all of its neighbors. If a neighbor of the
target cell repolarized in the previous time step then the potential of the neighboring
19
cell will be 10 mV lower than that of the target cell. In order to simulate a IIlead ECG
this change in potential is projected along the zaxis of the heart, so the contribution
to the repolarization array is .10 times the difference in the zcomponents of the
target cell and its neighbor cell. Neighbors of the target cell that did not repolarize
in the previous time step do not contribute to the repolarization array. The total
contribution of the target cell to the repolarization array in the current time step is
found by computing the contribution for all of its neighboring cells. The difference
in repolarization states between the target cell and its neighboring cells remains the
same during the entire time necessary to repolarize a cell which is equal to D or about
10 time steps. Thus the total contribution of each cell is added into the repolarization
array for 10 time steps after beginning of repolarization which occurs at the end of
the repolarization interval for that cell. In Figure 3.3 a schematic diagram of the
repolarization of the cell (10,11,4) at time step t = 34 is shown. The contribution
to the repolarization wave from each of the neighbors of this cell in the x z plane
is described.
Figure 3.3 Schematic representation of repolarization of a target cell located at (10,11,4).
The u values which tell how long a cell has been repolarizing are shown on the grid for each
cell. Each of the cells above the target cell began repolarizing one time step before the target
cell, so the charge contribution for each is. 1(4 5) = .1.
Using this algorithm with D set to 10, the Twave will be three times as long as
the Rwave for a given depolarization speed. However, using this algorithm without
allowing for variation of the repolarization intervals produces a Twave which is of
20
the correct duration, but is negative instead of positive. Many experiments have tried
to resolve the fact that the Twave has the wrong sign and some of these will be
explored in the next chapter.
The repolarization algorithms used in the MillerGeselowitz [9][10] and Japanese
[11] models are similar to this repolarization algorithm. However, as was described
in Section 3.2, the Japanese model works with an anisotropic heart which is layered
from epicardium to endocardium. In this model, the repolarization intervals for each
cell are found by first determining the general repolarization interval for the entire
layer and then varying it slightly for each individual cell based on its position in the
layer. A modification of this technique is applied to our model in Chapter 4.
3.4 Output
In most CA models, the output is in the form of an ECG which is an electronic
record of the activity of the heart as described in Chapter 2. In our model, only the
II Lead ECG is used although the other major leads can easily be implemented by
projections. The depolarization wave values formed using the algorithm discussed in
Section 3.2 are computed for each neighbor, C, of C' by forming a vector C' C
= (i il.j jl,k kl). The sum of the projection of these vectors onto the zaxis
which runs parallel to the II Lead is recorded as the pulse height at each time step t.
Since this is an integer grid and since any neighbor of C' must be in the unit cube
centered at C', the depolarization of each neighbor of C' adds 1 (the distance in the
z direction) to the pulse height which is recorded in the depolarization array, hd(t).
In the repolarization algorithm, the portion of repolarization that occurs at each cell
during one time step is recorded in the repolarization array, hr(t). Each cell con
tributes A = array jn any time step. The total pulse height at any time t
is the sum of the depolarization and repolarization heights and is recorded in a third
array h(t) = hr(t) + hd(t). This pulse height is then graphed with the corresponding
time step to form the QRST complex of the ECG scan. Since there are no atria, there
are no Pwaves in the ECG output of the model.
21
4. Experiments on the Twave
In a normal ECG recording, the Twave, which corresponds to the repolarization
of the heart is upright, positive and has about one third the height and three times
the length of the Rwave (Figure 4.1). The Rwave represents the depolarization of
the ventricles in which the potential of the individual cells changes from 90 mV
to 10 mV. Since the depolarization wave travels down the ventricles, the cells in the
top of the heart will be depolarized first causing the potential to be positive there
while the nondepolarized cells in the bottom of the heart will cause the potential
to be negative there. This difference accounts for the positive Rwave in the II
Lead ECG which measures the difference in potential along the zaxis or between
the top and the bottom of the heart. If the repolarization wave is also assumed to
travel down the heart, then the Twave which represents this wave will be negative.
This is because during repolarization, the individual cells change from a potential of
10 mV to 90 mV. Thus the potential top of the heart where the cells have begun
repolarization and have a potential of 0 to 90 mV, will be more negative than in
the bottom of the heart where the cells have not begun repolarization and still have
a potential of 10 mV. Thus the positive polarity of both the Rwave and the Twave
seems to indicate that the repolarization wave does not travel down the heart like the
depolarization wave does. In this section, four different modifications of the model
which explore the spread of the repolarization wave and additional features of the
myocardial cells will be used to attempt to explain the positive Twave.
22
Figure 4,1 ECG for a normal SA pulse without Pwaves.
4.1 Model A
The first model is the nonmodified model described in Chapter 3. As noted
before, if the repolarization intervals of all cells are the same, then the Twave of a
normal SA pulse is inverted, as shown in Figure 4.2. This is because the first cells
to repolarize will be those that depolarized first which are the six targeted cells in
the septum. The repolarization wave, although prolonged, will spread in the same
manner as the depolarization wave which is essentially down the septum from the
endocardium to the epicardium. Since the wave spreads predominantly down the
heart, at any given time step, most cells will have more cells with larger negative
charge (more repolarized) above them than below them. Thus, when the II lead ECG
is formed from the difference in potential between the top and the bottom of the heart,
the Twave will be negative. Figure 3.3 shows a cell with coordinates (10,11,4) and
the potential of its neighbors at time step t = 34. The contribution of this cell to the
height of the repolarization array hr (34) is 0.3 mV. In general, the total contribution
of all the cells in a given time step will be negative, leading to the inverted Twave.
Thus, using the same repolarization interval for all of the cells does not provide an
accurate ECG output.
23
Figure 4.2 ECG of an SA pulse for Model A The Twave is of correct duration and size, but
it is inverted. The minimum value of the Twave is 18 mV
4.2 Model B
In several cardiology texts [3], it is suggested that the positive Twave is formed
because of the rule that if a pulse originates from the SA node, then the first cells to
depolarize will be the last cells to repolarize. Although this rule will account for the
positive Twave, it is presented without an explanation in terms of the ionic processes
of depolarization and repolarization. A cell in the bottom of the heart depolarizes
because of an influx of sodium. In order to distinguish whether the depolarization
wave affecting this cell originated at the SA node or from an ectopic focus, there
must be a difference in either the amount of sodium in the extracellular fluid or in the
amount of sodium allowed to pass into the cell based on the origination site. Since
the CA does not include extracellular fluid or the ionic process, only the end results
of this process can be incorporated into the model. The end result of the rule is that
if the pulse originates from the SA node, the entries in the repolarization array which
form the Twave should be positive. In the model, if a pulse originates in the SA
node, the maximum height of the depolarization array hd(t) is generally much higher
than if the pulse originates from an ectopic focus. Thus in order to implement the rule
into the model, a trigger (s2) is set to 1 if the maximum value of the depolarization
array is in the range of that of an SA pulse. Then, when the output of the model is
24
formed, the entries of the repolarization array are:
hr(t + j) = .06hd(t r + 1) + hr(t + j) if s2 = 1 (pulse from SA node)
hr(t + j) = .llhd(t r + l) + hr(t + j) if s2 = 0 (pulse from ectopic focus)
This produces a positive Twave if the pulse originates from the SA node and a
negative Twave if the pulse originates from an ectopic focus. The values of .06 and
.11 were chosen so as to make the output match that of a normal heart [13], The
ECG produced for a normal SA pulse in this model is the same as that of a working
heart (Figure 4.1).
The main problem with this model is in identifying an SA pulse. It is possible for
some ectopic pulses to produce a maximum value of the depolarization array that is
in the normal range for an SA pulse. If this occurs, the output for the model will be
flawed. Also, the lack of biological support for this repolarization algorithm reduces
the predictive capability of this model. In the next two models, the biological basis
for the inverted Twave will be explored.
4.3 Model C
Several studies on canine and human hearts have found that in fact the repo
larization intervals are not the same for all the cells in the heart [4], In a majority
of the studies it has been found that the repolarization interval of ventricular cells is
longer in the endocardium than in the epicardium. This repolarization variation was
included in Model C by layering the ventricles into parabolic shells representing the
epicardium (outer shell), middle section, endocardium (inner shell) and septum. The
repolarization intervals were then made the shortest in the epicardium and the septum
(r 4 time steps) where r is the normal repolarization interval, slightly longer in
the middle shell (r 2 time steps) and the longest in the endocardium (r). These
intervals were picked after experimenting with several different interval times and
finding that these worked as effectively as any other variation. Implementing the
activation sequence with this repolarization method, the ECG for a standard SA pulse
is shown in Figure 4.3.
25
Figure 4.3 ECG for a normal SA pulse from Model C. The Twave is of correct duration, but
it is inverted. The minimum value of the Twave is 10 mV which is 8 mV higher than that
produced in Model A
Although the Twave is still negative in this output, it is less negative and shorter
than the Twave for model A (Figure 4.2). So, varying the repolarization intervals
from epicardium to endocardium does affect the Twave, but it does not appear to be
enough to make the Twave positive. Thus, it seems that there must be other factors
that affect the repolarization process.
4.4 Model D
In a study that considers the control of the action potential and how it relates
to ionic features, Nobel [19] hypothesized that in order to obtain a positive Twave,
the repolarization intervals should be shorter in the bottom than in the top of the
ventricles. In order to implement this effect in Model D, repolarization intervals for
each cell are based on epicardium to endocardium position as in Model C and then
this value is added to  where z is the zcoordinate of the cell and ranges in
value from 1 to 15. Thus if a cell is located in the epicardium and has coordinates
(4,5,9) then,
z 9
2 2
= 4 and its repolarization interval is (r 4 + 4) where
r is the general repolarization interval set by the user. The value
is based on
26
experimentation with model output. The repolarization intervals for an x z cross
section of the model are shown in Figure 4.4.
30 28 28 28 28 28
24 28 28 28 28 28
27 31 27 27 27 27 27
27 31 27 27 27 27 27
26 28 30 26 26 26 26 26
26 26 30 26 26 26 26 26
25 27 29 25 25 25 25 25
25 25 29
30 28 28
28 28 28
31 27 27 27
31 27 27 27
30 28 26 26 26
30 26 26 26 26
29 27 25 25 25
29 25 25 25 25
25 25 25 25
24 24(g) 28 24 24 24 2T24 28 26 24 24 24
24 24 28 24 24 24 24 24 28 24 24 24 24
23 23 25 27 27 27 23 27 27 27 25 23 23 23
23 23 27 27 27 23 27 27 27 23 23 23 23
22 24 26 26 22 26 26 24 22 22 22
22 24 24 24 24 24 22 22 22
21 21 23 21 21 21 21
This cell is in the septum with a zco ordinate of 8. The
repolarization interval is (user set4+jjJrJ)=(254+4)=25.
Q This cell is in the middle layer with a zcoordinate of 7. The
repolarization interval isi (user set2+ y )=(252+3)=26.
Figure 4.4 Eepolarization intervals on an xz cross section of Model D. The intervals increase
from bottom to top and from epicardium to endocardium. The value of each interval = (user
set value layer value + where the layer value is.epicardium and septum=4, middle
layer=2, endocardium=0.
The ECG for a standard SA pulse for model D is shown in Figure 4.5. In
this output, the Twave is of the correct height, length and direction and the ECG
corresponds very well with that of an actual SA pulse (Figure 4.1). However, when an
ectopic pulse is used in this model, the Twave remains positive instead of inverting
as it would in a working heart. Thus, this model works for a normal pulse but does
not explain the formation of the Twave for all possible activation sequences in the
heart.
27
Figure 4.5 ECG for a normal SA pulse from Model D. The Twave is of correct duration and
size, and it is positive
Of the four models, Model B and Model D both produce a correct ECG output
for the standard SA pulse. However, while the repolarization algorithm in Model B
is merely contrived in order to make the Twave work, the hypothesis used to form
the repolarization algorithm for Model D has a biological basis. The results of these
experiments with the Twave indicate that further research should be conducted on
the theories of apex to base repolarization interval variation and on determining the
biological basis for the rule used in model B.
28
5. Experiments
As cardiac computer models become more complex and more accurate in simulat
ing the functions of the heart, their use in cardiology is also increasing. A computer
model can be used in cardiology in two ways. The first method is to use the model to
evaluate a complex hypothesis about some biological phenomenon. By running the
computer simulation several times it is possible to support a hypothesis or perhaps
find cases in which the hypothesis cannot hold. The second method is to use the
computer model to predict what will happen in a particular biological situation. This
approach relies on the accuracy of the computer model not only in simulating what is
known, but in containing enough biologically significant properties to make accurate
predictions. Cellular automaton models have been used in both of these approaches
with substantial success. This chapter will look at several such experiments using the
models described in the previous chapters.
5.1 SelfSustained Ventricular Tachycardia (SSVT)
Self sustained ventricular tachycardia (SSVT) can be caused when there is some
type of block in the cells, such as scar tissue, that causes the depolarization wave to
change its normal course. If the geometry of the block is right, the depolarization wave
can actually circle around this blockage. This happens because cells that depolarize
on one side of the blockage have completely repolarized by the time the depolarization
wave hits them from the other side of the blockage [3], The blockage causing SSVT
does not have to be damaged tissue. If a group of cells are depolarized slightly before
the depolarization wave from the SA node arrives, then these cells form a blockage
because they cannot be depolarized again at that time. In the model, a single ectopic
pulse applied prior to a single SA pulse can cause such a blockage and initiate SSVT.
In Model B, a blockage occurs when the ectopic focus at coordinates (16,11,9),
29
which is in the interior of the left ventricle, is pulsed once at time step 12 and the SA
node is pulsed once at time step 19 with a repolarization interval of 120 ms (Figure
5.1). The heart beat is irregular until the SA pulse and then it settles down into a
regular pattern of ventricular tachycardia which continues without any stimulus from
the SA node or an ectopic focus. In models A and C which have an inverted Twave
for the SA pulse, SSVT can be initiated in a similar manner. This suggests that even
though the Twave for SA pulses is incorrect in these models, they are still effective
in modeling essential features of the activation process.
Figure 5.1 SSVT for Model B one ectopic and one SA pidse. The pulse is erratic for the first
600 ms and then it settles into steady VT without any additional pulsing of either the SA node
or the ectopic focus.
In Model D the apex repolarization interval must be set to less than 120 ms in
order to obtain a selfsustaining rhythm. The rhythm observed here is much more
irregular than in the other models. Although the tachycardia continues without any
additional stimulus, it does not seem to settle into a regular rhythm (Figure 5.2). This
suggests that the blocking factor in the interaction of the pulse from the ectopic focus
and the SA node shifts around the ventricle instead of being confined to one region.
Thus, the high variability of the repolarization interval seems to make it more difficult
to stimulate the rhythmic SSVT behavior. The fact that Model D does not include
the anisotropy of the myocardial fibers may also be a reason for the difficulties in
obtaining rhythmic SSVT beats. In the Japanese model [11] introduced in Chapter
3, the authors found that in the anisotropic model it was easier to initiate ventricular
30
tachycardia and ventricular fibrillation than in the isotropic model.
Figure 5.2 Nonrhythmic SSVT for Model D. Using the same parameters as in the SSVT for
Model B, Model D does not settle down into a rhythmic VT as Model B does.
In all of the models it was possible to stop the SSVT and return the heart to a
normal sinus rhythm using capture beats. A capture beat is either an external shock
or rhythmic firing of the SA node that disrupts and stops abnormal beats. Using the
ectopic focus and SA node interaction that produced SSVT in the previous example,
a capture beat occurs when the SA node pulses rhythmically eveiy 860 ms instead
of just once. In the ECG shown in Figure 5.3, the first SA firing, which occurs 90
ms after the firing of the ectopic focus causes SSVT. The second SA firing occurring
at 1000 ms disrupts the SSVT and the tachycardia ends at 1500 ms. The third SA
firing occurring at 1860 ms is normal, showing that the SSVT has been captured and
the heart is beating normally. This demonstrates that SSVT caused by the firing of
a single ectopic pulse can be stopped by external (or SA node) impulses occurring
at the correct time. However, application of these external impulses at the wrong
time can further disrupt the beating process making the SSVT worse, sometimes
causing ventricular fibrillation which is often fatal. The use of CA models allows for
experimentation in order to understand when an external pulse should be applied.
31
Figure 5.3 ECG of SSVTfor Model B with a capture beat occuring at I860 ms.
5.2 Ventricular Tachycardia
In addition to the abnormal rhythms that can be caused by just one ectopic pulse,
ventricular tachycardia (VT) can also be caused by a constant rhythmic pulsing of an
ectopic focus. In this type of VT, external pulses are still effective in arresting the
abnormal beating, but do not stop the automaticity of the ectopic focus. Thus, the
external pulse only interrupts the tachycardia which can start again if the ectopic focus
resumes its automaticity. In order to prevent recurrent tachycardia, the automaticity
of the ectopic focus must be permanently eliminated. This is done through a process
called ablation in which the focus is located by probing the tissue of the heart with
an electronic pulse from a surgically inserted catheter which initiates the automaticity
of the abnormal focus when it is probed. Once the abnormal focus is found it is
cauterized and thus can no longer fire. This procedure is effective in stopping the
automaticity, but it is both timeconsuming and technically difficult to locate the focus
through probing methods. A CA model can be used to help determine the location of
such a focus by experimentally choosing foci in the CA model and then seeing how
the ECG output from each focus matches the ECG of the patient.
To show that this type of differentiation is possible, VT is initiated from the
32
upper left ventricle at (16,11,11) (Figure 5.4) and from the lower right ventricle at
(5,11,6) (Figure 5.5). These ECG scans are significantly different, showing the ability
of a CA model to determine where to search for an ectopic focus. In addition, it
is possible that specific modifications could be made to the CA model to produce
greater differentiability in the ECG scans of VT.
Figure 5.4 Ventricular tachycardia initiated in the upper left ventricle. The Rwave has a
maximum value of 102 mV and the Twave has a minimum value of 70 mV.
Figure 5.5 Ventricular tachycardia initiated in the lower right ventricle. The Rwave has a
maximum value of 52 mV and the Twave has a minimum value of101 mV. The Rwave is
more rounded than in figure 5.4.
33
5.3 WolffFarkinsonWhite Syndrome (WPW)
Another area where CA models can be used to increase the precision of cardiac
surgery is in correcting WolffParkinsonWhite Syndrome (WPW). WPW is caused
by conduction from the atria to the ventricles through an accessory pathway in ad
dition to the conduction through the AV node. The abnormal depolarization cycle
in the ventricles can lead to various disorders such as supraventricular tachycardia
or ventricular tachycardia. WPW is very distinguishable on the ECG scan because
the depolarization through the accessory pathway produces a delta wave which is a
slurring in the initial portion of the QRS wave [3], There are eight main areas in
the base of the ventricles where these accessory pathways are usually formed [2]
(Figure 5.6). As with automaticity of an ectopic focus, ablation is the main technique
for eliminating these pathways. The ablation of the accessory pathway forces the
depolarization wave to travel only through the AV node and thus returns the heart to
normal beating.
Figure 5.6 Lateral slice across the top of the ventricles showing the eight main areas for
accessory pathways in WPW.
In a study using a modified MillerGeselowitz model [10], WPW was successfully
simulated in each of the eight sites shown in Figure 5.6. Since this model has no atria,
the accessory pathway of WPW was simulated by initiating an ectopic pulse on the
34
base of the ventricles in each of these sites. The ectopic pulse was fired approximately
40 ms before the start of normal excitation from the endocardium. The Lead I, II,
and HI ECGs for the simulation from site 8 is shown in Figure 5.7. The polarities of
the delta waves for all eight sites were then compared with each other and with actual
patient data. The results of this experiment showed that the model was reasonably
effective in predicting the site of the accessory pathway although the ECGs for some
sites were closer to those of actual patients than other sites. This discrepancy can
be explained partially by the sensitivity of the delta wave to the exact initiation site.
The fact that the model has only 50,000 cells may make the propagation too coarse
to reproduce the delta wave polarity accurately in all cases. Improvements both in
the precision and in the depolarization algorithms of the model may help to eliminate
these discrepancies.
Figure 5.7 Lead I, II, and III ECGs for WPW in a modified MillerGeselowitz model. The
Lead II and III ECGs show delta waves in the QRS complex [10].
5.4 Right and Left Bundle Branch Block
In addition to being able to simulate ectopic foci and accessory pathways, another
feature of CA models is the ability to simulate ischemia or damaged tissue and
conduction system disorders. One very common type of conduction system disorder
is right bundle branch block which occurs when damage to the conduction system
stops the depolarization wave from traveling down the right bundle branch to the
Purkinje fibers. This means that the right ventricle will be depolarized after the left
since the depolarization wave cannot reach the right ventricle directly, but must first
35
travel through the bottom of the heart. This disorder is recognizable on the ECG by
a prolonged QRS complex which indicates the additional time required to depolarize
the right ventricle. Since this is a depolarization disorder, and it is not affected by the
repolarization sequence, Models A, B, C, and D are all able to simulate this disorder
correctly. In Figure 5.8 the normal SA pulse ECG for model D is shown along with
the right bundle branch block ECG which has a QRS that is approximately 10 ms
longer than in the normal ECG. The CA model can also be used to simulate other
ischemic disorders such as scar tissue caused by a heart attack or other conduction
system blocks [12].
Figure 5.8 ECG for simulated right bundle branch block and a normal SA pulse in Model D.
The QRS of the right bundle branch block ECG is 10 ms longer than for the SA pulse.
36
6. Conclusions
The CA models explored in this thesis are effective in modeling many of the
functions of the heart. The models can be used to explore the interaction of pulses
originating at an ectopic focus, producing an ECG output that simulates that of an
actual heart. The models correctly reproduce ECGs for the depolarization of the
heart and provide insight into the repolarization process. By using a CA model it is
possible to understand exactly how the depolarization and repolarization waves are
formed based on the travel of action potentials through the cells. This understanding
provides a tool for improving analysis of ECG output.
The experiments on the inverted Twave show how CA models can be used to
support and encourage research in this area. In Chapter 5, two of the models pre
sented suggest specific approaches to analyzing the repolarization intervals in order
to determine why the Twave is inverted. Resolving the Twave question for CA
models would also answer many other unresolved questions about anatomy and cel
lular interactions. If the CA produced a working Twave based on some set of logical
rules, then it is possible that the biology supporting these rules may be applicable to
the heart.
The use of CA models as a diagnostic tool in disorders such as selfsustained ven
tricular tachycardia, ventricular tachycardia, WolffParkinsonWhite syndrome, right
bundle branch block and others may prove to be very useful. The CA can be a
powerful tool in improving surgical procedures by eliminating error and increasing
success in finding the source of disorders through simulations and comparisons of
model and actual ECG output. It is possible that the user interface of this model
could be improved so as to make it a viable model for doctors to use in diagnosis
and correction of cardiac disease. In addition, the CA model is a valuable tool for
learning and experimenting with the activation sequence and cellular interactions in
37
the heart.
There are many areas in which the CA models can be improved. The accuracy of
the activation sequence can be increased by including more of the essential character
istics of the action potential such as a relative refractory period. The number of cells
in the CA model can be increased to produce a finer resolution on the ECG output,
and provide more variability in modeling applications. Finally, with the increasing
power of computers it is possible to combine the ionic and CA models to produce
an accurate and efficient model which can precisely model the ionic interactions on
a cellular level.
38
7. References
1. G. Beeler, H. Reuter, Reconstruction of ventricular myocardial fibers, Journal of
Physiology, Vol. 268, (1977), pp. 177210.
2. D. Benson D, J. Gallagher, M. Spach, R. Sterba, Localization of the site of
ventricular preexcitation with body surface maps in patients with WolffParkinson
White Syndrome, Circulation, Vol. 65, No.6, (1982), pp. 12591268.
3. R. Berne, M. Levy, Cardiovascular Physiology, Mosby Year Book, (1992).
4. M. Burgess, P. Ershler, K. Spitzer, B. Steinhaus, Nonuniform epicardial activation
and Repolarization properties of in vivo canine pulmonary conus, Circulation
Research, Vol. 62, No. 2, (1988), pp. 233246.
5. A. Camm, M. Malik, Cardiac electrophysiological experiments in numero, parts
/////, PACE, Vol. 14, (1991), pp. 16481671 and 21672186.
6. A. Camm, M. Malik, Computer simulation of electronic interactions during exci
tation and repolarization of myocardial tissue, Medical and Biological Engineer
ing, (1991), pp. 425432.
7. D. Davis, Differential Diagnosis of Arrhythmias, WB Saunders Co, (1991).
8. A. Despopoulos, S. Silbemagl, Color Atlas of Physiology, Thieme Medical Pub
lishers, (1991).
9. D. Geselowitz, W. Miller, Simulation studies of the electrocardiogram. /. The
normal heart, Circulation Research, Vol. 43, (1978), pp. 301315.
10. R. Gulrajani, M. Lorange, Computer simulation of the WolffParkinsonWhite
preexcitation syndrome with a modified MillerGeselowitz heart model, EEEE
Transactions on Biomedical Engineering, Vol. BME33, No. 9, (1986), pp. 852
873.
39
11. E. Harasawa, K. Harumi, H. Hosaka, D. Okazaki, D. Wei, Comparative simulation
of excitation and body surface electrocardiogram with isotropic and anisotropic
computer heart models, IEEE Transactions on Biomedical Engineering, Vol. 42,
No. 4, (1995), pp. 343357.
12. K. Harumi, T. Musha, G. Nishiyama, Y Okamoto, H. Tusunakawa, D. Wei, G.
Yamada, Clinical application of electrocardiographic computer model, Journal of
Electrocardiology, Vol. 22 Supp., (1989), pp. 5436.
13. H. Hermes, A computer model for analysis and control of ventricular arrhythmias,
Unpublished, (1995).
14. A. Hodgkin, A. Huxley, A quantitative description of membrane current and its
applications to conduction and excitation in nerve, Journal of Physiology, Vol.
117, (1952), pp. 500544.
15. B. Horacek, L. Leon, Computer model of excitation and recovery in the anisotropic
myocardium, Journal of Electrocardiology, Vol. 24, (1991), pp. 115.
16. M. Kawato, K. Okazaki, R. Suzuki, S. Urushibara, A. Yamanaka, Feconstruction
of electrocardiogram using ionic current models for heart muscles, Japanese Heart
Journal, Vol. 27 Supp., (1986), pp. 185193.
17. R. McAllister, D. Nobel, R. Tsien, Feconstruction of electrical activity of cardiac
Purkinje fibers, Journal of Physiology, Vol. 251, (1975), pp. 159.
18. D. Nobel, A modification of the HodgkinHuxley equations applicable to Purkinje
fiber action and pacemaker potentials, Journal of Physiology, Vol. 160, (1962),
pp. 317352.
19. D. Nobel, Ionic mechanisms controlling the action potential duration and the
timing of repolarization, Japanese Heart Journal, Vol. 27 Supp., (1968), pp.
319.
20. J. Van der Mark, B. Van der Pol, The heartbeat considered as a relaxation oscil
lator and an electrical model of the heart, Philos Mag, (1928), pp. 763775.
40
grids. Chang, Clark and Hare [4] use the algorithm of Hare, Hedetniemi and
Hare to generate minimum domination patterns of 5 x r and 6. x r grids, and
what they conjectured to be minimum dominations of 7 x r, 8 x r, 9 x r and
10 x r grids for all r. Theorem 3.12 shows these dominations are minimum.
They also found formulas for 7nir when r < 122 and for 712,r when r < 33 (the
formula for 711 _r holds for all r, while the formula for 712,r does not).
The program in this paper was modified to store backtracking information
to give minimum dominations of 11 x r, 12 x r, 13 x r, 14 x r, and 15 x r grids
for all r up to a certain point. These were used to find minimum dominations
for n x r grids for a given n and for all r (this is actually quite difficult because
only one of what is usually many minimum dominations is generated). Chang,
Fisher and Hare [5] report these minimum dominations.
For n,r > 8, Chang [2] gives a construction which dominates an nxr grid
with [(n + 2)(r + 2)/5j 4 vertices. This is a generalization of a construction
of Cockayne, Hare, Hedetniemi and Wimer [7] for r x r grids. Theorem 3.12
shows this construction is minimum for 16 x r, 17 x r, 18 x r, and 19 x r grids
when r > 16. This causes one to suspect the following conjecture (originally
due to Chang [2]) might be true.
Conjecture 3.13 For all n,r > 16,
7(P x Pr) =
(n + 2)(r + 2)
5
38
If Conjecture 1 is true, dominating an n x r grid would be analogous
to 2packing (Fisher [8]j in that a simple formula holds for n and r large
enough. Further, it would answer a question of Hedetniemi and Laskar [12].
In commenting on [6], they asked is there a polynomial algorithm for find
ing 7(Pn x Pr)? Conjecture 3.13 together with Theorem 3.12 would give an
affirmative answer to the question.
3.4 Domination of G x Pr: The At Least Algorithm
We will now use the At Least Algorithm to find the domination number
of G x Pr where G is any graph on n < 6 vertices.
Livingstone and Stout [14] develop a dynamic programming algorithm for
finding 7(G x Pr) using an algorithm that is similar to the Exact Algorithm.
We are able to significantly reduce the complexity of this algorithm by stating
in terms of the At Least Algorithm.
To illustrate the At Least Algorithm for 7{G x Pr) we will let G be the
house graph on five vertices as shown in Figure 3.3.
39
Figure 3.3 The house graph on five vertices G and G x P4.
There are 77 possible states for the domination of the house graph on five
vertices. These states are shown below and the vertices are ordered v4... v5 as
labeled in Figure 3.3.
Si : = 00000 521 = 01221 541 = 11121 $61 = 21111
52 '  00001 S22 = 10000 542 = 11122 562 = 21112
S3 ' = 00010 S23 = 10001 543 = 11210 563 = 21121
S4 = = 00011 S24 ~ 10010 S44 = 11211 564 = 21122
s5  = 00100 525 = 10011 s45 = 11212 565 = 21211
s6  =.00101 526 = 10101 546 = 11221 566 = 21212
s7  = 00110 527 = 10 1 01 S47 = 11222 567 = 21221
Ss  = 00111 528 = 10110 548 = 12101 568 = 21222
s 9 = = 00121 529 = 10111 549 = 12111 569 = 22101
SlO = 01000 530 = 10121 550 = 12112 570 = 22111
Su = 01001 531 = 11000 551 = 12121 571 = 22112
512 = 01010 532 = 11001 552 = 12122 572 = 22121
5l3 = 01011 533 = 11010 553 = 12211 573 = 22122
Si4 = 01100 534 = 11011 s54 = 12212 574 = 22211
5l5 = 01101 535 = 11012 555 = 12221 575 = 22212
516 = 01110 536 = 11100 5 56 = 12222 576 = 22221
Sl7 = 01111 537 = 11101 s57 = 21001 577 = 22222
5l8 = 01121 538 = 11110 558 = 21011
Sl9 = 01210 539 = 11111 5 sg = 21012
S20 = 01211 540 = 11112 560 = 21101
These states form a partial ordering were s; follows Sj in the ordering if
Si is at least Sj which means the labelings sj(p) < 5;(p) for the corresponding
40
vertices p in each state. The partial ordering is used to form the At Least
Matrix A and the Minimum Previous Matrix B as described for Pn x Pr.
Next, we form the states vector tr where tr(st) is the number of vertices
needed to dominate G x Pr and leave the last copy of G in at least state s, and
the ending vector y which has a 0 for the state 1 and infinity otherwise. For
this example, 1 is state S39. The vectors for the house graph are shown below.
y =
t3 =00
t2 28
t3 =38
[00000000000000000000000000000000]
[21111110111111110010011111100011111110000000011100010112100001010112011121223]
[11111111111111111111010101000000000000000000001001010112000000010112001121222]
tis =17 0
tig = 18 0
t20 = 19 0
t2i = 20 0
t22 = 21 0
t23 = 22 0
[00000000000000000000100000000000000000001101112001121223001011121222112232323]
[11000000010000000000100000000000000000001101112001121223001011121223112232334]
[11000000010000000000100000000000000000001101112001121223001011121222112232323]
[11000000010000000000100000000000000000001101112001121223001011121222112232323]
[11000000010000000000100000000000000000001101112001121223001011121222112232323]
[11000000010000000000100000000000000000001101112001121223001011121222112232323]
Notice that t2o = 1 0 tig and t21 = 1 0 t2o Thus after 20 iterations,
periodicity occurs and we have tr = 1 0 tr_i for all r > 20. This allows us to
find a closed form for the domination number of the house graph G. For all
r > 1, we have
I 1 if r =1, 2, 3, 4, 6, 7, 8, 11, 12, 16
7(
Now to find the domination number of G x fi, first note that if G is
disconnected then so is G x Pr. Thus, if G\, G2,... are the components of G,
then the domination number of G is just the sum of the domination numbers of
41
its components. Thus we only need to find the domination number of connected
graphs. To further reduce the number of graphs we need to consider, the
following theorem gives the domination number for certain graphs with high
degree.
Theorem 3.14 If a graph on n vertices G has a set of vertices ai Ã‚Â£ A with
deg(ai) > n 3 and for each a, Ã‚Â£ A, all vertices not in the closed neighborhood
of ai are also in A, then 7(G x Pr) r or r A 1 or r A 2.
Proof: Let a1,a2) be the vertices in set A as described above. To
dominate (G x Pr) start with the k~ copy of G = Gk where k ^ l,r. Place
aik in the dominating set. Now there are at most two undominated vertices
in G both of which are also in A, call these vertices ajk and a^. Place ajk1
and aik+i in the dominating set. Now ahki is adjacent to a^k by definition
of G x Pr, so all the vertices in Gk are now dominated. Next we have that
in Gk1, aiki and the closed neighborhood of ajk1 are dominated and only
one vertex aqk~\ remains undominated. Place aqk2 in the dominating set and
then Gk1 is dominated and Gk2 has only one undominated vertex. The same
process is complete for Gk+2, Gkz until G\ and Gr are reached. At each step
only one vertex is added to the dominating set and both G\ and Gr have one
undominated vertex so adding these vertices to the set gives a domination of
Gx Pr with r + 2 vertices. In the diagram below, dark circles represent vertices
42
in the dominating set, all other vertices are dominated by the corresponding
vertex in an adjacent copy of G.
0 ajk1 ^qk 1 Cljk Ik O O K
\/ 0 O
Gi a%k1 Gkl dik Gk G
If one vertex ap in A has degree n 2 then it may be possible to rearrange
the order in which the vertices are place in the dominating set so that apl is is
used when domination C?2 If this is possible, then no additional vertex will be
necessary to dominate G\ and 7(G X Pr) = r + 1. Likewise it may be possible
to order the vertices in the dominating set so that both Gi and Gr use vertices
of degree n 2 and then 7(G x Pr) = r. If two vertices ap and aq of A have
degree n 2 then place apk and aqk+i in the dominating set for k odd. This
gives a domination number of r. Finally if any vertex of G has degree n 1
place one copy of this vertex in the dominating set for each copy, of G and this
gives 7(G x Pr) = r.
Next, for the connected graphs on n < 6 vertices that are not covered
by Theorem 3.14 the domination number is computed using the At Least
Algorithm.
43
There are 156 graphs on 6 vertices, 44 of these graphs are disconnected
and their domination number is found by taking the sum of the domination
number of there components which are graphs previously discussed. There are
112 nonisomorphic connected graphs on 6 vertices, 47 of them are covered by
Theorem 3.14. The remaining 65 connected graphs fall into 28 different groups
and the domination numbers of these graphs are given below.
Theorem 3.20 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
1
0
Theorem 3.21 Let G be one
if r =2 (mod 4)
otherwise.
of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
1
0
Theorem 3.22 Let G be one
ifr =0 (mod 2)
otherwise.
of the following graphs:
Then for all r > 1, we have
l{G x Pr) =
lOr
+
1 ifr =2 (mod 7), or r =1, 6
0 otherwise.
45
Theorem 3.23 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
lOr
+
1 if r =0, 2, 3, 4, 6 (mod 1) except r = 3
0 otherwise.
Theorem 3.24 Let G be one of the following graphs:
Then for all r >1, we have
l{G x Pr) =
Theorem 3.25 Let G be one of the following graphs:
7r
5
1 if r =0, 2, 4 (mod 5) except r = 4> 7> 9,
+ < u
0 otherwise.
Then for all r >1, we have
7(G x Pr) =
(r
T
+
1 ifr =0, 2, 4> 5, 6, 7, 9 (mod 10)
0 otherwise.
Theorem 3.26 Let G be one of the following graphs:
Then for all r >1, we have
7(G x Pr) =
llr
+
1
0
ifr =0, 2, 4> 5, 7 (mod 8) except r = 4, 7,
or r 6, 9
otherwise.
46
Theorem 3.27 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
'Ar
3~
+
1 if r =0, 2 (mod 3)
0 otherwise.
Theorem 3.28 Let G be one of the following graphs:
Then for all r > 1, we have
2 if r =0 (mod 3) except r = 3
1 ifr =1, 2 (mod 3) except r 1, f, orr =3
0 otherwise.
Theorem 3.29 Let G be one of the following graphs:
t{G x Pr) =
4r
Y
+
Then for all r > 1, we have
l{G X Pr) =
1 ifr =0, 2 (mod 3)
0 otherwise.
47
Theorem 3.30 Let G be one of the following graphs:
Then for all r > 1, we have
l{G x Pr) =
13r
To"
+
2
1
0
if r =0 (mod 10) except r = 10
if r =1, 2, 3, 4, 5, 6, T, 8, 9 (mod 10)
except r = 1, 4> 7, or r = 10
otherwise.
Theorem 3.31 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
9r
y
+
1 if r =7 (mod 14)
0 otherwise.
Theorem 3.32 Let G be one of the following graphs:
Then for all r > 1, we have
l{G x Pr) =
lAr
ll
+
1 if r 2, 3, 7
0 otherwise.
Theorem 3.33 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr)
14r
TT
+
1 if r =0, 2, 3, 5, 6, 7, 8, 9, 10 (mod 11)
except n = 5, 8, 9
0 otherwise.
48
Theorem 3.34 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
5 r
T
+
1 if r =0 (mod 4)
0 otherwise.
Theorem 3.35 Let G be one of the following graphs:
Theorem 3.36 Let G be one of the following graphs:
Then for all r >1, we have
7{G x Pr) =
1 if r =4 (mod 8)
0 otherwise.
Theorem 3.37 Let G be one of the following graphs:
Then for all r > 1, we have
l(G x Pr) =
1 if r =0, 2, 3 (mod 4) except r
0 otherwise.
Theorem 3.38 Let G be one of the following graphs:
Then for all r > 1, we have
l(G x Pr)
2 if r =0 (mod 4) except r = 4>
1 ifr =1, 2, 3 (mod 4) except r
6, 9, or r =4, 8, 12
0 otherwise.
2, 3, 6
12
1, 2, 3, 5,
49
Theorem 3.39 Let G be one of the following graphs:
Then for all r > 1, we have
l{G x Pr) =
1
0
Theorem 3.40 Let G be one
ifr =2, 3, 4, 7
otherwise.
of the following graphs:
Then 7(G X Pr) = y for all r > 1.
Theorem 3.41 Let G be one of the following graphs:
Then for all r > 1, we have
7(G x Pr) =
6 r
y
+
1 ifr4
0 otherwise.
Theorem 3.42 Let G be one of the following graphs:
Then for all r > 1, we have
7(G X Pr) =
6 r
T
+
1 if r =0, 4 (mod, 5) except r = 4
0 otherwise.
50
Theorem 3.43 Let G be one of the following graphs:
Then for all r > 1, we have
1 if r =0, 2, 3, 4 (mod 5) except r = 2, 3
0 otherwise.
Theorem 3.44 Let G be one of the following graphs:
Then for all r > 1, we have
1 if r =0, 4, 5 (mod 6) except r = 4> 5, 10
0 otherwise.
Theorem 3.45 Let G be one of the following graphs:
l{G x Pr) =
(r
T
7(G X Pr) =
or
T
Then for all r >1, we have
l{G X Pr) =
1 ifr=0,2, 3, 4, 5 (mod 6) except r = 2, 3
0 otherwise.
Theorem 3.46 Let G be one of the following graphs:
Then for all r >1, we have
7{G X Pr)
1 if r =6
0 otherwise.
51
Theorem 3.47 Let G be one of the following graphs:
Then for all r > 1, we have
l{G x Pr) =
Theorem 3.48 Let G be one of the following graphs:
8 r
y
+
1 if r =0, 2, 3, 4, 5, 6 (mod 1) except r = 2,
3, 5, 9, 10
0 otherwise.
Then for all r >1, we have
7(G x Pr)
lOr
+
1 if r =8
0 otherwise.
52
4. Domination of Circulant Graphs
In this chapter we will find the domination number of circulant graphs.
Definition 4.1 Given an integer r > 0 and subset of the positive integer S,
the circulant graph C^S1] is a graph on vertices 0, 1, r 1 with vertex i
adjacent to vertex j if i j S' or j i S (mod r).
Let Cr[l,3] be shorthand for Cr[{l,3}], where [1,3] is a graph on vertices 0,
1, ..., r 1 with vertex i adjacent to i 1 and i 3 (mod r). To illustrate
the algorithm we will find 7(C'r[l,3]) which has domination number 4 (see
Figure 4.1). To find the domination number for all r, we formulate the problem
as a dynamic programming algorithm similar to the algorithm used in Hare and
Hedetniemi [10] and Fisher [8] and Chapter 3 to find the domination number
of complete grid graphs. Then, like is done in [8], we iterate the algorithm
until a periodicity is detected. When this happens, we declare that we know
the answer for all r, and stop iterating. This allows us to find (7r[l,3] for all r
in a finite amount of time.
53
Figure 4.1 The circulant graph C7i4[l, 3]. The circled vertices are a
minimum domination.
4.1 The States and Transition Matrix
We can think of ^[1,3] as a sequences of overlapping paths on three
vertices (P3) where adjacent P3s share two vertices and there is an edge be
tween the nonshared vertices. Connecting together r of these p3s so the r
P3 overlaps the first builds CV[1,3].
We will build a domination by appending P3s clockwise from the existing
sequence. The vertices of P3 are either in the domination, dominated but not
in the domination, or not yet dominated (and thus need to be dominated by a
vertex in the clockwise direction). The states will be all valid combinations of
these three attributes.
Figure 4.2 shows the 15 valid states. This is less than the 27 ways in which
3 vertices can be assigned 3 attributes. This is because 10 potential states have
a dominated vertex next to an undominated vertex, while two others have the
54
first vertex undominated and the last vertex is dominated and these are not
valid since the last vertex cannot be dominated without dominating the first
vertex.
1=
2=
3=
4=
6= X"
7=
8=
9= KXX
13=
Figure 4.2 The fifteen states of P3. Circled vertices are in the domi
nation. Dominated vertices have a x on them.
Going clockwise, we will refer to the vertices of a state as being vertex a,
b, and c. When can two states overlap? Let vertices k, k + 1, and k + 2 be in
state j and k + 1, k + 2, and k + 3 be in state i. To be consistent, we must
have the following.
(1) Vertex b of state j must be the same as vertex a of state i.
(2) Vertex c of state j is the same as vertex b of state i unless vertex c
of state i is in the domination. In this case, if vertex c of state j is
undominated, then vertex b of state i is dominated but not in the
domination.
(3) If vertex a of state j is in the domination, then vertex c of state i is
either be in the domination, or dominated but not in the domination.
If vertex a of state j is dominated but not in the domination, then
55
vertex c of state i can either be in the domination, or not dominated.
Finally, if vertex a of state j is not dominated, then vertex c of state i
must be in the domination.
The last condition forces vertex a of state j to be dominated by vertex c in
state i if it has not already been dominated. Thus states that have vertex a
dominated match to two states, while states in which vertex a is not dominated
match to only one state.
We will represent which states match to other states by a transition matrix
A. The i,j entry of A is:
0 if state j matches to state i
and vertex c of state i is not in the domination
Aij = 1 if state j matches to state i
and vertex c of state j is in the domination
. oo if state j does not match to state i.
56
Then the transition matrix is:
1 00 oo 00 oo 1 oo CO oo oo oo CO oo oo oo
0 oo 00 00 oo 0 oo CO oo oo oo oo 00 00 00
CO 1 oo oo oo oo 1 oo oo CO CO oo oo oo 00
CO 0 oo oo oo oo CO CO oo oo oo oo CO oo CO
oo oo oo oo CO oo 0 CO oo oo oo oo OO 00 CO
oo oo 1 oo CO oo oo 1 oo CO oo oo 1 oo oo
oo oo 0 oo oo 00 oo 0 oo oo 00 oo oo oo oo
CO oo oo 1 1 CO oo oo 1 1 oo oo CO 1 oo
oo oo 00 0 CO CO oo oo oo oo oo oo oo oo 00
CO oo oo oo oo oo oo oo 0 oo oo oo oo oo oo
oo CO 00 00 0 oo oo oo 00 oo oo oo oo oo oo
oo oo oo oo oo 00 oo oo oo 0 CO oo oo oo oo
CO oo oo oo oo oo oo oo oo oo 1 1 oo oo 1
oo oo oo 00 oo oo oo oo oo oo 0 00 00 oo oo
CO oo 00 00 oo oo oo CO oo oo oo 0 CO CO oo
The information in the matrix can also be represented in the precedence
digraph Da which is shown in Figure 4.3.
Figure 4.3 Precedence Digraph for Circulant graphs. The unlabeled
arcs have weight zero. The arrowed arcs show a cycle whose cycle mean
is minimum.
57
Since the states in a domination must match, a domination of Cr[l,3] corre
sponds to a path (not necessarily simple) in Da Further a domination must
be a cycle because the first state and the r state must match to build Cr[l,3].
So a minimum domination of Cr[ 1, 3] corresponds to a cycle in Da of length r
with minimum weight.
Theorem 4.2 Let A be an 15 x 15 matrix shown above. Then for all r > 0,
we have that
7(Cr[l,3]) = trace(Ar)
where matrix exponentiation and trace are in MinPlus algebra.
Proof: By Theorem 2.7, (Ar),it is the minimum weight cycle of length r
from i to i in Da By the definition of Da, this is the number of vertices in a
domination starting and ending in state i. Thus trace(/T) = min(j4r),t, is the
minimum size of a domination of (7r[l, 3].
4.2 Periodicity of the Transition Matrix
Theorem 4.2 gives us a method for finding the domination number of
Cr[ 1,3] for any r. However, we cannot use it to find 7(C'r[l,3]) for all r. This
depends on the periodic property of Min Plus algebra. In this example, this
property can be shown by looking at A19 and A24.
58
A19 =
A24
6 6 6 6 6 6 6 6 6 7 6 7 7 7 7
5 5 5 5 5 5 5 5 5 6 5 6 6 6 6
5 5 5 5 6 5 5 5 5 6 5 6 6 6 6
5 5. 5 5 5 5 5 5 5 5 5 6 5 5 6
4 4 4 4 5 4 4 4 4 5 4 5 5 5 5
5 5 5 5 5 5 5 5 6 6 6 6 6 6 7
4 4 4 4 4 4 5 4 5 5 5 5 5 5 6
4 4 5 5 5 4 4 5 5 5 5 6 5 5 6
4 5 4 5 5 4 5 4 5 5 5 6 5 5 6
4 4 4 5 5 4 4 4 5 5 5 6 5 5 6
4 4 4 4 4 4 4 4 4 4 5 5 4 4 5
4 4 4 4 4 4 4 4 5 5 5 5 5 5 6
4 5 4 5 5 4 5 4 5 5 5 6 5 5 6
3 4 3 4 4 3 4 3 4 5 4 5 4 5 5
4 4 4 4 4 4 4 4 4 5 4 5 5 5 5
7 7 7 7 7 7 7 7 7 8 7 8 8 8 8
6 6 6 6 6 6 6 6 6 7 6 7 7 7 7
6 6 6 6 7 6 6 6 6 7 6 7 7 7 7
6 6 6 6 6 6 6 6 6 6 6 7 6 6 7
5 5 5 5 6 5 5 5 5 6 5 6 6 6 6
6 6 6 6 6 6 6 6 7 7 7 7 7 7 8
5 5 5 5 5 5 6 5 6 6 6 6 6 6 7
5 5 6 6 6 5 5 6 6 6 6 7 6 6 7
5 6 5 6 6 5 6 5 6 6 6 7 6 6 7
5 5 5 6 6 5 5 5 6 6 6 7 6 6 7
5 5 5 5 5 5 5 5 5 5 6 6 5 5 6
5 5 5 5 5 5 5 5 6 6 6 6 6 6 7
5 6 5 6 6 5 6 5 6 6 6 7 6 6 7
4 5 4 5 5 4 5 4 5 6 5 6 5 6 6
5 5 5 5 5 5 5 5 5 6 5 6 6 6 6
Note every entry in A24 is one more than the corresponding entry in A19. This
can be notated as A24 = 1 B A19. Then
A25 = Ah A24 = A b (1 b A19) = 1 b (A b A19) = 1 b A20.
Continuing in this manner, we get the following theorem.
59
Theorem 4.3 Let A be the 15 x 15 matrix given above. Then for all r > 24,
we have
Ar = 1 a Ar~5.
Proof: Since A24 = 1 B A19, the result holds for r = 24. For n > 24,
assume for induction purposes that Ar~l = 1 b Ar~6. Then
Ar = A b Ar~1 = A ei (1 0 Ar~6) = 1 a (A e Ar~6) = ls Ar~5.
Thus, we can use this repetition to find the domination number of Cr[l,3]
for all r. By observing the domination numbers when r < 24, we find an
explicit formula for Cr[l, 3].
Theorem 4.4 For all r >1, we have
7(Cr[l,3]) = {
ftl+1
if n = 4 (mod 5)
[j] otherwise.
Using this method, we let S be any set from the integers {1,2,... 9} and
we find the domination number of (7r[5] for all r. There are 512 different
circulant graphs on S and the domination numbers of these graphs is given in
Appendix A.
60
5. Knights Domination of k x r Chessboards
In this chapter we use MinPlus Algebra to develop a dynamic program
ming algorithm which solves the Knight domination problem for k x r chess
boards. The game of Chess is played on an 8 x 8 board on which a number of
pieces are placed. One piece, a Knight, can move (or attack) from corner to
diagonallyopposite corner of a rectangle three squares by two (Dr. Lasker as
quoted in [12]). See Figure 5.1.
X X
X X
X X
X X
Figure 5.1 Knight Moves: A Knight (Ã‚Â£}) attacks the squares marked
with x.
Definition 5.1 A Knight domination ofakxr board is a placement'of Knights
so each square either has a Knight on it or attacking it. The k x r Knight
domination number is the minimum number of Knights in a Knight domination
of a k X r board.
Figure 5.2 shows a domination of a 3 x 14 board with 11 Knights. Since 11 is
the minimum number, the 3 x 14 Knight domination number is 11.
61
Ã‚Â£ ft ft ft ft ft
ft ft ft ft
ft
Figure 5.2 A Minimal 3 x 14 Knight Domination.
The k x r Knight Graph, denoted Nk,r (N is the standard initial for
Knight, as K is used for King), has vertices {1,2, ...,&} x{1,2,..., r} with
an edge between vertices (g, h) and (i,j) if i g\\j h\ 2 (see Figure 5.3).
Then the k x r Knight domination number is 7( ,r).
Ns,8
Figure 5.3 The 8x8 Knight Graph. Vertices represent squares of
an 8 x 8 board with edges showing moves that a Knight can make. Its
domination number is the Knight domination number of an 8 x 8 board.
Gardner [9] generated considerable interest in Knight domination of square
boards (when k = r). Gardner and his readers found small Knight dominations
of a k x k board for k < 15. These are known to be optimum when k < 10.
k 12 3 4 5 6 7 8 9 10 11 12 13 14 15
l{Nk,k) 14 4 4 5 6 10 12 14 16 < 21 < 24 < 28 < 32 < 37
Hare and Hedetniemi [10] developed a dynamic programming algorithm for
the k x r Knight domination number which is exponential in k, but linear in r.
This allowed them to conjecture formulas for q( ,r) for k 3,4, and 6 (k 1
62
and 2 are trivial; k = 5 was not fully analyzed). They also found values
for k < 10 for some values of r. In this. chapter, we restate the algorithm in
MinPlus algebra and looks for periodic solutions to the dynamic programming
algorithm. By finding periodicity, we prove the conjectures in [10], and go on
to find the k x r Knight domination number for all k < 8 and for all r. To
illustrate the dynamic programming algorithm, we will use it to find the 2 x r
Knight domination number.
5.1 The States
A state will be a configuration of two columns of a k x r board. A square
either has a Knight on it (shown by Ã‚Â£}), is dominated by a Knight from the left
or within the two columns (shown by x), or to be dominated from the right
of the two columns (shown by a blank). Since a blank cannot be a Knight
move from a Knight, states correspond to labellings of AT)2 with 0, 1 and 2 (for
squares with a squares with a x, and blank squares, resp.) where 0 and 2
cannot be adjacent.
How many such labellings are possible? Let ar be the number of such
labellings on a path on r vertices. It is easy to show that a* = 3, a2 = 7, and
ar = 2a,i + ar_2 for all r > 2. So 17, = 41, as = 99, etc. We can then
63
solve this recursion to show that
CL*'
(l + n/2)"+1 + (l V2)
r+1
Since Nk,2 consists of two paths on [A:/2J vertices and two paths on \k/2]
vertices, the number of states needed to find the k x r Knight domination
number for all r is (Figure 5.4 shows the 81 states when k = 2):
. 2k+4
Number of States = a[fc/2ja[fc/2]
(l + v/2)
16
(5.1)
1.
10
19
28
37
46
55
64
73
s
a
X
X
X
X
2.
11
20
29
38
47.
56
65
74,
X
a a
X X
a a
X
X a
a X
X a
X X
X a
X
a
ax
a
X X
a
X
3.
12
21
30
39
48
57
66
75
a
aa
aa
a
X
X
a
X
a
a
a
a
X
4.
13
22
31
40
49
58
67
76
ax
as
ax
a
X
X
X
X
X
X
aa
X
aa
x
x
6'Sd
IX
XX
ax
X
X
X
X
X
X
X
X
X
X
X
ax
X
X
X
X
X
15
24
33
42
51
60
69
78
ax
X
a
X
X
X
X
X
X
X
X
X
X
a
X
a
a
a
X
aa
X
X
a
X
a
an oiiaa
X
a
a
8.
17
26
35
44
53
62
71
80
X
a
X X
a
X
X
a X
X
X X
X
X
a X
X X
X
9.
18
27
36
45
54
63
72
81
a
X
a
X
a
X
X
X
a
X
Figure 5.4 The 81 Knights configurations of a 2 x 2 board.
5.2 The Transition Matrix
Of the 36 possible 2x3 boards, only 81 need to be considered for domi
nation patterns. This is because each square in the third column of this board
64
is a Knight move from a square in the first column (square (1,1) and (2,3) are
a Knight move apart as are square (2,1) and (1,3)). The possible 2x3 boards
are created by matching two of the 2x2 boards shown in Figure 5.4 together
so that the middle columns overlap and the following conditions are satisfied.
(1) Any Knights in the overlapping columns must be in the same square
in each pattern.
(2) All dominated (x) squares in the last column of the second k x 2 board
must be dominated by Knights within the resulting k x 3 board.
(3) All undominated (blank) squares in the first column of the first k x 2
board must be dominated by Knights in the resulting k x 3 board.
(4) When matching a k x r board to a k x 2 board, condition 2 guar
anties that all dominated squares in the last column of the k xr board
are dominated by Knights within previous columns. Thus, an undom
inated square in the first column of the second pattern can match
with a dominated square in the last column of the first pattern since
the undominated square will become dominated in the match. (This
matching does not apply to 2 x 2 boards since the pattern is not pos
sible).
All possible domination patterns of a k X 3 board are found similarly by match
ing together two of the domination patterns of a k x 2 board following the rules
65
listed above. In Figure 5.5, possible matchings of k x 2 boards are shown to
illustrate the matching requirements.
Ã‚Â£>
X
X
X
Figure 5.5 Possible Matchings of k x 2 Boards
To find the domination number of a 2 x r board, a column is added on
to all of the possible patterns on 2 x (r 1) boards and then each of these
new patterns is checked to see if it is a domination. The minimum number of
Knights in the patterns which are dominations is then the domination number
of a 2 x r board. The r column that is being added onto the 2 x (r 1)
board is dependent on the last two columns of the 2 x (r 1) board. Thus,
the column is added by matching the last column of the 2 x (r 1) pattern to
the first column of one of the 2x2 boards following the matching constraints
listed above. To keep track of which patterns can be added onto other patterns
by using this matching, a transition matrix A is formed with rows and columns
66
being represented by the initial 2x2 boards. If the i k X 2 board matches
the j 2x2 board then the (i, j) entry in the transition matrix is the number
of Knights in the last column of the j 2x2 board. If the t 2 x 2 board
does not match the j board then the (i,j) entry in the transition matrix is
infinity or () (see below).
67
68
5.3 The States Vector
Next, a states vector xr is formed to keep track of the 81 possible minimal
states of a 2 xr board. The i entry in this vector is the number of Knights
in a 2 x r pattern whose last 2 columns are the same as the i2x2 board in
Figure 5.4. If one of these 2x2 boards does not form the last 2 columns of
any of the possible 2 x r patterns then the entry in xr is infinity or (). Now,
to find the possible patterns of a 2 x (r + 1) board, the transition matrix is
multiplied by xr using MinPlus matrix multiplication to form a new vector
xr+i. The i entry in xr+J represents the number of Knights in a 2 x (r + 1)
pattern whose last 2 columns are the same as the i 2 x 2 board. For each of
these 2 x (r +1) patterns, the first r 1 columns are dominated. If the f 2x2
board, which forms the r, and r + 1 columns of the pattern, is also dominated,
then the entire 2 x (r + 1) pattern will be dominated. Out of these dominated
states, the one with the least number of Knights in it (the lowest entry in xr+i)
forms a minimal domination of a 2 x (r + 1) board and the domination number
of a 2 x (r + 1) board is the number of Knights in this pattern.
This iterative process is started with an initial vector Xi which is used to
find the domination number of a 2 x 1 board. To do this, the first column of
the 81 2 x 2 boards is considered column 0 of the 2x1 board. Only the 2x2
boards in which the first column is dominated and contains no Knights are
69
possible patterns for a 2 x 1 board. Thus, the f entry in the initial vector Xi
is either the number of Knights in the last column of the i 2x2 board if the
first column of this board is dominated and contains no Knights, or infinity
() if the first column is not dominated or has Knights in it. Using this initial
vector, the domination number of a 2 x r board is found by iterative use of
MinPlus matrix multiplication of the states vectors and the transition matrix
A with
xr = Axr_i.
The states vectors X*, X2,..., Xi8 are shown below.
Xj = O0[2110]
X2 = 0 0 [4332322132 212110]
x3 = 2 0 [2222111111110000]
X4 = 4 0 [0000000000000000]
xs = 4'0 [211100100211100100211100100211100100]
x6 =40 [433322322322211211322211211322211211211100100211100100322211211211100100211100100]
X7 =40 [443443332332332221332332221332332221221221110221221110332332221221221110221221110]
x8 =40 [443443332443443332332332221443443332443443332332332221332332221332332221221221110]
x9 =60 [22222222222222222211 111 111 1222222222222222222111 111 111 111 111 11 111 111111100000000C]
Xio = 8 0 [000000000000000000000000000000000000000000000000000000000000000000000000000000000]
Xn = 80 [211100100211100100211100100211100100211100100211100100211100100211100100211100100]
X12 =80 [433322322322211211322211211322211211211100100211100100322211211211100100211100100]
5.4 Detecting Periodicity
To find the domination number of a 2 x r board and the pattern that
gives this domination number, repetition can be used. Looking at the 18
70
states vectors, note that x^ = X6 + 4. Since the values in xr only depend on
the values in Ã‚Â£r_i and the transition matrix, we have xr+6 = xr + 4 for all
r > 6. So the domination pattern for a 2 x (r + 6) board is found by finding
the domination pattern for a 2 x r board and adding onto this a dominated
pattern with 6 columns and 4 Knights. The domination pattern for a 2 x r
board for all r is shown in Figure 5.6.
ft a a & ft ft
ft a ft ft ft
Figure 5.6 Domination Pattern for a 2 x r board.
Periodicity is detected by storing each xr and comparing its entries with
the entries in all previous state vectors.
The algorithm for finding the domination number of a k x r board is
essentially the same as for a 2 x r. In each case, all possible fcx2 boards are
formed. The number of possible k x 2 boards is equal to 32k minus the number
of impossible patterns which is approximated by the formula 5.1. A transfer
matrix A and and initial vector Xi are then formed from these k x 2 boards.
MinPlus matrix multiplication is then used to find xr and the domination
number of a k xr board. As with 2 x r the domination number is the smallest
entry in xr which corresponds to a dominated k x r board. Successive xr
vectors are computed until a repetition is found and this is used to determine
the pattern corresponding to the domination number for all r. The results for
71
k < 7 are given in the following sections.
5.5 Domination Number for a 3 x r board
For a 3 x r board, the periodicity is given by
xr+6 = xr + 4.
This means that once an initial pattern is established, the pattern can be
extended by adding on a pattern of 6 columns which is dominated by 4 Knights.
The domination patterns are shown in Figures 5.7 5.8 5.9.
4) 4l 41 41 41 41
& 41 4l 41 41
Figure 5.7 Domination Pattern of 3 x r for r = 0, 3, 4, and 5. For
r = 0, 3, 4, and 5 (mod 6) r > 4. For r = 3 use columns 3, 4, and 5 of
this pattern. For r = 8, remove columns 4 and 5 of this pattern.
41 41 41 41 41 41
41 41 41 41 41 41 41 41
Figure 5.8 Domination Pattern of 3 x r for r 1 (mod 6) r > 4
41 41 41 41 41 41 41 41
41 41 41 41 41 41
41
Figure 5.9 Domination Pattern of 3 x r for r = 2 (mod 6) r > 14
72
5.6 Domination Number for a 4 x r board
For a 4 x r board, the periodicity is given by
xr+6 = xr + 4.
This means that once an initial pattern is established, the pattern can be
extended by adding on a pattern of 6 columns which is dominated by 4 Knights.
The algorithm detects periodicity at r = 27. The domination patterns are
shown in Figures 5.10 5.11 5.12
Ã‚Â£ & to to to to
Ã‚Â£ to to to to to
Figure 5.10 Domination Pattern of 4 x r Board, for r = 0, 3, 4, and 5
(mod 6) when r > 4. For r = 8, remove columns 4 and 5 of this pattern.
to to to to to to to
to to to to to to to
Figure 5.11 Domination Pattern of 4 x r Board, for r = 1 (mod 6)
r > 4.
to to to to to to to to to to
to to to to to to to to to to
Figure 5.12 Domination Pattern of 4 x r Board for r = 2 (mod 6)
r > 14.
73
5.7 Domination Number for a 5 x r board
For a 5 x r board, the periodicity is given by
xr+i8 = xr + 14.
This means that once an initial pattern is established, the pattern can be
extended by adding on a pattern of 18 columns which is dominated by 14
Knights. The algorithm detects periodicity at r = 39. The domination patterns
can be described in terms of the separate components shown in Figure 5.13.
to
Ã‚Â£ to Ã‚Â£ Ã‚Â£ Ã‚Â£
to to
to to
to to
E =
L =
B =
to to
to to
Figure 5.13 Components in Domination Patterns of 5 x r Boards.
These components are combined together to form the patterns for a 5 x r
board. The table below shows the patterns with the repeating 18 columns of
74
each pattern in parentheses.
n 7(^5,r) Pattern for r = 0 Pattern for r > 0
18r+ 0 14r + 1 Meaningless (.PBE4B)r'1PBE5
18r+ 1 14r + 2 See 1x5 result (.PBE4B)r1PBE6
18r + 2 14r + 3 See 2x5 result (.PBE4B)r1PBE7
18r + 3 14r + 4 See 3 x 5 result lPBE4B)r'1PBE8
18r+ 4 14r + 4 (.PBE4B)rU
18r + 5 14r + 5 E5 (.PBE4B)r~1PBE5BU
18r + 6 14r + 6 E6 (.PBE4B)r~1PBE6BU
18r + 7 14r + 7 E7 (PBE4B)r~1PBE7BU
18r + 8 14r + 8 E8 (PBE4B)r~1PBE8BU
18r + 9 14r + 8 (PBE4B)rUBE4
18r + 10 14r + 9 (.PBE4B)rUBE5
18r + 11 14r + 10 {PBEAB)rUBE6
18r + 12 14r + 10 (PBEAB)rP
18r + 13 14r + 11 UBE8 (PBE4B)r~1PBE4EBP
18r + 14 14r + 12 (PBE4B)rUBE4BU
18r + 15 14r + 13 (PBE4B)rUBE5BU
18r + 16 14r + 14 (PBE4B)rUBE8BU
18r + 17 14r + 14 (.PBE4B)rPBE4
5.8 Domination Number for a 6 x r board
For a 6 x r board, the periodicity is given by
xr+4 = xr + 4.
This means that once an initial pattern is established, the pattern can be
extended by adding on a pattern of 4 columns which is dominated by 4 Knights.
The algorithm detects periodicity at r = 15. The domination patterns are
shown in Figures 5.14 5.15 5.16
75
to to to to to
Ã‚Â£ to to to to to
Figure 5.14 Domination Pattern of 6 x r Board, for r = 0, and 3 (mod
4) r > 4.
to to to to to to to
to to to to to to to
Figure 5.15 Domination Pattern of 6 x r Board, for r 1 (mod 4)
r > 4
to to to to to to to to
to to to to to to to to
Figure 5.16 Domination Pattern of 6 x r Board, for r = 2 (mod 4)
r > 4
5.9 Domination Number for a 7 x r board
For a 7 x r board, the periodicity is given by
xr+5 = xr + 6.
This means that once an initial pattern is established, the pattern can be
extended by adding on a pattern of 4 columns which is dominated by 4 Knights.
76
The algorithm detects periodicity at r = 32. The domination patterns can be
described in terms of the components shown in Figure 5.17.
to
ft
&
& to
to to
to to
to to
to
to to
to
to
to to
to to
to to
to
to to
to to
to to
to
to to to
to v = to w = Y = Z =
to to to
Figure 5.17 Components in Domination Patterns of 7 x r Boards.
These components are combined together to form the patterns for a 7 x r board.
The table below shows the patterns with the repeating 6 indicated by raising
the pattern to the exponent r.
77
n 7(^7 ,r) Pattern for r > 0
16 See Previous results
7 10 ZVZWZVZ
8 11 ZVZWZUZZ
9 12 ZZUZWZUZZ
10 14 ZVZWYYWZVZ
11 15 ZVZWYYWZUZZ
12 16 ZZUZWYYWZUZZ
5 r + 13 6r + 18 T{pys
5r + 14 6r + 19 Rxipys
5 r + 15 6r + 20 Ri(P)rQ
5r + 16 6r + 21 T(P)rR2
5r + 17 6r+ 22 Ri(P)rR2
5.10 Conclusions
The algorithm developed in this paper can be used to find the domination
number of a k x r board for all r. The algorithm is not linear in k so the size
of k is limited by computer time.
Theorem 5.2 The domination number of akxr board 7(iVfc,r) for k < 6 and
for all r is:
7(AW
n
4 [
^ 6
2r+4
3
2r+4
3
2r+5
3
r
6
2r+4
r
if k
if k
if k
ifk
if k
ifk
ifk
ifk
ifk
1
2
2,
3,
3;
3,
3;
4,
5,
r = 1 (mod 6) and r > 6
r = 1 (mod 6), and r > 6
r 2 (mod 6), and r > 6
r ^ 1, 2 (mod 6) and r > 6
r = 1 (mod 6), and r >7
r/ 1 (mod 6) and r >7
r or 17 (mod 18), and r > 13
7r 9  1 ifk
+ 1 ifk
r 4 ifk
6r+12 5 ifk
5, r 4,12,13, or 17 (mod 18), and r > 13
6, r = 1 (mod 4), and r >4
6, r 7^ 1 (mod 4), and r >4
7, and r > 13
78
6. Time Complexity of The Power Method
All of the applications presented in the previous chapters rely on the pe
riodic property of dynamic programming algorithms in MinPlus algebra. The
algorithms terminate when periodicity occurs, and this allows us to find in
finitely many solutions in a finite number of iterations of the algorithm.
Definition 6.1 The preperiodic interval of a matrix A denoted Ro(A), is the
smallest R0 such that Ar q s Ar~p for all r > Ro.
In Theorem 2.3 we stated that if A is the irreducible transfer matrix for the
dynamic programming algorithm, then periodicity implies that for some Ro, p
and q, we have that Ar = q 12 Ar~p for allr > Ro.
Although some matrices that are not irreducible are periodic, this is not
the case in general. Consider a matrix whose precedence digraph has two
disjoint cycles of length one with different weights as shown below.
t CO A100 = ' loot CO
oo s oo 100s
In this case A[x = rt and A22 = rs and thus if t ^ s there is no r where
Ar = q la Ar~p and Ro(A) is undefined. Thus when using the periodicity of
MinPlus algebra to solve problems, we must make sure that the matrix is
irreducible.
79
6.1 Finite Termination of the Power Method
In Chapter 2, we stated Theorem 2.3 without proof. We have referred to
this theorem throughout all of the applications presented in this thesis because
it guarantees that at some point the dynamic programming algorithms will
be periodic and the power method can terminate. We will now restate the
theorem and give a proof based on minimum cycle means and the diameter of
A.
Definition 6.2 The diameter of A denoted diam(A) is the maximum weight of
all minimum weight paths in A.
Theorem 6.3 If A is an irreducible matrix on m vertices, then there exists a
p, q and Rq so that
Ar q & Ar~p for all r > Ro
Proof: Let p(A) be the minimum cycle mean of A. Since A is irreducible,
we know that the precedence digraph of ^4m2m+x has paths between all pairs
of vertices. Let k be the smallest common divisor of the lengths of the critical
cycles of A, with the property that k > m2 m + 1. Now, consider Ak and let
S be the subgraph of the precedence digraph of Ak, (DAk) whose vertices are
contained in at least one of the critical cycles in A, and let T be the subgraph
of DAk on the remaining vertices. In Da, every vertex of S has a minimum
weight path of length k to itself using the critical cycle that contains it. The
80
weight of this path is kp(A), so the minimum cycle mean of Ak is kp(A) and
the critical cycles are 1loops on every vertex in S.
We will now look at possible paths in the precedence digraph of Ak. Let
P be any minimum weight path of length r between V{ and Vj that uses no
vertices of S. This path consists of cycles and simple paths in T. Let b(Ak) be
the minimum entry in Ak, then the weight of each arc in the simple paths is at
least b(Ak). If there are h such arcs, the weight of the simple path portion of P
is at least hb(Ak). Next if p{T) is the minimum cycle mean of T, then the arcs
in every cycle in T have weight of at least p{T). P contains r h arcs that are
part of some cycle, so the weight of these cycles is at least (r h)p(T). Thus
w{P) > (r h)p(T) + hb(Ak)
This is a nonnegative value since p(T) > b(Ak), and since h < m 1 we have
to(P) > rp(T) + (m l)(i>(<4) p{T))
Next consider all minimum weight paths Q from u, to vj of length r that use
vertices of S. These paths consist of a path from u, to S, loops on the vertices
of 5, and a path from S to Vj. The paths to and from S have weight at
most diam(v4A:) and the loops in S each have weight kp(A) = p(Ak). Since Q
contains at most r loops in S, we have
w(Q) < 2diam(Afc) + rp(Ak).
81
Now we want to show that all minimal paths of length r > R0 must use vertices
of S. This is true if w(P) > w(Q) for all paths of type P and Q of length r.
Since p(T) > p{Ak) we have w(P) > w(Q) for all r > R0 when
2diam(J4fc) (m l)(t(A) p(T))
p(T) p(A)
Now for any minimal path of length r > i?0 a minimal path of length r +1
is obtained by looping on some vertex of S that was in the minimal path of
length r, and the weight of these paths differs by p{Ak) in every case. Thus
(Ak)r = p(Ak) (Afc)r+1 for allr > R0
This proves that Ak is periodic and hence A is periodic with Ro(A) = kRo(Ak).
The construction of the proof also gives a bound for Ro defined in terms of the
entries of Ak.
We will now illustrate this method by finding the bound on i?0(A) where
A is the 15 x 15 transition matrix used in Chapter 4 to find the domination
number of (7,. [1,3]. The precedence digraph in Figure 4.3 shows that there is
only one critical cycle C for this matrix and that p(A) = . To find the
bound for Rq using the method described above, we need to find the smallest
k > m? m + 1 with k being a common divisor of the lengths of the critical
cycles in A. Thus since m = 15, and there is only one critical cycle of length
5, we have k = 215. This bound is large because we are required to have
6.1)
82
k > m2 m + 1 to guarantee that the matrix has paths of any length greater
than k. (In this example, we have that A8 has no infinity enties so k could be
reduced to 10). We will illustrate the bound using k 215 and considering
possible paths in A215 which is shown below.
A215 = 39 a A20 39 B
6 6 6 6 6 6
5 5 5 5 5 5
5 5 5 5 5 5
5 5 5 5 5 5
4 4 4 4 4 4
5 5 5 6 6 5
4 4 5 5 5 4
4 5 4 5 5 4
5 5 5 5 5 5
4 5 4 5 5 4
4 4 4 4 5 4
4 4 4 5 5 4
5 5 5 5 5 5
4 4 4 4 4 4
4 4 4 4 4 4
6 6 7 7
5 5 6 6
6 5 6 6
5 5 5 6
5 4 5 5
5 5 6 6
4 5 5 5
5 4 5 6
5 5 5 5
5 4 5 5
4 4 4 5
4 4 5 5
5 5 5 5
4 4 4 4
4 4 5 5
7 7 7 7 8
6 6 6 6 7
6 6 6 6 7
5 6 6 6 6
5 5 5 5 6
6 7 6 6 7
5 6 5 5 6
5 6 5 6 6
5 6 5 5 6
5 6 5 5 6
4 5 5 5 5
5 6 5 5 6
5 6 5 5 6
5 5 4 4 5
5 5 5 5 6
Using this matrix and the fact that S = {us, V7, vs, vu, U14}, we have that
07 1
b(Ak) = 43, p(T) = ^ and p(Ak) = kp(A) = (215) = 43.
A 0
Next, for any matrix, its diameter is less than the number of vertices times the
maximum entry. Using this bound we have
diam(AA:) < 690.
This is considerably larger than the actual diameter of A which is 269 and thus the
bound on Rq will be significantly effected by using the bound instead of the actual
83
diameter. Putting all of this together in equation 6.1 we have
Rq =
2 x 705 (14) x (43 43.5)
43.5 43
= 2806
Thus, the bound is Ro(Ak) = 2806 and the bound for A is i?o(A) = (215)(2806) =
603290, whereas the actual value of iEo(4) was 24. In computing this bound, we
used k = 215 as the smallest value where Ak had no infinity entries instead of the
actual value of k = 20. If we compute the bound with k 10 and diam(A) = 269
it gives Rq(A) (10)(1062) = 10620 which is still quite large. Thus although the
bound shows that Rq is finite, it is not a very good approximation of Rq(A) in most
cases.
The large value of k and the large approximation of diam(A) both contribute
to the difference between the bound and the actual value of i2o(^4) In addition, if
the growth of the entries in successive powers of A is large, then the fact that the
bound requires us to consider paths in Ak where k is large may have a much more
significant effect. In the next section we will restrict A and find a bound for Rq that
is strong.
6.2 Exact Bound for Rq when A is MCMHamiltonian
In this section, we will determine the bound for Rq when A is restricted to a
special case.
Definition 6.4 A minimum cycle mean (MCM) Hamiltonian digraph is a digraph
where the minimum cycle mean occurs on a Hamiltonian cycle.
Definition 6.5 A matrix is MCM Hamiltonian if its precedence digraph is.
84
Note that if A is MCM Hamiltonian then the subgraph T as described in the
previous section does not exist. Thus, we cannot use the bound from that section
and we must look at other' ways to bound Rq in this case.
Definition 6.6 If B is an MCM Hamiltonian matrix, define the following.
(1) Let A = 1 be the unique eigenvalue of B.
(2) Let P be the permutation matrix that permutes the rows and columns of B so
that F = PT (A)) El P is a matrix with minimum cycle mean 0 and
the minimum cycle C oriented on the vertices 1, 2, ..., m where vertex 1 is
represented by row 1 of F. Thus the arcs of C are (i, i + 1) for all i < m, and
the arc (m, 1).
(3) Let D be a diagonal matrix with Djj = (Fi ,i+1) + Fm,l
Definition 6.7 If B is an MCM Hamiltonian matrix on m vertices, then the nor
malized form of B is A = (D) El PT El (5 S (A)) El P 13 D. A normalized MCM
Hamiltonian matrix is a matrix that is equal to its normalized form.
Theorem 6.8 Let B be an MCM Hamiltonian matrix with minimum cycle mean X,
and normal form A = (Z?) El PT El (B El (A)) El P El D. Then
(1) j4i,2 = ^2,3 = = Am1,771 = <4.771,1 = 0
(2) The minimum cycle mean of A is 0.
(3) A is nonnegative.
(4) Br = AAr G3 (D) 0 PT El Ar E! P El D and
(5) Br = q Br~P then Ar = Ar~r.
85
Proof: Let Cb, Cr, Ca be the cycles with minimum cycle means in F, R, and
A. Then F = PT El B 12 (A) El P so the arcs of Cb have been permuted to the arcs
of Cf + {F\,2i ^2,3; Fmi 1} and there weights have been reduced by A. Also for any
arc in A,
Ai,j = Fij Dlti + Djj
Fi,j [^i,t+i + . + Fmi>m + Fm> 1] + [Fjtj+i + ... + Fmii?n + Fmt{\.
which gives a direct correspondence between the cycles in A and the cycles in B.
1) By construction,
Ai,i+\ = Fij+x [Fi,t+1+. +F1)7il,m + ^m,l] + [F1i+l1i+2 + +Pml,m+^77i,l] = 0
for all i < m, and
Am,l = Fmii [^771,1] + [Fl,2 + ^2,3 + + Fmltm + Fmii] 2 0.
Thus all the arcs in C = {Ai^ ALmi.mj Am,i} have weight 0 and form a
hamiltonian cycle of weight 0.
2) For any cycle C>i; of length i in A which has weight w(CAi), there is a corre
sponding cycle Cb{ of length % in B which has weight w(Cb;) = w(Ca;) + Ai.
By part 1 C is a hamiltonian cycle with cycle mean 0 so w(Ca) < 0. If
w(Ca) < 0 then vj(Cb) + w(Ca) + Am has cycle mean less than A which is a
contradiction. Thus the minimum cycle mean of A is 0.
3) Suppose A has some entry A, j < 0 Consider the cycle
{Aij, Ajj+i,... Amii,Ai Since all the arcs of this cycle ex
cept A{j are contained in the cycle with minimum cycle mean, they have weight
86
0 and the weight of the cycle is A{j < 0 which is a contradiction since the min
imum cycle mean of A is 0.
4 & 5) The results are found by using the commutative and associative properties
in MinPlus algebra.
To illustrate the normalization process consider the following matrix which has the
arcs in the cycle with minimum cycle mean underlined.
B
7 9 2 8
18 6 6
9 8 6 5
3 4 7 8
which has A = 3
Then using the definitions above,
P =
oo 0 00 oo ' 5 ^2 3 3 0 oo oo oo
0 oo 00 oo 6 4 1 5 D = oo 2 oo 00
F =
oo 00 0 oo 5 6 3 2 00 oo 3 oo
0 00 oo 0 _ I 0 4 5 . 00 oo 00 1
and finally.
A =
5 0 6 4
4 4 ol 4
2 5 3 Ol
0 1 6 5
Thus any minplus problem with an MCM Hamiltonian transfer matrix B can
be solved using the normalized MCM Hamiltonian matrix A and then converting
the solution using parts 4 and 5 of Theorem 6.8
To find a bound on the number of iterations needed for the transfer matrix to
repeat, we will consider normalized MCM Hamiltonian matrices. The precedence
87
digraph for a normalized MCM Hamiltonian matrix has a Hamiltonian cycle with
all arcs having weight zero and other arcs of nonnegative weight. The bounds on
repetition will be formed by considering paths in the precedence digraph since if
Ar = Ar~m where m is the size of the matrix, then in the precedence graph the
minimum weight of a path of length r from i to j is the same as the minimum weight
of a path from i to j of length r m for all i and j.
Definition 6.9 Let A be a normalized MCM Hamiltonian matrix and let hi be the
i arc in the minimum cycle of A which by definition has weight 0.
Definition 6.10 Let ak be the number of arcs from i to j used in the path where
i j + 1 = k (mod m).
Theorem 6.11 A path P of length r from vertex i to vertex j in a MCM Hamilto
nian digraph on m vertices satisfies the following:
ci + 2a2 +b (to l)am_i r = (i j) mod m (6.2)
Proof: First we will prove this for r = 1. There are two cases for a path of length
1. Case 1: j = i + 1. In this case, the arc from i to j is one of the b{ arcs from the
minimal cycle. Thus no ai arcs are used and equation 6.2 becomes 1 = (i (i + 1))
(mod m) which holds. Case 2: j ^ i + 1 In this case, if j = i + k the arc used is a;
where l = i j + 1 = k + 1 and equation 6.2 becomes l 1 = (i j) (mod to) or
(i (i + k) + 1) 1 = (i (i + k)) (mod to) which also holds.
88
For induction purposes, we will assume that a path of length r 1 satisfies
( 6.2) Now any path of length r contains a path of length r 1 plus one arc. Let
j be the vertex at the end of the path of length r 1 and k be last vertex in the
path. Thus the last arc in this path of length r goes from vertex j to vertex k. Let
a[, a'2,..cC_i be the coefficients that satisfy ( 6.2) for the path of length r 1.
Then ( 6.2) is satisfied for the path of length r if a'x + 2a'2 Hb (to T)a'm + arc
from j to k r = (i k) (mod m). Using the induction hypothesis, this will hold if
i j (mod to) + arc from j to k 1 = i k (mod to) .
Case 1: The additional arc is one of the b{ arcs on the minimal cycle. This
gives k = j + 1 (mod to) by the labeling of the minimal cycle, and ( 6.2) becomes
i j 1 = i (j + 1) (mod to) since there is no additional a; arc added to the
equation.
Case 2: The additional arc is a\ where l = j k + 1. Then ( 6.2) becomes
i j + (j k + 1) 1 = i k (mod to) which holds.
Definition 6.12 Given a path P the arcs can be divided into three exclusive classes:
Type 1 Arcs used in a minimal path from the initial vertex to the final vertex in the
path. There are < m of these.
Type 2 Arc used in any minimal cycle containing an a,. There are Y^=i^ai f
these.
Type S Arcs that are in the minimal cycle.
Let N{ be the number of arcs of type i used in a given path.
Each type is exclusive so once an arc has been designated as being in type one it
cannot be in type two or three. Also the determination of the arc type must be done
in the order they are described. For example using the graph in Figure 6.1, consider
the following path of length 13 from vertex 1 to vertex 3.
P = b1b2a3bI 62b3b4b5a^ b4 b5b1b2
In this path, the bold arcs are designated type one, the italic arcs are designated
type two and the remaining arcs are type three.
1
02
Figure 6.1 Precedence digraph of a normalized MCM Hamiltonian
matrix. The b{ arcs represent the arcs in the minimal cycle and have
weight 0. The weight of the a, arcs is nonnegative.
Now we need to show that for some r, any minimum weight path in the prece
dence digraph of length r has the same weight as a minimum weight path of length
r m. This will be possible if there are always at least m type 3 arcs in a minimum
weight path of length r. Since m of these arcs form a Hamiltonian cycle and we
assumed the type 3 arcs have weight zero by theorem 3.2, their removal will not
change the weight of the path and will give a path of equal weight and length r m.
90
To find the value of r for which we are guaranteed to have N3 > m, we must
first show that we need only consider paths which use less than m of the a* arcs. To
do this we will show that for any path of length r from i to j that uses more than m
of the di arcs there is a path of length r from i to j that uses less than m of these
arcs. This is done by showing that such a path satisfies 6.2 in the same way that
the original path does.
Definition 6.13 Let a = (a1?..., am_i) and c = (ci,..., cm_i) and let /(a) =
a\ + 2ci2 + + (m 1 )ami. Next define a < c if ai < C{ for all 1 < i < m.
Now a and c are the coefficients of two paths described by 6.2. We want to
show that if a satisfies 6.2 then there is some c that also satisfies 6.2 and uses less
than m of the a, arcs.
Theorem 6.14 If a satisfies a* = m> then there is a c with 0 < c < a and
/(c) = /(a) (mod m).
Proof: Construct a chain of Cj, j = 1,2,..., m with a = co > ci > C2 > ... >
cm by recursively defining cj+i by decreasing one nonzero entry of cj by 1. After
ai + + + 1 = m we reach 0 so cm 0
Now consider the values of /(ci),...,/(cm) (mod m). If they are all distinct,
then /(a) is congruent mod m to a unique /(cj), satisfying the conclusion of the
theorem. Otherwise, for some j,k with 1 < / < k < m, we have /(cj) = /(ct)
(mod m), and cj > c^. Then a > cj Ck > 0 and a > a (cj Ck) with
91

PAGE 1
I I CELLULAR AUTOMATON MODELS OF THE HEART by Anne Spalding A thesis submitted to the F acuity of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Applied Mathematics 1996
PAGE 2
This thesis for the Master of Science degree by Anne Spalding has been approved by William Briggs Stephen Billups Weldon Lodwick Date
PAGE 3
Spalding, Anne (M.S., Applied Mathematics) Cellular Automaton Models of the Heart Thesis directed by Professor William Briggs ABSTRACT Understanding the mechanisms of the heart is essential to the diagnosis and cure of many cardiac disorders. Cellular automaton models represent the cells of the heart with an integer grid and model the ionic interactions of these cells. Cellular automaton models provide a tool for understanding the correlation between an ECG (electrocardiogram) and the actual movement of charge through the heart. These models are used to gain a better understanding of the mechanisms of activation of the cells of the heart. They can also be used as a diagnostic tool in modeling and predicting cardiac disorders. This thesis develops a cellular automaton model of the ventricles of the heart. Discrete rules are developed to model cellular interaction and recreate the activation sequence of the heart. The model is then used to produce ECGs corresponding to several known cardiac disorders in order to show the effectiveness of the model in representing actual cardiac phenomena. This abstract accurately represents the content of the candidate's thesis. I recommend its publication. Signed William Briggs
PAGE 4
ACKNOWLEDGMENTS I would like to take this opportunity to thank Professor Briggs and Professor Hermes for their guidance in completing this research. I would also like to thank my friends and my family for their constant support and encouragement in all of my pursuits.
PAGE 5
CONTENTS 1. Introduction ... 1 2. Physiolot,>y of the Heart ............. 3 2.1 Anatomy .. 3 2.2 Electrophysiology 4 2. 3 Modeling Ionic Interchanges 8 3. The Cellular Automaton Model 13 3.1 Model Setup 13 3.2 Depolarization Algorithm 16 3.3 Repolarization Algorithm 18 3.4 Output 21 4. Experiments on the Twave 22 4.1 Model A 23 4.2Model B 24 4.3 Model C 25 4.4 ModelD 26 5. Experiments 29 5.1 SelfSustained Ventricular Tachycardia (SSVT) 29 5.2 Ventricular Tachycardia 32 5.3 WolffParkinsonWhite Syndrome (WPW) 34 5.4 Right and Left Bundle Branch Block 35 6 Conclusions 37 7 References 39 v
PAGE 6
FIGURES Figure 2.1 Anatomy and conduction system of the heart 4 2.2 Ionic currents and channels in an action potential [3] 5 2.3 Standard 12Lead ECG placement[?] 6 2.4 Formation of the ECG during the activation sequence of the heart 7 2.5 The normal IlLead ECG for a SA pulse .. 8 2.6 Representation of the septum and the left ventricle of the heart used in an ionic model l 0 3.1 Cross section of the CA model and the model setup for one cell l 5 3.2 Elliptical wavefronts on a cross section of a modified MillerGeselowitz model (9] 3.3 Schematic representation of repolarization of a target cell in theCA model 4.1 ECG for a normal SA pulse without Pwaves ...... ( 4.2 ECG for a normal SA pulse from Model A: inverted Twave Vl 18 20 23 24
PAGE 7
4.3 ECG for a normal SA pulse from Model C: inverted Twave. 26 4.4 Repolarization intervals on an xz cross section of Model D 27 4.5 ECG for a normal SA pulse from Model D: T wave is correct ...... 28 5.1 SSVT for Model A with one ectopic and one SA pulse 5.2 Nonrhythmic SSVT for Model D .. 30 31 5.3 ECG ofSSVT for Model D with a capture beat occurring at 1860 ms. 32 5.4 Ventricular tachycardia initiated in the upper left ventricle 5. 5 Ventricular tachycardia initiated in the lower right ventricle 5.6 Lateral slice of the ventricles showing areas for accessory pathways inWPW .. 5.7 Lead I, II, and Ill ECGs for WPW in a modified MillerGeselowitz model [9] 5. 8 ECG for simulated right bundle branch block and a normal SA pulse in Model D Vll 33 33 34 35
PAGE 8
1. Introduction Cardiac disorders such as infarction, tachycardia, and fibrillation are the cause of thousands of critical care situations every year. One of the main methods of diagnosing and analyzing these disorders is through an electrocardiograph (ECG) scan which records the electronic changes in the heart tissue. The effectiveness of the ECG in diagnosis is hindered by the fact that there is still an incomplete understanding of the relationship between the output of the ECG and the changes that actually occur on a cellular level. This thesis presents a mathematical model called a cellular automaton (CA) that helps to explore this relationship. This cellular automaton model, represents the ventricles of the heart by 2605 points on a threedimensional integer grid. These points are given attributes and interaction rules so that they resemble the actual cells of the heart. The interaction of the grid points in the cellular automaton model is represented by a simulated lead II ECG scan which can be compared to ECG scans from the actual heart. The main advantage of this model is the ease of simulation and the ability to relate ECG output with the point of origin and path of electrical pulses. The cellular automaton is effective in modeling several disorders of the heart and reproducing the ECG scans for irregular heart activations. Analysis of the ECG output of theCA model in comparison with that of patients may prove to be a valuable tool in diagnosis and treatment of these disorders. The second chapter of this paper discusses the underlying anatomy necessary for understanding the functions of the heart in order to develop an accurate cellular automaton model. This includes a brief description of the action pQtential and the propagation of action potential through the conduction system of the heart. In the next chapter one particular CA model is developed in detail. This model simulates both the depolarization and repolarization sequences in the ventricles. In the fourth
PAGE 9
chapter, some of the limitations of various CA models such as the inverted Twave are discussed. Several possible solutions to the discrepancies between the models and the actual heart are addressed and evaluated. In the last chapter experiments are presented that show the use of the CA models in simulating several cardiac conditions. These experiments include ventricular tachycardia, capture beats, WolffParkinsonWhite syndrome and right bundle branch block. The CA model developed in this thesis was originaliy created by Professor Henry Hermes of the University of Colorado at Boulder [13]. The research contribution of this thesis is to modify the code and algorithms of the original model and use this modified model in simulations of cardiac disorders. This thesis also uses at new methods in the repolarization algorithm to analyze the inverted Twave which is a limitation in the original model.
PAGE 10
2. Physiology of the Heart 2.1 Anatomy The human heart has four chambers. The right and left atria are the receiving chambers; the right atrium receives blood from the body and the left atrium receives blood from the lungs. From the atria the blood flows into the right and left ventricles which are separated by a thick muscle wall called the septum. The right ventricle pumps blood out to the lungs while the left ventricle pumps blood to the entire circulatory system of the body. Thus the left ventricle walls, which contract to pump the blood, are three times as thick as in the right ventricle. The wall of each chamber consists of three layers. The inner layer is the endocardium, the middle layer is the myocardium and is composed of cardiac muscle, and the outer layer is the epicardium [3]. The electrical conduction system of the heart begins at the sinoatrial (SA) node which is located above the right atrium; it then follows internodal pathways through the muscles of the atria to the atrioventricular (AV) node which is located on the right side of the interatrial septum. In the ventricles, the conduction system consists of the bundle of His, which begins at the AV node, and travels down the septum until it branches into the right and left bundle branches [7]. These branches pass through the right and left ventricles until they branch out into the Purkinje fibers which spread along the muscle walls at the bottom of the ventricles (Figure 2.1). 3
PAGE 11
BUNDLE Of HIS R. BUNDLE BRANOI Figure 2.1 Anatomy and Conduction system of the Heart f3]. 2.2 Electrophysiology The pulsing of the heart is controlled by an electrochemical wave which passes from cell to cell within the heart. In a normal beating heart, the electrical stimulus originates in the SA node at the top of the atria and propagates through conductive pathways to the Purkinje fibers. The electrochemical wave then spreads through the myocardial muscle of the ventricles. This wave causes the sequential and periodic contraction of the muscles in the heart which forces blood from the atria to the ventricles and out into the circulatory system of the body. This electrochemical wave is formed by the propagation of an action potential. An action potential is the changing of polarity or charge of a cell due to ionic exchanges between the extracellular fluid and the cell. The action potential of a cell is defined by five phases. In phase 0 sodium enters the cell through fast sodium channels which causes the polarity of the cell to change from its resting potential of 90 mY to a potential of +10 mY. This phase is referred to as the depolarization of the cell. In phase 1 potassium leaves the cell and leads to a slight decrease in the potential to about 0 m V. This stage is referred to as early In phase 2, which is the plateau, there is a balance between potassium leaving the cell and calcium entering it. In phase 3, rapid repolarization, the cell continues to repolarize as potassium leaves the cell and the potential decreases back towards 4
PAGE 12
the resting potentiaL Finally in phase 4, return to resting potential, the amount of potassium leaving the cell is only slightly higher than the amount entering it and the potential levels off at 90 mV During this phase the ionic balance of the cell is actively restored by the sodium potassium pumps in the cell [3] (Figure 2.2). From the end of phase 1 and through most of phase 3 the potential of the cell is too large to allow the cell to depolarize again; this is called the absolute refractory period. During the end of phase 3 and the beginning of phase 4, the cell can depolarize but the height of the action potential will be less than normal; this is called the relative refractory period. Figure 2.2 Ionic curents and channels in an action potential [3]. The chemical and electro static movement o(potassium (K+), sodium(Na+), and calcium(Ca+ +) are shown .fhr each phase (J( the action potential. Any cell can begin and maintain a rhythmic firing of action potentials which is referred to as automaticity. Certain cells of the heart, such as the SA node, have an inherent automaticity while other cells develop automaticity causing abnormal conditions in the heart. In a normal heart, the SA node produces an acti?n potential at regular intervals. The depolarization of the SA node changes the concentration of sodium and potassium in the extracellular fluid around the SA node causing the sequential depolarization of neighboring cells. Action potentials propagate quickly along the conduction pathways and more slowly through the cardiac muscle cells,
PAGE 13
so that the electrochemical wave will depolarize the entire heart in about l 00 ms. In addition, propagation through the cardiac muscle is anisotropic so that the action potentials travel faster along the muscle fiber than across muscle fibers. If a cell other than the SA node begins firing or producing action potentials, then this cell is called an ectopic pulse or focus. A depolarization wave will spread through the heart from an ectopic focus only by cell to cell propagation instead of traveling along the normal conduction system as an SA pulse does. Depending on the position of the ectopic pulse, the time to depolarize the entire heart can be longer or shorter than the 80100 ms required when the depolarization wave originates at the SA node. One of the main methods of observing the electrophysiology of the heart is through an electrocardiogram (ECG) which is an electronic record of the activity of the heart. This activity is recorded by placing up to twelve electrodes on the surface of the body and recording the difference in electrical potential between various combinations of these electrodes. Figure 2.3 The 12/.J!ad ECG placement and the standard placement of the three main leads. The bottom enlargement shows the placement of the six precordial leads [7/. The specific placement of electrodes on the body (Figure 2.3) gives a standard twelve lead ECG. The three standard leads (I,II,III) view the heart in the frontal plane and are generally used for monitoring purposes. The three augmented leads 6
PAGE 14
I (a VR,a VL,a VF) use different combinations of the three standard electrodes in order to record small deflections in the potentials. Finally, the six precordial leads (vl...v6) view the heart in the horizontal plane. The pulsing of the SA node is not detectable in the ECG, since the atria depolarize almost instantaneously after the firing of the SA node. The average direction of the depolarization wave in the atria is down towards the patients left side and is recorded in the Pwave of the ECG [7]. The repolarization of the atria is usually not visible in the ECG because it takes place at the same time as the depolarization of the ventricles. The ventricular depolarization is represented by the QRS complex in the ECG. In different leads this complex can contain any combination of the Qwave, Rwave, and Swave. In all the Qwave is the negative wave preceding the positive Rwave, and the Swave is the negative wave following the Rwave. Finally the Twave represents the repolarization of the ventricles (Figure 2.4). The PR interval represents the time for atrial depolarization and the pause at the AV node before ventricular depolarization. The QRS interval represents the time for ventricular depolarization. The QT interval represents the time between depolarization and the beginning of repolarization in the ventricles. This interval varies inversely with heart rate and its length must be compared with a normal heart to detect abnormalities. j> (mVJ EGG Figure 2.4 Formation of the ECG during the activation sequence of the heart. lhe shaded areas represent depolarized cells. The 11!J:ad ECG jiJr each step in the depolarization process is shown below the figure that shows which cells are depolarized [8}. 7
PAGE 15
Each of these waves and intervals has a different normal state for each lead and the comparison of all twelve leads is used to gain a complete picture of the electrical activity and abnormalities in the heart The II lead for a normal SA pulse has positive Rand Twaves with the Twave having about 1/3 the height and 3 times the length of the Rwave (Figure 2.5). R \ T I Q Figure 2.5 The normal IIDcad ECG for a SA pulse. The ECG has a positive Rwave and a positive 2.3 Modeling Ionic Interchanges An ECG is one method used to understand the ionic interchanges between cells based on the changes in polarity of portions of the heart, but this process can also be understood through models that replicate cellular interactions. The first work towards modeling the electrochemical interchange in cells was done by Hodgkin and Huxley in the 1960's [14]. By performing voltage clamp experiments on a squid giant axon, they were able to develop a set of differential equations that model the action potential within the cell and activation conduction between cells. These equations for a neuron were then adapted in various ways to fit the action potentials of Purkinj e fibers and cardiac cells. In a study by McCallister, Nobel and Tsien, the adaptions were formed by using similar voltage clamp techniques on the heart. The Purkinje fjbers were found to have nine different ionic currents [17][18] in contrast to the three currents described in (14] for the neuron. These nine currents are shown in Table 2.1; the specific time and voltage dependence of each current is described in detail in (17]. 8
PAGE 16
Table 2.1: Nine variables describing ionic currents in Purkinje Fibers Type of Channel Variable Description POTASIUM 'tKz outward currents changing 'lx1 exponentially with step changes in membrane potential SODIUM ZNa fast outward current slow outward current .. .... CHLORIDE 'tqr variable activation and inactivation BACKGROUND CURRENT potasium outflow sodium influx during rest (Phase 4) 1:Na,b __ I. .. "..... 'Lcz,b chloride flow during plateau ''"The main equation which relates the interaction of these ionic currents to the overall membrane current is d.B =, dl; Cm (21) where E is the membrane potential, i1 is the sum of the nine Purkinje currents and Crn is the membrane capacitance. In a similar study conducted by Beeler and Reuter [1], the action potentials of myocardial cells were found by defining equations for five ionic channels. The ionic channels for myocardial cells are shown in Table 2 2, the specific equations for voltage and time dependency of the ionic exchanges can be found in [I]. 9
PAGE 17
Table 2.2: Five variables describing the ionic currents in myocardial cells. of Channel Variable I Description ] .. SODIUM Fast sodium channel Slow sodium channel POTASSIUM Time/voltage independant cunent ZK! "=Zxl Time/voltage independant cunent BACKGROUND Zextern Externally applied cunent during experiment The main equation relating the interactions of these ionic currents is similar to that for the Purkinje fibers and is given as dJ!J {iNa+ is+ iKI + ixl iextern) = .. (2.2) dt In [16] the action potential and excitation sequence in the ventricles are modeled by combining these two ionic current models. The model in [16] is a Jshaped two dimensional cut of the septum and left ventricle that has 730 cells (Figure 2.6). The Purkinje fibers are in the uppermost layer and the remaining nine layers are detlned as myocardial cells. The following set of differential equations is used to describe the ionic state of both Purkinje and myocardial cells located at coordinates (i, j), as well as the propagation of an action potential starting at a stimulus point and moving through the cells of the model [16]. Figure 2.6 Representation of the septum and the leji ventricle of the heart used in an ionic model. PoPurldnje fiber, muscle fJ 6}.
PAGE 18
dt dEi,[ dt t:,.i,j f) where f C' ' D f(Ei,j, Xi,j) + D !::,.i,j (change in membrane potential in Purkinje cells) G(Ei,j, Xi,j) (change in membrane potential in myocardial cells) Bi_ ..... _t,J + IEi!l,J + Ei,J 1 + .Ei,J+l Membrane potential of unit cell with coordinates (i, j) State variables for Purkinje and myocardial cells from tables 2.1 and 2.2 Functions describing relationship of Xi,j from eq. 2.1 and 2.2 respectively Electrical resistance of gap junction between cells Surface area of membrane of unit cell (larger in Purkinje) Membrane capacitance Coupling coefficient between unit cells. The coupling coefficients and the initial conditions for the variables are assigned numerical values estimated to correspond with experimental data [1)[16)[17]. The coupling coefi!cient D depends on the membrane capacitance Cm which is 12 times higher in Purkinje fibers than in myocardial cells. Thus, in order to simulate symmet rical interaction between cells, the value of D from Purkinje fibers to myocardial cells is set to 1/12 of the value of D from myocardial cells to Purkinje fibers. Simulation of one cycle of the entire model is done using numerical integration by Euler and Runga Kutta methods. The surface potential for the simulation is then calculated and an ECG of the simulation is formed. The ECG scans produced by the simulation cor responded closely to actual ECGs indicating that the ionic equations reasonably accurate in describing the action potential propagation. The main limitation of this model is the computational complexity of the numer ical integration of the equations. Because of this complexity it is only possible to II
PAGE 19
complete one cycle of simulation, which greatly restricts the possibilities for model simulation of activation phenomena of the heart. One way to avoid this computa tional complexity is to work with the essential features of the differential equations instead of the equations themselves. While this type of modeling may lose some of the accuracy and detail of the ionic model, it is much more efficient in modeling wholeheart phenomena and in analyzing the activation sequence over several periods. The remainder of this thesis will focus on ways to simulate the characteristics of the activation sequence without actually solving the differential equations which produce these characteristics.
PAGE 20
3. The Cellular Automaton Model A cellular automaton (CA) model is a set of discrete points which interact based on a given set of rules. The aim of the CA model of the heart is to reproduce accurately the depolarization and repolarization waves of the working heart by using a set of biologically based rules on an anatomically correct discrete representation of the cells of the myocardium. This is done by letting the discrete points of an integer grid represent the "cells" of the myocardium. The rules governing the interaction of these discrete points are derived from the end results of the equations used in the ionic models. Using the differential equations of the ionic model, the polarization state of each cell can be determined at any time step throughout the activation sequence. In a CA model, the underlying equations are ignored and only the timedependent states (repolarized or depolarized) of each cell are considered. Likewise, the ionic model can determine the cell which will depolarize in each time step in relation to the cells that have depolarized previously. Thus, since the ionic model determines that a depolarization wave will travel from one cell to a certain number of surrounding cells in 10 ms, the corresponding .rule for the CA is that this number of cells surrounding a depolarized cell will depolarize 10 ms or one time step after the initial cell. Thus the time step used in the model is 10 ms. The precise interactions between the cells of the CA model are described in the following sections. 3.1 Model Setup Models of the activation sequence of the heart were first used in the 1920s (20]. Since this time, the introduction of CA models has improved' the original models so that currently, CA models can produce very accurate simulations of the activation waves in the heart. The majority of these models are ventricular models only. Although some detail is lost by ignoring the atria, the nature of the activation
PAGE 21
process and the insulation of the ventricles from the atria allow ventricular models to retain many of the qualities of a full heart model while reducing the computational demands. The main differences among ventricular models is the accuracy of the geometrical description of the ventricles and the rules of interaction both within and between the cells. The main model we will be looking at was developed by Hermes [13] in order to analyze the interaction between ectopic and SA pulse firing. This model deals with a ventricular heart in which the ventricles are assumed to be electrically insulated from the atria. The epicardium is formed by the ellipsoid (x11)2 (y11)2 (z 16)2 100 + 100T 225 = 1. (3 .1) The endocardium is fanned by the paraboloid (x IW + (yn? = 4.9(z4j2 for z ::; 15. (3 2) The apex of the paraboloid is offcenter of the apex of the ellipsoid to make the left ventricular wall about three times thicker than the right. The septum is formed by the intersection of the inside of this paraboloid and the cylinder (3 3) The entire model lies in the first octant of a Cartesian coordinate system where 1 ::; :r; ::; 21, 1 ::; y ::; 21, and 1 ::; z ::; 15. A cross section of the model in the xz plane is shown in Figure 3 .1. These curves overlay an integer grid so that 2605 points with integer coordinates lie in the myocardial fiber between the ellipsoid and the paraboloid or in the septum; these points represent cells of the heart. Each cell is considered to be a neighbor of another cell if it has a Euclidean distance of V3 or less from another cell; thus each cell has at most 26 neighbors. The geometry is such that if a wave passes from a cell to its neighbors in one time step, the QRS complex will have a duration of 8 time steps (80 ms) when the heart is beating at 60 beats per minute [13]. Each cell has a record in a threedimensional array, whose indices correspond to the integer coordinates of the cell. In this record, cell type (epicardium, Purkinje,
PAGE 22
conduction system, etc.), the layer position of the cell, the polarization state of the cell and various other characteristics are recorded. The array and record for a sample cell is shown in Figure 3.1. Endocurdium ./ Cdl Posltln11 :',Se.Pl"um , , &. . a:::1' 'Epkl'lrjlum Cell Attributes Record p=U: io:::i : timi: after depolarization purk=l Purtclnje Fiber u;:;O rwl repuiMizlny r:::user set t 5: repolarlzation Figure 3.1 Cross section of" the CA model and the model setup for one cell. Each cell is represented by a position in an array of" records. The indices array are the coordinates of" the cell and the record keeps track t!t" the cell attributes. The general shape of this heart model is similar to most other CA heart models. The main differences between this and other ventricular models are the definition of a neighborhood of a cell and the number of cells in the model. In the MillerGeselowitz model (9] the general outline of the model is the same, but there are 50,000 cells. This increase in the number of cells leads to more realistic output and gives a better resolution of small changes in the depolarization and repolarization of the heart. In the Japanese model [11 ], the atria are included and the number of celis is 75,000. The cells in the Japanese model are also on a tilted plane so that the number of planar neighbors of each cell is increased from eight to twelve. Each cell in this model has 16 characteristics associated with it, 10 of which define the shape and duration of the action potential at any time. This provides a more accurate picture of the activation sequence because it allows partial depolarization to occur during the relative refractory period as was described in section 2.2. Also, the addition of the
PAGE 23
I atria in this model provides a better simulation of the depolarization wave from the SA node and allows for investigation of atrial anythmias and other disorders. 3.2 Depolarization Algorithm In a normal heart, the depolarization of the ventricles starts with the firing of the AV node. The pulse then travels down the bundle oflIis, separates into the right and left bundle branches and travels to the Purkinje fibers. In the ventricular model, there is no AV node so normal depolarization is started by firing six cells in the septum, placed experimentally to obtain the best agreement with the wave travel in the actual heart. In order to simulate the conduction system, the Purkinje fibers are set to fire four time steps after the firing of the initial cells in the septum. This makes the firing of the Purkinje fibers coincide with their firing if the depolarization wave traveled through the conduction system, as it does in an actual heart. The firing or depolarization of a cell with coordinates (i, j, k) is recorded by setting the polarization variable P(i,j,k) equal to 0 and a time variable t(i,j,k) equal to 1. In this model, the depolarization is isotropic so there is no preferred direction for the depolarization wave to traveL In the next time step, any of the 26 neighbors of the original cell that have not yet been depolarized (P(il,_jl, kl) = 1) may be depolarized. Thus in general, a cell C with coordinates (il, jl, kl) can be depolarized at a time step t if it is not yet depolarized and if its neighbor C' with coordinates (i,j, k) was depolarized in the previous time step. Once a cell is depolarized, t(i,j, k) is incremented by one for each time step after depolarization. When t(i, j, k) reaches the repolarization interval (r), the polarization variable P(i, j, k) is reset to 1 and the cell can be depolarized again. During depolarization, current is assumed to flow instantly from the cell depolarized in the last time step to the cell being depolarized in the current time step. The model also alJows for automaticity of ectopic cells located anywhere in the heart. The depolarization algorithm for the firing of these cells uses the celltocell propagation described above, but the depolarization wave spreads from the ectopic site instead of from the six targeted cells in the septum. Also, the Purkinje fibers are only depolarized through cell to cell propagation and are not manually activated four 16
PAGE 24
time steps after the firing of the ectopic focus. The model also allows for a combination of ectopic and SA firing so that waves can spread from two different starting points. In this case each of the waves will depolarize part of the heart and meet somewhere between the two starting points. The waves can overlap only if a cell depolarized by one wave has repolarized before being hit by the second wave. If a cell is already depolarized when it is hit by another depolarization wave, the wave will not spread to this cell. This interaction of the two waves is what causes many of the irregular heart beats, arrythmias and tachycardias that will be discussed in Chapter 5. In several CA models, the anisotropic nature of the heart muscles is modeled [9] [10] [ll] [15]. This is done by giving cells a preferred direction along the muscle fiber so that the depolarization wave wiil spread faster in that direction than across the muscle. The MillerGeselowitz model [9][1 0] approximates the anisotropic properties of the cells in each layer from endocardium to epicardium by using a block of anisotropic tissue. In this block there is an axis ofiow conductivity oriented normal to the epicardium which represents the slow flow of charge across muscle fiber. In all directions parallel to the epicardial surface there are axes of high conductivity representing the fast flow of charge along the muscle fiber. Thus a wave starting at a point within this anisotropic tissue spreads out in an elliptical wave front, depolarizing more of its neighbors within the same muscle fiber than in adjacent muscle fibers (Figure 3 .2). The magnitudes of the slow (v,) and the fast (vp) velocities for the elliptical wave are determined using (3.4) where 9it., get, gip, and 9ep are constants of intercellular and extracellular conductivity determined from experiments on the cardiac tissue. Using these values, the isochrones of the depolarization wave were found to match those of a working heart The values of v, and vP are constant throughout the heart. VI
PAGE 25
n "T .... .. .. + "''!I .... ... ........ .. .. ,. ..... ... ,, ' ... z .. ... .... 1 ,.. ... '.r > < T "' .,.,. +r ... ,.. ..... ... .... ' ........ . ' 'I .. .,.. .. .. : >rT ":r . ,.,. .,.. .. .. ............ ... ........... ... ...... ... ,, '.,. .... .. + ....... ......... ........ .. .. .. ,," @" .. ?_ ..,. .......... ,,. .,....,.._..... .,... Vr . ..... ' ' .. .......... ... .............. ... Figure 3.2 Elliptical wavejfonts on a cross section of a modified MillerGeselowitz model. A sample ellipse of activation is shown with the axis of slow and fast propagation flO}. In the atrial and ventricular Japanese heart model [11], anisotropy is also used in the depolarization algorithm. The method of determining the anisotropy in this model involves layering the heart like an onion and then twisting each layer a specific angle from apex to base. The elliptical wave fronts of the depolarization wave found with this model are similar to those of the MillerGeselowitz model. When comparing the global excitation sequences of isotropic and anisotropic models, the authors of the Japanese model found that there were no significant differences. This is attributed to both the high conduction velocities of the isotropic Purkinje fibers and to an averaging effect of the rotation of the fibers over the entire heart. The anisotropic heart did produce significant improvement over the isotropic heart in various experiments such as initiation of ventricular fibrillation. The algorithm used for incorporating anisotropy in either the MillerGeselowitz or the Japanese model could be added into our model in a similar manner The main advantage of incorporating anisotropy into our model would be to obtain a closer correlation between actual and simulated ECGs in the cardiac experiments which are discussed in Chapter 5. 3,2 Repolarization Algorithm As explained in Section 2.2, the action potential of a myocardial cell has five 18
PAGE 26
main phases. The third phase, rapid repolarization, is primarily responsible for the T wave on the ECG because most of the energy transfer between the interior and exterior of the cell occurs in this phase. However, the repolarization of a cell is not instantaneous and in fact occurs over a time period which is about two to three times as long as the duration of the Rwave. The Twave also has an absolute height of one third of the height of the Rwave. In our model, the repolarization wave is based on a parameter D which is the number of time steps required to depolarize all the cells in the heart; usually D is 8 to 10 time steps which corresponds to 80 to I 00 ms. Since the repolarization of a cell is gradual, in the model it is necessary to make each cell repolarize over several time steps instead of instantly. In order for the duration and height of the Twave to correspond correctly to that of the Rwave, the time needed for a cell to completely repolarize should be the same as the time required to depolarize the entire heart. Thus in each time step, a cell wili become repolarized and it will take D time steps for complete repolarization. If a cell depolarizes at time step t then it will begin repolarizing after the repolarization interval r or at time t + r. The repolarization interval corresponds to phase 1 through the beginning of phase 3 of the action potential and is usually 80 to 250 ms long. This interval represents the absolute refractory period and is set by the user. In the next section which deals with the Twave, the repolarization interval of each cell will depend both on the user set interval and on the location of the cell in the heart. Once a cell begins repolarizing, a repolarization parameter u is used to keep track of how long the cell has been repolatizing and when v. reaches 1 + D, the cell is totally repolarized. The value of u corresponds to the potential of the cell as was discussed in Chapter 2; so that v. = 1 corresponds to a potential of 0 m V and u = 10 corresponds to a potential of 90 mV. The height of the repolarization wave is recorded in a repolarization array hr(t) at each time step t. Thus, i( a target cell repolarizes in the current time step, the cont1ibution of this cell to the repolarization array is found by analyzing the potential of all of its neighbors. If a neighbor of the target cell repolarized in the previous time step then the potential of the neighboring
PAGE 27
cell will be 10 m V lower than that of the target cell. In order to simulate a IIIead ECG this change in potential is projected along the zaxis of the heart, so the contribution to the repolarization array is .10 times the difference in the zcomponents of the target cell and its neighbor cell. Neighbors of the target cell that did not repolarize in the previous time step do not contribute to the repolarization array. The total contribution of the target cell to the repolarization array in the current time step is found by computing the contribution for all of its neighboring cells. The difference in repolarization states between the target cell and its neighboring cells remains the same during the entire time necessary to repolarize a cell which is equal to D or about to time steps. Thus the total contribution of each cell is added into the repolarization array for 10 time steps after beginning of repolarization which occurs at the end of the repolarization interval for that cell. In Figure 3.3 a schematic diagram of the repolarization of the cell (10, 11, 4) at time step t = 34 is shown. The contribution to the repolarization wave from each of the neighbors of this cell in the xz plane is described. Neighborhood ol ce!i ai j11,1 iAJ z on crass section y:::ll 4 :1 Hi 1i L ll" CELL\ U jcn>rgel CoohiouiiVf< io I Reason 1 [J lhr{K] I I g h 140 .i difference !s r 140 .. 1 10anddlfferenceof .l toflrdinates in z, d!redlon !s 1. 1Jn 1 30 30 30 sam!': z value repolarized at same time step 1 taryet cell I Figure 3.3 Schematic representation of repolarization of a target cell located at (10, 11, 4). The v. values which tell how long a cell has been repolarizing are shown on the grid/or each ceil Each of the ceiis above the target ceil began repolarizing one time step before the target cell, so the charge contribution for each is.l(45) = .1. Using this algorithm with D set to 10, the Twave will be three times as long as the Rwave for a given depolarization speed. However, using this algorithm without 1 allowing for variation of the repolatization intervals produces a Twave which is of
PAGE 28
the correct duration, but is negative instead of positive. Many experiments have tried to resolve the fact that the T wave has the wrong sign and some of these will be explored in the next chapter. The repolarization algorithms used in the MillerGeselowitz [9][10] and Japanese [ 11] models are similar to this repolarization algorithm. However, as was described in Section 3.2, the Japanese model works with an anisotropic heart which is layered from epicardium to endocardium. In this model, the repolarization intervals for each cell are found by first determining the general repolarization interval for the entire layer and then varying it slightly for each individual celi based on its position in the layer. A modification of this technique is applied to our model in Chapter 4. 3.4 Output In most CA models, the output is in the form of an ECG which is an electronic record of the activity of the heart as described in Chapter 2. In our model, only the II Lead ECG is used although the other major leads can easily be implemented by projections. The depolarization wave values formed using the algorithm discussed in Section 3.2 are computed for each neighbor, C, of C' by forming a vector C'C = (i il,jjl,k kl). The sum of the projection of these vectors onto the zaxis which runs parallel to the II Lead is recorded as the pulse height at each time step I. Since this is an integer grid and since any neighbor of C' must be in the unit cube centered at C', the depolarization of each neighbor of C' adds (the distance in the z direction) to the pulse height which is recorded in the depolalization array, hd(t). In the repolarization algorithm, the portion of repolarization that occurs at each cell during one time step is recorded in the repolarization array, hr(t). Each cell con tlibutes i = to the array in any time step. The total pulse height at any time i is the sum of the depolalization and repolalization heights and is recorded in a third array h(t) = hr(t) + hd(t;). This pulse height is then graphed with corresponding time step to form the QRST complex of the ECG scan. Since there are no atria, there are no Pwaves in the ECG output of the model.
PAGE 29
4. Experiments on the Twave In a normal ECG recording, the Twave, which corresponds to the repolarization of the heart is upright, positive and has about one third the height and three times the length of the Rwave (Figure 4.1 ). The Rwave represents the depolarization of the ventricles in which the potential of the individual cells changes from 90 m V to 10 mV Since the depolarization wave travels down the ventricles, the cells in the top of the heart will be depolarized first causing the potential to be positive there while the nondepolarized cells in the bottom of the heart will cause the potential to be negative there. This difference accounts for the positive Rwave in the IT Lead ECG which measures the difference in potential along the zaxis or between the top and the bottom of the heart. If the repolarization wave is also assumed to travel down the heart, then the Twave which represents this wave will be negative. This is because during repolarization, the individual cells change from a potential of 10 m V to 90 m V Thus the potential top of the heart where the cells have begun repolarization and have a potential of 0 to mV, will be more negative than in the bottom of the heart where the cells have not begun repolarization and still have a potential of 10 mV Thus the positive polarity of both the Rwave and the Twave seems to indicate that the repoiarization wave does not travel down the heart like the depolarization wave does. In this section, four different modifications of the model which explore the spread of the repolarization wave and additional features of the myocardial cells will be used to attempt to explain the positive Twave. 22
PAGE 30
Fr j \ ,I ,j i\ li "'""" 10 Figure 4.1 ECGjor a normal S/l pulse without Pwaves. 4.1 Model A The first model is the nonmodified model described in Chapter 3. As noted before, if the repoiarization intervals of all ceils are the same, then the Twave of a normal SA pulse is inverted, as shown in Figure 4.2. This is because the first cells to repolarize will be those that depolarized first which are the six targeted cells in the septum. The repolarization wave, although prolonged, will spread in the same manner as the depolarization wave which is essentially down the septum from the endocardium to the epicardium. Since the wave spreads predominantly down the heart, at any given time step, most cells will have more cells with larger negative charge (more repolarized) above them than below them. Thus, when the II lead ECG is formed from the difference in potential between the top and the bottom of the heart, the Twave will be negative. Figure 3.3 shows a cell with coordinates (10, 11,4) and the potential of its neighbors at time step t = 34. The contribution of this cell to the height of the repolarization array hr(34) is 0.3 mV. In general, the total contribution of all the cells in a given time step will be negative, leading to the inverted Twave. Thus, using the same repolarization interval for all of the cells does not provide an accurate ECG output.
PAGE 31
; i!J msx10 Figure 4.2 ECG of an SA pulse for Model A The Twave is of correct duration and size. but it is inverted The minimum value of the is 18 mV 4.2 Model B In several cardiolOb'Y texts [3], it is suggested that the positive Twave is formed because of the rule that if a pulse originates from the SA node, then the first cells to depolarize will be the last cells to repolarize. Although this rule will account for the positive T wave, it is presented without an explanation in terms of the ionic processes of depolarization and repolarization. A cell in the bottom of the heart depolarizes because of an influx of sodium. In order to distinguish whether the depolarization wave affecting this cell originated at the SA node or from an ectopic focus, there must be a difference in either the amount of sodium in the extracellular fluid or in the amount of sodium allowed to pass into the cell based on the origination site. Since the CA does not include extracellular fluid or the ionic process, only the end results of this process can be incorporated into the model. The end result of the rule is that if the pulse originates from the SA node, the entries in the repolarization array which form the Twave should be positive. In the model, if a pulse originates in the SA node, the maximum height of the depolarization array hd( t) is generally much higher than if the pulse originates from an ectopic focus. Thus in order to implement the rule into the model, a trigger (s2) is set to I if the maximum value of the depolarization array is in the range of that of an SA pulse. Then, when the output of the model is 24
PAGE 32
formed, the entries of the repolarization array are: hr(t + j) .06hd(tr + 1) + hr(t + j) if .s2 = 1 (pulse from SA node) hr(t + j) = .llhd(t r + 1) + hr(t + j) if s2 = 0 (pulse from ectopic focus) This produces a positive Twave if the pulse originates tram the SA node and a negative Twave if the pulse originates from an ectopic focus. The values of .06 and .11 were chosen so as to make the output match that of a normal heart (13]. The ECG produced for a normal SA pulse in this model is the same as that of a working heart (Figure 4.1 ). The main problem with this model is in identifying an SA pulse. It is possible for some ectopic pulses to produce a maximum value of the depolarization array that is in the normal range for an SA pulse. If this occurs, the output for the model will be flawed. Also, the lack of biological support for this repolarization algorithm reduces the predictive capability of this model. In the next two models, the biological basis for the inverted Twave will be explored. 4.3 Model C Several studies on canine and human hearts have found that in fact the repo larization intervals are not the same for all the cells in the heart [ 4]. In a majority of the studies it has been found that the repolarization interval of ventricular cells is longer in the endocardium than in the epicardium. This repolarization variation was included in Model C by layering the ventiicles into parabolic shells representing the epicardium (outer shell), middle section, endocardium (inner shell) and septum. The repolarization intervals were then made the shortest in the epicardium and the septum (r 4 time steps) where r is the normal repolarization interval, slightly longer in the middle shell (r2 time steps) and the longest in the endocardjum (r). These intervals were picked after experimenting with several different interval times and finding that these worked as effectively as any other variation. Implementing the activation sequence with this repolarization method, the ECG for a standard SA pulse is shown in Figure 4.3.
PAGE 33
m&xlD Figure 4.3 ECCi fi>r a normal SA pulse finm Model C The is of correct duration, but it ts inverted The minimum value of the Twave is 10 mV which is 8 mV higher than that produced in Model A Although the Twave is still negative in this output, it is Jess negative and shorter than the Twave for model A (Figure 4.2). So, varying the repolarization intervals from epicardium to endocardium does affect the Twave, but it does not appear to be enough to make the Twave positive. Thus, it seems that there must be other factors that affect the repolarization process. 4.4 Model D In a study that considers the control of the action potential and how it relates to ionic features, Nobel [19] hypothesized that in order to obtain a positive Twave, the repolarization intervals should be shorter in the bottom than in the top of the ventricles. In order to implement this effect in Model D, repolarization intervals for each cell are based on epicardium to endocardium position as in Modet C and then this value is added to where z is the zcoordinate of the cell and ranges in value from 1 to 15. Thus if a cell is located in the epicardium and has coordinates (4, 5, 9) then, = = 4 and its repolarization interval is (T4 + 4) where r is the general repolarization interval set by the user The value is based on
PAGE 34
experimentation with model output. The repolarization intervals for an x z cross section of the model are shown in Figure 4.4. 3Q 28 28 28 28 28 30 28 28 24 28 28 28 28 28 28 28 28 27 31 27 27 27 27 27 31 27 27 27 27 31 27 27 27 27 27 31 27 27 27 26 28 30 26 26 26 26 26 30 28 26 26 26 26 26 30 26 26 26 26 26 30 26 26 26 26 25 27 29 25 25 25 25 25 29 27 25 25 25 25 25 29 25 25 25 251251 29 25 25 25 25 24 24@ 28 24 24 24 24 24 28 26 24 24 24 24 24 28 24 24 24 24 24 28 24 24 24 24 23 23 25 27 27 27 23 27 27 21 25 23 23 23 23 23 27 27 27 23 27 27 27 23 23 23 23 22 24 26 26 22 26 26 24 22 22 22 22 24 24 24 24 24 22 22 22 21 21 23 21 21 21 21 D This cell is in the septmn with a of 8. The repolalizationinterval is (user 0 This cell is in the middle layer with a of 7. The repolalization interval is1 (user Figure 4.4 fie polarization intervals on an xz cross section ofModel D. The intervals increase fiom bottom to top and from epicardium to endocardium The value of each interval c (user set value layer value 1 where the layer value is:epicardium and septum4. middle layer2, endocardiztm''O. The ECG for a standard SA pulse for model D is shown in Figure 4. 5. In this output, the Twave is of the correct height, length and direction and the ECG corresponds very well with that of an actual SA pulse (Figure 4.1 ). However, when an ectopic pulse is used in this model, the Twave remains positive instead of inverting as it would in a working heart. Thus, this model works for a normal pulse but does not explain the formation of the Twave for all possible activation sequences in the heart.
PAGE 35
mV rc' ji' Figure 4.5 ECG jiJr a normal SA pulse from Model D_ The Twave is of correct duration and size, and it is po:..'itive Of the four models, Model B and Model D both produce a correct ECG output for the standard SA pulse_ However, while the repolatization algorithm in Model B is merely contrived in order to make the Twave work, the hypothesis used to form the repolarization algorithm for Model D has a biological basis. The results of these experiments with the Twave indicate that further research should be conducted on the theories of apex to base repolarization interval variation and on determining the biological basis for the rule used in model B. 28
PAGE 36
5. Experiments As cardiac computer models become more complex and more accurate in simulating the functions of the heart, their use in cardiology is also increasing. A computer model can be used in cardiology in two ways. The first method is to use the model to evaluate a complex hypothesis about some biological phenomenon. By running the computer simulation several times it is possible to support a hypothesis or perhaps find cases in which the hypothesis cannot hold. The second method is to use the computer model to predict what will happen in a particular biological situation. This approach relies on the accuracy of the computer model not only in simulating what is known, but in containing enough biologically significant properties to make accurate predictions. Cellular automaton models have been used in both of these approaches with substantial success. This chapter will look at several such experiments using the models described in the previous chapters. 5.1 SelfSustained Ventricular Tachycardia (SSVT) Self sustained ventricular tachycardia (SSVT) can be caused when there is some type of block in the cells, such as scar tissue, that causes the depolarization wave to change its normal course. If the geometry of the block is right, the depolarization wave can actually circle around this blockage. This happens because cells that depolarize on one side of the blockage have completely repolarized by the time the depolarization wave hits them from the other side of the blockage [3]. The blockage causing SSVT does not have to be damaged tissue. If a group of cells are depolarized slightly before the depolarization wave from the SA node arrives, then these cells 'form a blockage because they cannot be depolarized again at that time. In the model, a single ectopic pulse applied prior to a single SA pulse can cause such a blockage and initiate SSVT In Model B, a blockage occurs when the ectopic focus at coordinates (16,11,9),
PAGE 37
which is in the interior of the left ventricle, is pulsed once at time step 12 and the SA node is pulsed once at time step 19 with a repolarization interval of 120 ms (Figure 5.1). The heart beat is irregular until the SA pulse and then it settles down into a regular pattern of ventricular tachycardia which continues without any stimulus from the SA node or an ectopic focus. In models A and C which have an inverted Twave for the SA pulse, SSVT can be initiated in a similar manner This suggests that even though the Twave for SA pulses is incorrect in these modeis, they are still effective in modeling essential features of the activation process. Figure 5.1 SSVTfor Model B one ectopic and one SA pulse. The pulse is erratic fhr the first 600 ms and then it settles into steady VI' without any odditional pulsing of either the SA node or the ectopic .fi1cus. In M.odel D the apex repolarization interval must be set to less than 120 ms in order to obtain a selfsustaining rhythm. The rhythm observed here is much more irregular than in the other models. Although the tachycardia continues without any additional stimulus, it does not seem to settle into a regular rhythm (Figure 5 .2). This suggests that the blocking factor in the interaction of the pulse from the ectopic focus and the SA node shifts around the ventricle instead of being confined to one region. Thus, the high variability of the repolarization interval seems to make it more difficult to stimulate the rhythmic SSVT behavior. The fact that Model D does not include the anisotropy of the myocardial fibers may also be a reason for the difficulties in obtaining rhythmic SSVT beats. In the Japanese model [11] introduced in Chapter 3, the authors found that in the anisotropic model it was easier to initiate ventricular
PAGE 38
tachycardia and ventricular fibrillation than in the isotropic model. I ,__j Figure 5.2 Nonrhythmic 5'S'VTfor Model D. Using the same parameters as in the for Model B. Model f) does not settle down into a rhythmic VT as Model B does. In all of the models it was possible to stop the SSVT and return the heart to a normal sinus rhythm using capture beats. A capture beat is either an external shock or rhythmic firing of the SA node that disrupts and stops abnormal beats. Using the ectopic focus and SA node interaction that produced SSVT in the previous example, a capture beat occurs when the SA node pulses rhythmically every 860 ms instead of just once. ln the ECG shown in Figure 5.3, the first SA firing, which occurs 90 ms after the firing of the ectopic focus causes SSVT The second SA firing occurring at 1000 ms disrupts the SSVT and the tachycardia ends at 1500 ms. The third SA firing occurring at 1860 ms is normal, showing that the SSVT has been captured and the heart is beating normally. This demonstrates that SSVT caused by the firing of a single ectopic pulse can be stopped by external (or SA node) impulses occurring at the correct time. However, application of these external impulses at the wrong time can further disrupt the beating process making the SSVT worse, sometimes causing ventricular fibrillation which is often fatal. The use of CA models allows for experimentation in order to understand when an external pulse should be applied.
PAGE 39
mV : '}$ msx10 Figure 5.3 ECG ofS;';'VT fix Model B with a capture beat occuring at 1860 ms. 5.2 Ventricular Tachycardia In addition to the abnormal rhythms that can be caused by just one ectopic pulse, ventricular tachycardia (VT) can also be caused by a constant rhythmic pulsing of an ectopic focus. In this type of v1: external pulses are still effective in arresting the abnormal beating, but do not stop the automaticity of the ectopic focus. Thus, the external pulse only interrupts the tachycardia which can start again if the ectopic focus resumes its automaticity. In order to prevent recurrent tachycardia, the automaticity of the ectopic focus must be permanently eliminated. This is done through a process called ablation in which the focus is located by probing the tissue of the heart with an electronic pulse from a surgically inserted catheter which initiates the automaticity of the abnormal focus when it is probed. Once the abnormal focus is found it is cauterized and thus can no longer fire. This procedure is effective in stopping the automaticity, but it is both timeconsuming and technically difficult to locate the focus through probing methods. A CA model can be used to help determine !he location of such a focus by experimentally choosing foci in the CA model and then seeing how the ECG output from each focus matches the ECG of the patient. To show that this type of differentiation is possible, VT is initiated from the 32
PAGE 40
upper left ventricle at (16,ll,ll) (Figure 5.4) and from the lower right ventricle at (5, 11,6) (Figure 5.5). These ECG scans are significantly different, showing the ability of a CA model to determine where to search for an ectopic focus. In addition it is possible that speciftc modifications could be made to the CA model to produce greater differentiability in the ECG scans of VT I L_'_ I Figure 5.4 Ventricular tachycardia initiated in the upper lefl ventricle. The Rwave has a maximum value of"l02 mV and the Twave has a minimum value ol mV r1 .. A 1 f\, r1 A ." .J ('\( i \I II II \: ::1,11 II II i II II i I I i _____ j Figure 5.5 Ventricular tachycardia initiated in the lower right ventricle. The Rwave has a maximum value
PAGE 41
5.3 WolffParkinsonWhite Syndrome (WPW) Another area where CA models can be used to increase the precision of cardiac surgery is in correcting WolffParkinsonWhite Syndrome (WPW). WPW is caused by conduction from the atria to the ventricles through an accessory pathway in addition to the conduction through the AV node. The abnormal depolarization cycle in the ventricles can lead to various disorders such as supraventricular tachycardia or ventricular tachycardia. WPW is very distinguishable on the ECG scan because the depolarization through the accessory pathway produces a delta wave which is a slurring in the initial portion of the QRS wave [3]. There are eight main areas in the base of the ventricles where these accessory pathways are usually formed [2] (Figure 5.6). As with automaticity of an ectopic focus, ablation is the main technique for eliminating these pathways. The ablation of the accessory pathway forces the depolarization wave to travel only through the AV node and thus returns the heart to normal beating. Figure 5.6 lt!teral slice across the top of the ventricles showing the eight main areas for accessory pathways in WPW In a study using a modified MillerGeselowitz model [10], WPW was successfully simulated in each of the eight sites shown in Figure 5.6. Since this model has no atria, the accessory pathway of WPW was simulated by initiating an ectopic pulse on the
PAGE 42
base of the ventricles in each of these sites. The ectopic pulse was fired approximately 40 ms before the start of normal excitation from the endocardium. The Lead I n and III ECGs for the simulation from site 8 is shown in Figure 5. 7. The polarities of the delta waves for all eight sites were then compared with each other and with actual patient data. The results of this experiment showed that the model was reasonably effective in predicting the site of the accessory pathway although the ECGs for some sites were closer to those of actual patients than other sites. This discrepancy can be explained partially by the sensitivity of the delta wave to the exact initiation site. The fact that the model has only 50,000 cells may make the propagation too coarse to reproduce the delta wave polarity accurately in all cases. Improvements both in the precision and in the depolarization algorithms of the model may help to eliminate these discrepancies. l Figure 5.7 fad I. II, and III ECCJs for WPW in a modified MillerCJeselowitz model. The !ad ll and Ill ECCJs show delta waves in the QRS complex [10}. 5.4 Right and l.eft Bundle Branch Block In addition to being able to simulate ectopic foci and accessory pathways, another feature of CA models is the ability to simulate ischemia or damaged tissue and conduction system disorders. One very common type of conduction 'system disorder is right bundle branch block which occurs when damage to the conduction system stops the depolarization wave from traveling down the right bundle branch to the Purkinje fibers. This means that the right ventricle will be depolarized after the left since the depolarization wave cannot reach the right ventricle directly, but must first
PAGE 43
travel through the bottom of the heart. This disorder is recognizable on the ECG by a prolonged QRS complex which indicates the additional time required to depolarize the right ventricle. Since this is a depolarization disorder, and it is not affected by the repolarization sequence, Models A, B, C, and D are all able to simulate this disorder correctly. In Figure 5.8 the normal SA pulse ECG for model D is shown along with the right bundle branch block ECG which has a QRS that is approximately 10 ms longer than in the normal ECG. The CA model can also be used to simulate other ischemic disorders such as scar tissue caused by a heart attack or other conduction system blocks [12]. ' mV ( 'Ulr:J '"' Fi!l"re 5.8 ECG for simulated right bundle branch block and a normal SA pulse in Model D. The QRS of the nght bundle branch block ECG is 10 ms longer than for the SA pulse.
PAGE 44
6. Conclusions The CA models explored in this thesis are effective in modeling many of the functions of the heart. The models can be used to explore the interaction of pulses originating at an ectopic focus, producing an ECG output that simulates that of an actual heart. The models correctly reproduce ECGs for the depolarization of the heart and provide insight into the repolarization process. By using a CA model it is possible to understand exactly how the depolarization and repolarization waves are formed based on the travel of action potentials through the cells. This understanding provides a tool for improving analysis of ECG output. The experiments on the inverted Twave show how CA models can be used to support and encourage research in this area. In Chapter 5, two of the models pre sented suggest specific approaches to analyzing the repolarization intervals in order to determine why the Twave is inverted. Resolving the Twave question for CA models would also answer many other unresolved questions about anatomy and cel lular interactions. If theCA produced a working Twave based on some set of logical rules, then it is possible that the biology supporting these rules may be applicable to the heart. The use of CA models as a diagnostic tool in disorders such as selfsustained ven tricular tachycardia, ventricular tachycardia, WolffParkinsonWhite syndrome, right bundle branch block and others may prove to be very useful. The CA can be a powerful tool in improving surgical procedures by eliminating error and increasing success in finding the source of disorders through simulations and comparisons of model and actual ECG output. It is possible that the user interface of this model could be improved so as to make it a viable model for doctors to use in diagnosis and correction of cardiac disease. In addition, the CA model is a valuable tool for learning and experimenting with the activation sequence and cellular interactions in
PAGE 45
the heart. There are many areas in which the CA models can be improved. The accuracy of the activation sequence can be increased by including more of the essential character istics of the action potential such as a relative refractory period. The number of cells in the CA model can be increased to produce a finer resolution on the ECG output, and provide more variability in modeling applications Finally, with the increasing power of computers it is possible to combine the ionic and CA models to produce an accurate and efficient model which can precisely model the ionic interactions on a cellular level. 38
PAGE 46
7. References 1. G. Beeler, H. Reuter, Reconstruction of ventricular myocardia/fibers, Journal of Physiology, Vol. 268, (1977), pp. 177210. 2. D. Benson D, J. Gallagher, M. Spach, R. Sterba, wcalization of the site of ventricular preexcitation with body surface maps in patients with WolffParkinson White Syndrome, Circulation, Vol. 65, No.6, (1982), pp. 12591268. 3. R. Berne, M. Levy, Cardiovascular Physiology, Mosby Year Book, (1992). 4. M. Burgess, P. Ershler, K. Spitzer, B. Steinhaus, Nonuniform epicatdial activation and &polarization properties of in vivo canine pulmonwy conus, Circulation Research, Vol. 62, No. 2, (1988), pp. 233246. 5. A. Camm, M. Malik, Cardiac electrophysiological experiments in numero, parts 11111, PACE, VoL 14, (1991), pp. 16481671 and 21672186. 6. A. Camm, M. Malik, Computer simulation of electronic interactions during exci tation and repolarization (![myocardial tissue, Medical and Biological Engineering, (1991), pp. 425432. 7. D. Davis, D!!forential Diagnosis of Arrhythmias, WB Saunders Co, (1991). 8. A. Despopoulos, S. Silbemagl, Color Atlas of Physiology, Thieme Medical Pub lishers, (1991) 9. D. Geselowitz, W. Miller, Simulation studies of the electrocardiogram. I. lhe normal heart, Circulation Research, Vol. 43, (1978), pp. 301315. 10. R. Gulrajani, M. Lorange, Computer simulation of the Wo!jfParkinsonWhite preexcitation syndrome with a modified MillerGeselowitz heart model, IEEE Transactions on Biomedical Engineering, VoL BME33, No. 9, PP 852873. :l9
PAGE 47
11. E. Harasawa, K Harumi, H Hosaka, D. Okazaki, D. Wei, Comparative simulation of excitation and body surface electrocardiogram with isotropic and anisotropic computer heart models, IEEE Transactions on Biomedical Engineering, VoL 42, No. 4, (1995), pp. 343357. 12. K Harumi, T Musha, G. Nishiyama, Y Okamoto, H. Tusunakawa, D. Wei, G. Yamada, Clinical application of electrocardiographic computer model, Journal of Electrocardiology, VoL 22 Supp., (1989), pp. 5436. 13. H. Hermes, A computer modeljbr analysis and control of ventricular arrhythmias, Unpublished, (1995). 14. A. Hodgkin, A. Huxley, A quantitative description of membrane current and its applications to conduction and excitation in nerve, Journal of Physiology, Vol. 117, (1952), pp. 500544. 15. B. Horacek, L. Leon, Computer model of excitation and recovery in the anisotropic myocardium, Journal of Electrocardiology, Vol. 24, (1991), pp. 115. 16. M. Kawato, K Okazaki, R Suzuki, S. Urushibara, A. Yamanaka, Reconstruction ()[electrocardiogram using ionic current models for heart muscles, Japanese Heart Journal, VoL 27 Supp., (1986), pp. 185193. 17. R. McAllister, D. Nobel, R. Tsien, Reconstruction of electrical activity (y, VoL 251, (1975), pp. 159. 18. D. Nobel, A modification (
