ABSTRACT
Title of Dissertation: VLSI DESIGN OF HEART MODEL
Li Wang, Doctor of Philosophy, 2007
Dissertation directed by: Professor Robert W. Newcomb
Department of Electrical and Computer Engineering
Heart disease is a leading cause of death in the United States and abroad.
Research interests arise in understanding the nature of the dynamics of the heart
and seeking methods to control and suppress arrhythmias. Simulation of the heart
electrical activity is a useful approach to study the heart because it yields some
quantities of interest that cannot practically be obtained in any other way. How-
ever, the complexity of the human heart leads to complicated mathematical mod-
els, and consequently, modeling arrhythmias of a whole heart with computers is
extremely data intensive and computational challenging.
In this dissertation, we introduce an analog VLSI design that simulates cardiac
electrical activities. The selected cardiac model is based on the Beeler-Reuter
equations and the continuous core-conductor model. The Beeler-Reuter equations
formulate the membrane ionic kinetics of ventricular cells, and the core-conductor
model describes the electrical signal conduction on cardiac tissues. We discuss
the design  ows of mapping equations into circuits and present a set of circuit
blocks of basic mathematical function units. The transistor circuits for realizing
the ionic model of a single cell is introduced, and capacitors are used to calculate
time directives. A method of shifting the initial conditions of difierential equations
to zero is discussed for saving the circuit which sets up the initial voltages of
the capacitors. We also introduce a method of implementing reaction-difiusion
systems using non-linear RC networks, and present the circuit which simulates the
reaction-difiusion process, i.e. the electrical propagation, of the heart.
Erroranalysisiscarriedoutforthecircuit-realizedBeeler-Reutermodelbycom-
paring the simulated functions with the equation calculated values. The PSpice
simulation results show that the circuit created action potential is satisfactory. The
important reentry phenomena, the primary mechanism underlying flbrillation, is
presented, and an anatomical reentry in the 1-dimensional model and a functional
reentry (spiral wave) in the 2-dimensional model are successfully simulated in cir-
cuits.
The presented methods of implementing equations with analog VLSI circuit
contribute to the fundamentals for a novel technique of obtaining numerical solu-
tions and potential fast application-specifled analog computational devices if the
circuits are fabricated on chips. Unlike computing with digital computers, which is
mainly a serial process and needs to discretize the space and the time domain for
flnding numerical solutions of the discretization points one by one, computation
with analog VLSI relies on the physics of the electrical devices and takes advantage
of the integration properties of capacitors and, hence, computing in analog circuit
hardware is a parallel process and can be real-time, that is, the calculation time is
the time simulated by equations.
VLSI DESIGN OF HEART MODEL
by
Li Wang
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulflllment
of the requirements for the degree of
Doctor of Philosophy
2007
Advisory Committee:
Professor Robert W. Newcomb, Chairman
Professor Pamela Abshire
Professor William E. Kirwan
Professor Martin Peckerar
Professor Gang Qu
c Copyright by
Li Wang
2007
DEDICATION
To the people who can be potentially beneflted by this dissertation.
ii
ACKNOWLEDGEMENTS
I would like to thank the following people who made this work a reality. First,
I would like to thank my advisor, Professor Newcomb, for his constant and patient
guidance. He was always on campus hard-working on Sunday, and this made it
possible for me to discuss with him frequently and get quick suggestions when I
full-time worked ofi the campus. Professor Newcomb often read my dissertation in
his fully occupied schedule by scarifying his personal time at home and gave me
immediate helpful feedback. His valuable knowledge in the academic fleld and his
dedicated attitude to science not only made my graduate study very fruitful but
also gave me a life-model for my future career.
I would also like to thank my committee members, Dr. Abshire, Dr. Peckerar,
and Dr. Qu, for kindly consenting to join the defense and review this disserta-
tion; special thanks to Dr. Kirwan for managing to be on the committee from
his extremely heavy business for the University. I would also like to thank Ms.
Xinhua He, who gave me helpful discussions, and Ms. Yu Jiang, Ms. Koranan
Limpaphayom, and Ms. Yingying Zhao for contributing a friendly and pleasant
working environment in the Microelectronics Lab.
Finally, I would like to thank my families - my parents and sister, and special
thanks to my husband, Zhu Han, for having been always by my side when I had
di?culties, giving me great encouragement, and supporting me throughout these
years.
iii
TABLE OF CONTENTS
List of Tables vi
List of Figures vii
1 Introduction 1
1.1 Modeling Heart Electrical Activity . . . . . . . . . . . . . . . . . . 1
1.2 Review of Mathematical Models for Cardiac Cell Electrical Activity 4
1.3 Motivation for VLSI Implementation of Heart Model . . . . . . . . 8
1.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Organization of Dissertation . . . . . . . . . . . . . . . . . . . . . . 13
2 Heart Electrophysiology and Its Mathematical Models 15
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Cellular Electrophysiology . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.1 Cell Membrane . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.2 Ionic Mechanism of Action Potential . . . . . . . . . . . . . 20
2.3 Mathematical Description of Cardiac Electrophysiology . . . . . . . 24
2.3.1 Voltage Clamp . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.2 Ionic Membrane Models . . . . . . . . . . . . . . . . . . . . 26
2.3.3 Beeler-Reuter Model . . . . . . . . . . . . . . . . . . . . . . 35
3 VLSI Design of the Beeler-Reuter Model 42
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2 Model Reformulation . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2.1 Transformation of Mathematical Description . . . . . . . . . 44
3.2.2 Parameter Scaling for VLSI Design . . . . . . . . . . . . . . 48
3.2.3 Initial Value Shifting . . . . . . . . . . . . . . . . . . . . . . 51
3.2.4 Reformulated Beeler-Reuter Equations for VLSI Design . . . 54
3.3 Circuit Blocks of Function Units . . . . . . . . . . . . . . . . . . . . 54
3.3.1 Linear Function Circuits . . . . . . . . . . . . . . . . . . . . 59
3.3.2 Exponential Functions . . . . . . . . . . . . . . . . . . . . . 72
3.3.3 Multipliers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.4 VLSI Design of Beeler-Reuter Model . . . . . . . . . . . . . . . . . 89
iv
3.4.1 Time-Independent Potassium Current IK1 . . . . . . . . . . 92
3.4.2 Gating Variables . . . . . . . . . . . . . . . . . . . . . . . . 94
3.4.3 Time-Dependent Potassium Current Ix1 . . . . . . . . . . . 98
3.4.4 Time-Dependent Inward Sodium Current INa . . . . . . . . 103
3.4.5 Time-Dependent Slow Inward Calcium Current Is . . . . . . 109
3.4.6 Action Potential Vm . . . . . . . . . . . . . . . . . . . . . . 114
3.4.7 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4 Propagation of Electrical Activity in Cardiac Tissue 139
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.2 Electrical Propagation Model of the Heart . . . . . . . . . . . . . . 141
4.2.1 1-Dimension Cable Model in Cardiac tissue . . . . . . . . . . 142
4.2.2 Multi-Dimensional Cardiac Propagation Model . . . . . . . . 145
4.2.3 Summary of Cardiac Propagation Model . . . . . . . . . . . 146
4.3 Modeling Propagation of Cardiac Active Potential with VLSI . . . . 148
4.3.1 ComputationofDiscretizedReaction-DifiusionEquationwith
Circuits in Cartesian Coordinates . . . . . . . . . . . . . . . 150
4.3.2 TransistorCircuitSimulationof1-DimensionalCardiacProp-
agation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
4.3.3 2-Dimensional Cardiac Propagation: Spiral Reentry . . . . . 159
5 Conclusions and Future Work 168
5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
5.2.1 More - Application Perspectives . . . . . . . . . . . . . . . . 173
A Transistor Circuits of Beeler-Reuter Model Implementation 178
B PSpice Simulation of Anatomical Reentry Using Unidirectional
Block 213
Bibliography 219
v
LIST OF TABLES
2.1 Typical ion concentrations and equilibrium potentials in cardiac cell
[24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Units used in the Hodgkin-Huxley equations . . . . . . . . . . . . . 34
2.3 Initial conditions for the Beeler-Reuter Model . . . . . . . . . . . . 41
3.1 Statistics of goodness of flt for equations (3.5) and (3.6) . . . . . . . 47
3.2 Relations of original variables and the scaled and shifted variables,
and their electrical representative in circuits. . . . . . . . . . . . . . 55
3.3 Summary of function circuits. Ci (i = 1;2:::N) are constant. . . . . 56
3.4 Vbias for initial conditions of gating variable capacitors (unit: V). . . 97
3.5 Equations of the opening and closing rate variables for INa. . . . . . 103
3.6 Capacitances in the VLSI realization of Beeler-Reuter model (unit:
F). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
3.7 Monte Carlo analysis of threshold voltage variation on circuits . . . 131
A.1 List of schematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
vi
LIST OF FIGURES
1.1 Cross section of the heart [3]. . . . . . . . . . . . . . . . . . . . . . 2
2.1 Membrane potential Vm is the intracellular electrical potential with
respect to the extracellular media. . . . . . . . . . . . . . . . . . . . 16
2.2 Variable dependence in the Beeler-Reuter model. . . . . . . . . . . 17
2.3 Responses of membrane potential to various stimulus. (a)Four stim-
ulating currents with difierent strength. (b)Corresponding mem-
brane potentials in response of the stimulus. An action potential is
created when the membrane potential is charged by the 4th stimulus
to pass a threshold. . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Phases of ventricular action potential and active ionic currents in
each phase. (0)upstroke, (1)early re-polarization, (2)plateau, (3)late
re-polarization, (4)rest. . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.5 Measurement of ionic current using voltage clamp technique. . . . . 25
2.6 Equivalent circuits for Hodgkin-Huxley type of membrane models.
(a)Membrane is approximated by a circuit with a capacitor and
ionic currents. (b)Each ionic current represents a type of ionic  ow
and can be modeled as a nonlinear resistor connected with a Nernst
potential. RA = 1=(gAy). . . . . . . . . . . . . . . . . . . . . . . . . 27
2.7 Equivalent circuits for the squid axon model of Hodgkin & Huxley.
RNa = 1=(gNam3h); RK = 1=(gKn4); and Rl = 1=gl. . . . . . . . . . 32
2.8 Simulated ventricular action potential using the Beeler-Reuter model. 36
3.1 Design  ow of VLSI Realization. . . . . . . . . . . . . . . . . . . . . 43
3.2 Plot of transformed Ik1 equation and its original equation. . . . . . 47
3.3 Plot of transformed fim equation and its original equation. . . . . . 47
3.4 Initial value shifting for difierential equations. . . . . . . . . . . . . 53
3.5 CMOS difierential pair works as a linear voltage to current con-
verter. (a)Difierential pair circuit. (b)Altered difierential pair that
works as a linear resistor. . . . . . . . . . . . . . . . . . . . . . . . . 60
3.6 Linear current to voltage converter. . . . . . . . . . . . . . . . . . . 62
vii
3.7 Voltage bufier. (a)Traditional voltage bufier. (b)Proposed voltage
bufier for the heart implementation. (c)Another version of proposed
voltage bufier, modifled for reducing drain-source voltage difierence
of M1 and M2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.8 Voltage subtractor. (a)Maundy voltage subtractor. (b)Improved
voltage subtractor. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.9 Pool circuit works as voltage adder/subtractor. . . . . . . . . . . . 70
3.10 Improved pool circuit. . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.11 Emitter-coupled circuit. . . . . . . . . . . . . . . . . . . . . . . . . 73
3.12 Circuit that implements equation (3.49). . . . . . . . . . . . . . . . 74
3.13 Circuit implementation of equation (3.60). . . . . . . . . . . . . . . 76
3.14 Circuit implementation of logarithm function. . . . . . . . . . . . . 79
3.15 Tannocurrentmultiplier. (a)BasicblockofTannomultiplier. (b)Tanno
multiplier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.16 Current subtractor. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.17 Multiplier based on NPN transistors. . . . . . . . . . . . . . . . . . 85
3.18 Current amplifler/shrinker based on bipolar transistors. . . . . . . . 88
3.19 Top-level circuit diagram of the Beeler-Reuter model. . . . . . . . . 90
3.20 Diagrams of ionic current modules: (a)IK1, (b)Ix1, (c)INa, (d)Is. . . 91
3.21 IK1 moduleiscomposedofaconstantcurrentsourceandtwosigmoid-
function circuits. Refer to Appendix A, schematics page Sch-2 to
Sch-4 for transistor circuits. . . . . . . . . . . . . . . . . . . . . . . 92
3.22 Ideal and simulated curves: IK1 vs. Vm. . . . . . . . . . . . . . . . . 93
3.23 Detailed structure of gating variable circuits x1, m, h, j, d, and f. . 95
3.24 Ideal and simulated opening and closing rate of x1. (a)fix1, (b)flx1. . 99
3.25 Ideal and simulated curves: Ix1 vs. Vm. . . . . . . . . . . . . . . . . 100
3.26 PMOS current multiplier based on Tanno multiplier. . . . . . . . . 101
3.27 Simulation results of modifled Tanno multiplier. . . . . . . . . . . . 102
3.28 Ideal and simulated opening and closing rates of gating variables for
INa. (a)fim. (b)flm. (c)fih. (d)flh. (e)fij. (f)flj. . . . . . . . . . . . . 105
3.29 Circuit diagram of INa. Refer to Appendix A, schematics page Sch-
10 for transistor circuits. . . . . . . . . . . . . . . . . . . . . . . . . 107
3.30 Circuit of a multi-input adder/subtractor for INa module. Refer to
Appendix A, schematics page Sch-17 for transistor circuits. . . . . . 108
3.31 Ideal and simulated opening and closing rate variables for Is. (a)fid.
(b)fld. (c)fif. (d)flf. . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
3.32 Circuit diagram of Is. Refer to Appendix A, schematics page Sch-19
for transistor circuits. . . . . . . . . . . . . . . . . . . . . . . . . . . 112
3.33 Ideal and simulated action potential. . . . . . . . . . . . . . . . . . 115
3.34 Idealandsimulatedmembraneioniccurrents. (a)IK1. (b)Ix1. (c)INa.
(d)Is. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
viii
3.35 Comparison among experimental data, extracted equations, and cir-
cuit simulation results. . . . . . . . . . . . . . . . . . . . . . . . . . 118
3.36 Experimental data [14] (dots) and its fltting curve of steady-state
current-voltage relation of outward current in purkinje flbers. . . . . 119
3.37 [48] (a) Steady-state activation d1 of calcium current Is; data ob-
tained from flve cow ventricular trabeculae superfused with solution
containing 0.2(?), 0.45(?) or 1.8 mM CaCl2 (? ? 24). (b) Rate
constant for decay 1=?d, measured with solution containing 0.2(?),
0.45(?) or 1.8 mM CaCl2 (?). . . . . . . . . . . . . . . . . . . . . . 120
3.38 Current amplifler/shrinker circuit used for amplifying capacitance. . 124
3.39 Current amplifler/shrinker circuit using NPN transistors. . . . . . . 125
3.40 Simulation results of using current shrinker circuit to amplify ca-
pacitance. (a)Voltages of original capacitor and amplifled capacitor.
(b)Currents of original capacitor and amplifled capacitor. . . . . . . 126
3.41 Cascode structure can be used to minimize the Early efiect on Qn4
and improve the accuracy of Iout. . . . . . . . . . . . . . . . . . . . 127
3.42 Histogram and Gaussian flt for the NMOS threshold voltages mea-
sured from 62 sets of data (AMI 1.5 ?m technology). . . . . . . . . 130
3.43 Monte Carlo analysis of Current-Voltage Converter #1 (Figure 3.6). 131
3.44 Proposed linear current to voltage converter. (a) Circuit of linear
current to voltage converter. (b) Monte Carlo analysis: 20 runs. . . 133
3.45 Difierence of Iout in sigmoid function circuit caused by variation of
the threshold voltage. ?Vrandom=VT = 0:3. . . . . . . . . . . . . . . 135
3.46 (a)ModifledemittercoupledpairusingcascodeNPNdevices; (b)Cascode
NPN structure used in exponential function circuit. . . . . . . . . . 138
4.1 1-Dimensional cable model for cardiac tissue. (a)A 6-cell cardiac
muscle flber (top) and its representative (bottom) in cable model -
a continuous media that does not consider the boundaries of cells.
(b)Circuit diagram of cable model. . . . . . . . . . . . . . . . . . . 142
4.2 Simplifled structure of cardiac tissue and its continuous representa-
tive in cable-model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
4.3 Using an RC network to implement reaction-difiusion systems. . . . 151
4.4 Using an RC network to implement reaction-difiusion systems in the
case of (a) 1-dimensional, (b) 2-dimensional. . . . . . . . . . . . . . 154
4.5 Circuit of 1-dimensional cardiac propagation. . . . . . . . . . . . . . 156
4.6 PSpice simulation results of 1-dimensional cardiac propagation con-
structed by transistor circuit. . . . . . . . . . . . . . . . . . . . . . 158
4.7 Action potentialapplied with asecond stimulusatdifierentmoment:
(1) t = 200 ms, (2) t = 300 ms, (3) t = 350 ms, (4) t = 380 ms, (5)
t = 400 ms, (6) t = 410 ms, and (7) t = 420 ms. Stimulus strength:
300 mA=cm2 for 0.1 ms. . . . . . . . . . . . . . . . . . . . . . . . . 161
ix
4.8 Example of reentry. (a) Action potential travels normally. (b) Reen-
try is caused by a unidirectional block. . . . . . . . . . . . . . . . . 162
4.9 Circuit for 2-dimensional cardiac propagation. . . . . . . . . . . . 163
4.10 A spiral wave in 2-dimension cardiac model. . . . . . . . . . . . . . 167
A.1 Legends for hierarchical schematics. . . . . . . . . . . . . . . . . . . 179
B.1 A circuit ring composed of 120 node for modeling anatomical reentry.214
B.2 Implementation of bidirectional and unidirectional resistors using a
modifled difierential pair circuit. . . . . . . . . . . . . . . . . . . . . 215
B.3 Reentrywavecausedbyunidirectionalblockin1-dimensionalmodel.
217
B.4 The time-course of the membrane potential of node #1. . . . . . . 218
x
Chapter 1
Introduction
1.1 Modeling Heart Electrical Activity
The heart is one of the most essential organs in our body. The heart working
together with the blood vessels transports blood in the body and maintains the
circulation of our life. In this system the heart works as a reliable rhythmic pump
and as a blood reservoir. The heartbeat is so critical to our life and any interruption
in the heartbeat for more than a few minutes can cause circulatory collapse and
death. Heart disease is a leading cause of death in the United States and abroad.
In 2002, a total of 696,947 people died of heart disease in the United States[1],
accounting for 29% of all U.S. deaths. In 2005, heart disease is projected to have
cost $393 billion, including health care services, medications, and lost productivity.
Therefore, great research interests have arisen around people in investigating the
factors that cause and sustain these life-threatening conditions.
The pumping mechanical activity is driven by electrical activity of the heart.
The study of cardiac electrical excitation is a necessary and important step for
understanding the heart?s electrophysiology and its mechanical deformation to as-
1
Figure 1.1: Cross section of the heart [3].
sist investigating the methodologies of controlling and suppressing heart disease.
Shown in Figure 1.1 is a cross section of a heart. The heart is separated into
two similar structures, the right and left halves, which represents two functionally
difierent parts in the blood circulation system [2](pp. 98). The right half collects
the deoxygenated blood from the body and pumps it to the lungs. The left half
receives the oxygenated blood from the lungs and delivers it to the body. The
halves can be further divided into an upper and lower part, called atria and ven-
tricles respectively. The atria and ventricles are composed of walls surrounding a
cavity, which is normally fllled with blood. The walls consist primarily of a muscle
structure, the myocardium 1. The contraction of the myocardium causes the blood
 ow. The electrophysiology is tightly coupled with the pump function of the heart
1Myocardium is the muscular wall of the heart[4].
2
by controlling the development of tension. The origin of the electrical activity of
the heart are the cardiac cells, which show electrical excitability like nerve cells.
The electrical excitation of a cardiac cell is tightly coupled with its mechanical con-
traction. The electrical excitation is propagated from a cell to neighboring cells
via gap junctions by intercellular transport of ions.
The study on live models for cardiac electrophysiology, though fruitful and
necessary, have only provided partial insights into the mechanisms of arrhythmia
formation and termination [5]. One reason for this is that all quantities of inter-
est cannot be measured in vivo. An experimental investigation can be limited by
coarse resolutions, hard data collection during the critical time following shocks,
easy damage to tissue, and di?cult optical mapping for cells not on the surface. To
integrate and interpret the measurements, computational modeling of the system
can contribute greatly to our understanding of the heart. The simulation of the
heart yields information that cannot practically be obtained in any other way, and
makes it possible to predict prior unknown behavior through complex phenomena.
The simulation of the heart allows us to gain knowledge of the mechanisms of
heart failure, to simplify the development and validation of heart drugs and med-
ical devices, and to explore the side efiects of a product. Modeling of the heart
is an interdisciplinary research activity that includes the areas of anatomy, elec-
trophysiology, excitation propagation, force development and mechanics as well
as the coupling of these areas [5]. In this dissertation, we focus on modeling the
electrophysiology of the heart. We introduce a novel methodology of simulation
which utilizes analog VLSI circuits to model electrical activities of the heart. Com-
pared to running simulation of a heart model on digital computers, the introduced
method can be implemented on analog hardware which is able to obtain simulation
3
results fast. We present a set of circuits components that are useful for implement-
ing basic mathematical expressions in membrane ionic models, and perform the
PSpice simulation for the VLSI circuits which realize the Beeler-Reuter equations.
1.2 Review of Mathematical Models for Cardiac
Cell Electrical Activity
A large number of mathematical models have been constructed to describe the
electrophysiological behavior of cardiac cells. External electrical stimulus to a car-
diac cell may lead to a changing membrane potential, known as anaction potential.
The models that describe action potentials and the ion transport across the cell
membrane are also referred to as ionic models. In the following an overview of car-
diac electrophysiological models are presented. The principles and the advances of
the models are discussed.
Most ionic models of cardiac cells are based on the Hodgkin-Huxley formalism
[6], which flrst reported a quantitative description of the active and passive elec-
trical behavior of the axon membrane of giant squids [7]. The model was built
in 1952, a work that resulted in a Noble Prize for the authors. It described the
electrophysiology of the axon membrane using a RC (Resistor-Capacitor) circuit,
in which the resistor is signiflcantly nonlinear. The nonlinearities can be attributed
to the behavior of ionic channels [5]. Hodgkin and Huxley separated the membrane
current into three types of ionic current components: sodium ion  ow INa, potas-
sium ion  ow Ik and a small \leakage current" Il. The currents were characterized
as the product of conductances and the difierences of the driving forces for ions,
namely chemical gradient and electrical gradient. The conductances for INa and Ik
4
are time- and voltage-dependent, and modeled by using gating variables, which are
formulated by flrst-order ODEs (Ordinary Difierential Equations), and make the
whole model a system of degree-four difierential equations, i.e. a time-dependent
membrane potential and three gate variables. The Hodgkin and Huxley model
described the way in which ionic currents vary with membrane potential and time.
Its structure forms a basis for almost all models of excitable membrane behavior.
The FitzHugh-Nagumo model reduces the Hodgkin-Huxley model from degree-
four difierential equations to degree-two, for which phase plane analysis applies.
The FitzHugh-Nagumo model [8][9] was a modiflcation of the van der Pol equation .
The Van der Pol equation was intended to represent the qualitative rather than the
quantative properties of a wide class of excitable-oscillatory systems using simple
algebraic form [10]. The resulted BVP model (B. van der Pol model) has two
variables of state, representing excitability and refractoriness, and its properties
can therefore be visualized on a phase plane. For example, Glaze utilized BVP
model to perform simulations for heart in [11]. The FitzHugh-Nagumo model
used a single mathematical processes to represent multiple channel properties and
did not target to model the accuracy and complexity of cellular processes. It
preserved the essential behavior of the membrane and is still used, due to its low
computational cost, as the excitability component in some models of cardiac action
potential propagation [12].
In 1962, Noble published the formulation applicable to Purkinje flbers of the
heart [13]. It is one of the flrst mathematical models of a cardiac cell. The Purk-
inje flbers difier from squid nerve in that depolarization decreases the potassium
permeability of the membrane. During large depolarization, this decrease is tran-
sient and the potassium permeability increases slowly during the passage of the
5
depolarizing current. By modifying the Hodgkin-Huxley equations to take account
of this behavior in potassium current, Noble described the long lasting action and
pacemaker potentials of the Purkinje flbers of the heart. A novel quality of this
model was the decomposition of the sodium and potassium currents as well as the
reconstruction of the pace-maker property for the cardiac cells. The Noble model
adopted all the gating variables from Hodgkin-Huxley, and was also a degree-four
difierential system.
ThenextsigniflcantcardiacmembranemodelfollowingNoble?sistheMcAllister-
Noble-Tsien model published in 1975 [14]. By 1962, the theory of long-lasting
action potentials had developed much further than the experimental work and it
did not seem fruitful to explore further theoretical modiflcations due to the lack
of the quantitative information of the ionic currents as complete as that provided
by Hodgkin and Huxley in the case of squid nerve. The McAllister-Noble-Tsien
model, which was based on a wide range of experimental results for cardiac muscle
obtained using the voltage clamp technique, was not constructed until su?cient
experimental information was accumulated, about 10 years after the technique was
flrst successfully applied to cardiac membranes in 1964. The model formulated the
electrical activity of cardiac Purkinje flbers. In this model, the total ionic current
is broken down into nine discrete individual ionic  uxes. There are eleven gating
variables, among which nine are time-dependent and described by flrst-order dif-
ferential equations, and this makes the McAllister-Noble-Tsien model a degree-ten
ODE system.
The Beeler-Reuter model, published in 1977, is a pioneering efiort to describe
the cardiac ventricular active potential [15]. This work incorporated the majority
of the experimental data of ventricular myocardial action potentials by using the
6
voltage clamp method. The total ionic  ux is divided into only four individual
ionic currents in the Beeler-Reuter model: a time-independent outward potassium
current, a inward sodium current, which is primarily responsible for the rapid
upstroke of the action potential, a slow inward calcium current, and a outward
time-dependent potassium current. The main contribution of the Beeler-Reuter
model is that it includes the intracellular calcium ion concentration, which plays
a dominant role in the creation of the myocardial action potential plateau. The
model contains eight coupled, flrst order, non-linear difierential equations, which
are for the membrane potential, the intracellular calcium ion concentration, and
six gating variables for the various membrane conductances.
The experimental data used in the Beeler-Reuter model was subject to limita-
tions in available voltage clamp techniques and their application to multicellular
preparations of cardiac muscle [16]. These limitations were overcome with the
single-cell and single-channel recording techniques developed in the 1980s. The
Luo-Rudy I model [17], a system of degree-nine ODEs, was based on the data de-
rived from the new measurement techniques. The Luo-Rudy I model was published
in 1991, and described the electrophysiology of guinea pig ventricular cells. It is a
signiflcant update of the Beeler-Reuter mammalian ventricular model. The model
reformulated the fast inward sodium current and the outward potassium currents,
and investigated the phenomena dominated by these currents. The slow inward
current developed in the Beeler-Reuter model, which is to support the plateau of
the action potential, was retained.
The Luo-Rudy I model was further developed in 1994 to become the Luo-
Rudy II model [18]. The phase-II model reformulated the slow calcium current,
and incorporated a more thorough description of the processes which regulate
7
the intracellular calcium ion concentration and the movement of calcium ions.
This model is also based on experimental data of guinea pig ventricular cells from
single-cell and single-channel experiments. It provides a framework for future
development of models of the excitation-contraction coupling process in cardiac
cells. The model consists of eleven types of membrane currents and four types
of calcium currents that move to and from sarcoplasmic reticulum. There are six
gating variables and seven species of ion concentration formulated using difierential
equations, which make the Luo-Rudy II model a system of degree-fourteen ODEs.
Reviews of more cardiac cell ionic models are available in [5][12][19]. While the
experimental methods are becoming more sophisticated to obtain data for con-
structing more accurate ionic models that take into account more ion  ow mecha-
nisms, we adopt Beeler-Reuter?s ventricular cell model in our work for simulating
the action potential. The advantage lies in its less complex formulation, and hence
less simulation and hardware cost, while still keeping the description of the basics
of membrane ionic  ows between the intra- and extracellular media.
1.3 MotivationforVLSIImplementationofHeart
Model
Computer simulation of cardiac tissues is extremely computationally intensive.
Combining the equations for transmembrane ion  ows and binding processes pro-
duces a virtual cell as a complex nonlinear system of difierential equations. We can
take a relatively simple electrophysiology model, say, Beeler-Reuter as an example.
The model describes a system of eight coupled, flrst order, non-linear difierential
equations, which are for the membrane potential, the intracellular calcium ion
8
concentration, and six activation or inactivation parameters for the various mem-
brane conductances. Considering a sheet of 2-dimensional cardiac tissue of 1 cm2,
assuming each cell is 100 ?m? 100 ?m in size, there are totally 104 cells in this
tissue and hence there are 8?104 coupled difierential equations to represent the
virtual model. Furthermore, the demand of high spatial discretization and small
time-step sizes needed by the computational accuracy in turn requires the use of
huge memories and very fast processors. An approximate way to quantifying the
computational load is to estimate the number of  oating point operations that
would be necessary to compute one heartbeat [20]. Suppose a heartbeat takes
about 1 s, and the time step of the computation is 0.01 ms and the space step
is 100 ?m, then it would require some 1014  oating point operations to simulate
the cardiac electrical activity. This would take at least 30 hours of run-time on a
1GHz processor.
To handle the simulation request of huge memory and high computational
speed, many researchers have explored computing parallelism as the solver. In the
work presented in [21], the author, Pavarino, ran the simulation of the cardiac
reaction-difiusion system on two platforms. One platform employed an IBM SP
RS/6000 with Power 4 processors. The other consisted of an HP SuperDome 64000
with PA8700 processors. The work simulated a complete cardiac cycle (excitation-
recovery) in a slab of cardiac tissue of size 4?4?0.5 cm3 and mesh 400?400?50
using the Luo-Rudy I model as the membrane model and the Monodomain model
as the cardiac reaction-difiusion model. An adaptive algorithm was applied on the
time steps. To complete the simulation, the HP platform with 32 processors took
20 hours, and the IBM SP4 machine with 64 processors took about 2.5 hours.
However the performance does not scale up linearly as the number of parallel
9
processors increases. Pormann [22] presented the results of computational speed-
up vs. the number of parallel computing nodes. A cluster of 180 dual-CPU nodes
was used to do the investigation. Each node contained two 200 MHz Power3 CPUs
and 1 GB of memory. The results showed that to calculate 1.0 ms of electrical
activity of a 290?290?290 grid, with 0.001 ms in time-step, it takes 5386 seconds
(1.5 hours) using 1 CPU. The speed-up degraded with the increase of node number.
For a cluster of 32 CPUs, instead of achieving 32 times faster in speed, there was
only a speed-up of 21.97 due to the overhead of the communications between the
CPUs.
In this dissertation, we present an alternate but novel approach to simulate
cardiac electrophysiology using analog VLSI circuits. As stated in the previous
sections, the complexity of the heart leads to tremendous complexity in computer
simulations. The complexity arises from the nonlinearities caused by the extremely
highly coupled difierential equations. In order to solve the nonlinear system, time
is flnely discretized and the three dimensions are well meshed. Computation is
carried out for each time step and spatial grid. The introduced analog VLSI of the
heart model can be made on chips which serve as processing devices. Computation
with analog hardware has the advantage of obtaining fast solutions for time-based
difierential systems due to the difierential nature of capacitor components, i.e.
the derivative of a capacitor?s voltage with respect to time is proportional to its
current. By using analog VLSI hardware, the simulation of the heart model can
be performed with great speed-up. In the situations that people are attacked by
heart diseases, time is very critical to save lives. The VLSI implementation of the
simulation allows the fast analysis of fatal heart behavior and, thus, may allow for
obtaining real-time cure.
10
Compared to running simulation on digital computers, the simulation using
analog VLSI may require smaller hardware so that it is easy to carry it out. Digital
circuits operate based on the mathematics of Boolean logic. Their transistors work
as switches, and each wire represents a single bit with a value of either 0 or 1.
Analog circuits operate with continuous signals. Their primitives of computation
arise from the physics of the devices. One wire represents many bits, and hence
the amount of computation per a single transistor is usually higher than for digital
circuits. For example, an 8-bit current multiplier can be implemented with 8
transistors in analog circuits, whereas, digital implementation takes approximately
3000 transistors [23]. The area usage and power consumption demanded by an
analog system for a computation is also more e?cient than a digital system in
general, because more wiring and communication overhead is required in digital
circuits due to the greater number of transistors [23]. Therefore, using analog VLSI
to perform the simulations for a heart model can lead to lower power dissipation,
and less area and devices than using the general-purpose CPUs or other dedicated
digital systems.
In this dissertation, we introduce a methodology of simulating cardiac electro-
physiology using analog VLSI circuits. We consider one of the most used detailed
membrane models in the literature, the Beeler-Reuter model, introduced in Section
1.2. The Beeler-Reuter model mathematically reconstructs the action potential of
a mammalian ventricular cell. There are four ion currents and six gating variables
described in this model. More details are provided in the next chapter. It is worth
mentioning that though the introduced work targets the ventricular model devel-
oped by Beeler-Reuter, the proposed methodology of VLSI implementation can
be applied to any other cell model and can be used for implementing almost any
11
time-based difierential system.
1.4 Contributions
In this work, we investigate how to design an analog VLSI system to simulate
a cardiac electrical activity model which is usually done using high-performance
computers. We present the design methodology and the circuit components of
implementing the mathematical model using VLSI. The ventricular myocardium
model of Beeler-Reuter is realized with VLSI circuits to demonstrate the feasibility
of simulating complex numerical models using analog circuits.
The main contributions of this work are:
? A design of an analog VLSI circuit that simulates the Beeler-Reuter?s cardiac
electrophysiology model.
? An RC circuit architecture that calculates reaction-difiusion equations, and
the application of the RC circuit architecture to model the electrical propa-
gation process on the heart.
? A method to transform the original formulas to mathematical expressions
that are suitable to realize by VLSI circuits.
? Schemes for quantitatively scaling the equations to adapt to the computa-
tional circuits which are restricted by device parameters and working regions.
? An approach to change the initial conditions of difierential equation systems
to zero to save the elaborate hardware necessary to set the initial voltages of
integration capacitors.
12
? A set of transistor circuit blocks that are useful for implementing basic math-
ematical expressions, especially the equations in membrane ionic models.
? A methodology for implementing time-based difierential equations, numerical
models, and general reaction-difiusion systems using analog VLSI circuits
with special reference to difierential equations describing the heart.
1.5 Organization of Dissertation
The dissertation is organized in flve chapters.
Chapter 2 introduces the basics of cardiac electrophysiology related to this
work. It starts with presenting cellular physiology. A general description of cell
membrane is given, followed by the introduction of cardiac action potential and
its ionic mechanism. The formulation of the membrane ionic models is presented
next. Finally the ventricular myocardium model developed by Beeler-Reuter is
discussed in detail.
Chapter 3 flrst introduces the design  ow of implementing the mathematical
formulas using analog VLSI circuits. The transformation of mathematical expres-
sions to forms that are feasible for circuit implementation is discussed. Then, the
scaling of the equation parameters and the shifting of initial conditions are pre-
sented, followed by the introduction of a set of transistor circuit blocks used for
realizing basic mathematical functions and specialized for implementing the car-
diac model. An analog VLSI design of the Beeler-Reuter model is given next. The
circuits for realizing four ionic currents and the action potential are introduced
and the implementation accuracies of the circuits are discussed.
Chapter 4 provides the circuit simulation of the electrical propagation in cardiac
13
tissues. It starts with presenting the core-conductor model of the cardiac propa-
gation process, which is described by a reaction-difiusion equation. The method
of mapping a generalized reaction-difiusion system to an RC circuit is introduced
next. Then, the PSpice simulation of a transistor circuit for modeling a segment of
a 1-dimensional cardiac flber is discussed. The reentry phenomena and the spiral
wave simulated by the introduced RC network are presented flnally.
Chapter 5 concludes our work in this dissertation, and provides a discussion of
the future works at the end.
14
Chapter 2
Heart Electrophysiology and Its
Mathematical Models
2.1 Introduction
The contraction of the heart is controlled by the cardiac electrical conduction
system. This system generates electrical impulses, known as action potentials, and
conducts them rapidly throughout the heart, stimulating the heart to contract and
pump blood. The creation of the electrical impulses is tightly related to the ionic
mechanisms of the cardiac cells.
Each cell in our body is surrounded by a thin membrane [24](pp. 3). A cell
membrane, also known as a plasma membrane or plasmalemma, is a selectively
permeable bilayer, which mediates the difiusion of ions and molecules into and out
of the cell and establishes difierent ion concentrations between the intra and extra-
cellular space. The established ion concentration gradient also leads to an electrical
potential difierence, called the membrane potential, across the cell membrane due
to the charge of ions. The membrane potential is conventionally expressed as the
15
Intracellular
-+
Extracellular
Vm
Figure 2.1: Membrane potential Vm is the intracellular electrical potential with
respect to the extracellular media.
potential in the cell minus that on the exterior of the cell, and the extracellular
potential is taken as the ground, as illustrated in Figure 2.1. The membrane po-
tential is usually electrically negative when a cell is resting, i.e. when a cell does
not conduct any electrical signals. The membrane potential can be perturbed by
passing a positive current pulse from the exterior into the interior to charge a cell
and make it depolarized. Depolarization refers to the change in the membrane po-
tential to a more positive value [25][26](pp. 31). The larger the stimulating current
that passes into the cell, the larger is the change in the membrane potential. When
the applied current exceeds a certain strength, a threshold membrane potential is
reached, and this triggers the cell to flre an action potential, which is a much larger
response than the responses created from less stimulating currents.
The action potential created by a cardiac cell is described by models based on
the Hodgkin-Huxley formulation, and the Beeler-Reuter model is one of them, with
the latter being the one used in our work to simulate the cardiac action potential
with VLSI. The Beeler-Reuter model can be summarized with Figure 2.2, which
shows the dependence of the variables in the mathematical presentation. As we
will present in the last part of this chapter, the Beeler-Reuter model describes four
membrane ionic currents, six gating variables, and a varying calcium concentration,
16
Membranepotential
Vm
Sodiumcurrent CalciumcurrentTimeindependent potassium
current
Timedependent 
potassiumcurrent
IK1 Ix1 INa Is
mgate jgate hgate dgate fgate CalciumConcentration
[Ca]i
x1gate
Figure 2.2: Variable dependence in the Beeler-Reuter model.
and formulates a nonlinear system of degree-eight, flrst order difierential equations.
In the rest of the chapter, we will flrst introduce the cellular electrophysiology,
and explain membranes, ionic currents across membranes, and the cardiac action
potential. The mathematical models of ionic membranes are presented next. The
Hodgkin-Huxley model and the Beeler-Reuter model are described in detail flnally.
2.2 Cellular Electrophysiology
2.2.1 Cell Membrane
A cellular membrane is a lipid bilayer, which separates a cell from its extracellular
environment [24](pp. 3). It serves as a permeable barrier that permits the selective
passage of molecules across the membrane at a controlled rate. This selective
permeability allows a cell to maintain an interior composition difierent from the
extracellular media.
The direction of net ion movement across the membrane depends on both the
chemical force caused by the concentration difierence and the electrical force caused
17
by the electrical potential difierence. For a given type of ion A, the ion tends to
difiuse down the concentration gradient. Also, ion A tends to move from the
side with a higher electrical potential to the side with a lower electrical potential.
When the chemical and the electrical forces are equal in magnitude and opposite
in direction, no net movement of the ion occurs through the membrane, and the
ion is said to be in electrochemical equilibrium. The membrane potential that is
required to produce an electrical force that counterbalances the concentration force
is called the equilibrium potential or equivalently, the Nernst potential [26](pp. 22),
and can be calculated with the Nernst equation:
EA = RTzF ln
?[A]
e
[A]i
!
; (2.1)
where EA is the equilibrium potential of ion A, R is the ideal gas constant (8.3
J
mol?K), T is the absolute temperature in degrees K, z is the valence of the ion, F
is Faraday?s number, i.e., the charge of 1 mol of electrons (96,500 C/mol), [A]e is
the extracellular concentration of ion A, and [A]i is the intracellular concentration
of the ion. The polarity of the equilibrium potential is related to the direction of
a chemical force. A positive equilibrium potential indicates that the force caused
by the concentration gradient tends to make the ion  owing into the cell, and a
negative sign indicates the concentration force drives the ionic current out of the
cell. The magnitude of the equilibrium potential represents the strength of the
chemical force.
Ions in a cell are usually impacted by the same electrical potential, i.e., the
membrane potential, whereas difierent types of ions have various chemical poten-
tials due to difierent ion concentrations. Hence in most tissue, ions are not in
electrochemical equilibrium. Table 2.1 shows the typical ion concentrations of a
cardiac muscle cell, and the corresponding equilibrium potentials [24](pp. 3). For
18
Table 2.1: Typical ion concentrations and equilibrium potentials in cardiac cell
[24].
Intracellular Extracellular Equilibrium
Concentration (mM) Concentration (mM) Potential (mV)
Sodium(Na+) 10 145 +71
Potassium(K+) 140 4 -95
Calcium(Ca2+) 0.0001 1.5 +128
Na+ and Ca2+, the difiusion caused by the concentration gradient is inward, since
their intracellular concentrations are much less than the extracellular media. In
opposite, K+ tends to move outward in the efiect of its concentration gradient. A
ventricular cardiac cell normally has an action potential that ranges from approxi-
mately -90 to 40 mV [27]. When the action potential of a ventricle is negative, the
electrical force is in the same direction as the concentration force of Na+ and Ca2+,
and this causes a net inward movement of Na+ and Ca2+. The chemical force on
K+, though in the opposite direction of the electrical force, has a larger magnitude,
known by comparing the magnitude of K+ equilibrium potential, 95 mV, with the
maximal absolute value of the membrane potential, 90 mV. Consequently, K+
shows a tendency of inward movement. When the action potential is positive, the
direction of the ionic currents can be deduced in a similar way. As a result, the
Na+ and Ca2+  ows are always inward and the K+  ow is always outward during
the period of a ventricular action potential.
19
2.2.2 Ionic Mechanism of Action Potential
The permeability of the membrane is mediated through ion channels that act
as pathways for the passage of charged molecules. Ion channels are specialized
proteins embedded in the membrane that allow only desired ions to cross the
membrane [28][29]. Most ion channels are gated and are capable of  ipping between
a conducting (open) and non-conducting (closed) state. The transition of the state
is a random event [25],and hence it is not possible to predict whether a given
channel is open or closed. However, the laws of probability allow people to make
certain predictions of the average behavior of a channel. The probability of being
open can be estimated by calculating the fraction of time a channel spends in the
open state during a period of time. The open probability can be regulated by
chemical or electrical signals, temperature, or mechanical force [30].
The ionic currents that occur through protein ion channels make possible the
important electrical activities based on the voltage changes of the membrane poten-
tial. The rapid movement of ions via ion channels sometimes can create electrical
currents that are large enough to produce rapid changes in the membrane poten-
tial. The rapidly changed membrane potentials are termed action potentials, and
can be created by a stimulating current that exceeds a certain strength. An inward
current pulse applied to a cell membrane causes the membrane potential to depo-
larize, and the size of the depolarization depends on the magnitude of the stimulus.
When the current pulse reaches above a certain value to make the membrane po-
tential exceed a threshold, the cell flres an action potential, a much larger response
with a difierent shape compared to the potential changes resulted from a smaller
stimulus. This is illustrated in Figure 2.3. Figure 2.3(a) shows four stimulating
pulses with difierent amplitudes. Figure 2.3(b) shows four time courses of chang-
20
(a) (b)
Time [ms] Time [ms]
Istim
 [
a0 A
/cm
2 ]
Vm
[mV
] 4
32
1
4
3
2
1 threshold
Figure 2.3: Responses of membrane potential to various stimulus. (a)Four stim-
ulating currents with difierent strength. (b)Corresponding membrane potentials
in response of the stimulus. An action potential is created when the membrane
potential is charged by the 4th stimulus to pass a threshold.
ing membrane potentials, each of which represents the response to the stimulus
of the same number in Figure 2.3(a). Among the stimulus, the number 4 current
impulse is the only one that is able to charge the membrane potential to reach the
threshold at around 60 mV and make an action potential.
A cardiac action potential pulse is the electrical activity of a individual cardiac
muscle cell. The cardiac action potential voltage difiers in difierent portions in the
heart [31](chapter 6). This difierentiation of the action potentials allows the dif-
ference in the electrical characteristics and functions of the varied types of cardiac
cells. The numerical simulation models of the myocardial action potential have
been reported for a few difierent locations on the heart: sinoatrial node, Purkinje
flber, atrial myocardium and ventricular myocardium [14]-[18] [32]-[35].
Figure 2.4 shows a waveform of the action potential in a ventricular cell. The
21
Na+- -- - Na+ K++ + K+ Ca
2+ Ca2+ - -
- -+ + 
+ + 
+ + K
+ K+ K
+ K+ - -- -K+ K+
Figure 2.4: Phases of ventricular action potential and active ionic currents in each
phase. (0)upstroke, (1)early re-polarization, (2)plateau, (3)late re-polarization,
(4)rest.
ventricular action potential has 5 phases [36][37][38], established due to the fast
changes in membrane permeability to certain ions, mainly to sodium, potassium,
and calcium, resulting from the opening and closing of the ion channels. Figure
2.4 summarizes the ionic channel currents that occur during each phase of the
ventricular action potential. The net movement direction of the ions depends on
the polarity of the intercellular potential with respect to the exterior, denoted with
\-" and \+" inside the cells in the flgure, as well as the concentration difierences
between the interior and exterior of the cell, indicated with difierent sizes of the
fonts for the ion names. The arrows represent the directions of the ion  ows.
Phase 0 is the rapid depolarization phase. It arises from the sudden opening
of the sodium channels and the subsequent rapid in ux of the sodium current.
The rapid opening of the activation gate is responsible for the large and abrupt
increase in the sodium conductance in the membrane. The activation gate refers
22
to the gates that tend to open when the membrane potential becomes less nega-
tive. The sodium continues entering the cell when the channels are open and this
eventually reverse the polarity of the membrane potential. The inward sodium
 ux flnally stops and phase 0 is terminated due to the closure of the other type of
sodium gates, called inactivation gates, which are deflned as the gates that tend
to close when the membrane potential becomes less negative. Phase 0 is followed
by a slight repolarization period, called phase 1. The repolarization results from
the out ow of potassium. After the repolarization, a plateau stage follows, known
as phase 2. At this stage, the calcium channels start to open and this allows an
in ux of positive charge, Ca2+, under the in uence of the inward electrochemical
potential for calcium. The inward  ow of calcium and the outward  ow of potas-
sium currents roughly cancel each other. The membrane potential then reaches the
steady plateau. After the plateau, the calcium channels close, and the potassium
current dominates and repolarizes the cell back to the resting potential, and this
makes up phase 3. In phase 4, the membrane potential stays in a constant level
of approximately -85 mV, a little more positive than the equilibrium potential
of potassium (refer to Table 2.1). Hence, the chemical force that favors the out
 ow of potassium exceeds very slightly the electrical force that favors the in ux
of potassium. Consequently there is a tiny outward potassium current during the
resting state of phase 4.
The ionic concentrations, which largely in uence the membrane potential, are
maintained by the active transport [5](pp. 161) [26](pp. 14) [39](chapter 15). The
active transport is also protein-mediated. In addition, this process allows ions
to move through membranes against an electrochemical potential, and hence it
requires energy. The concentration gradient of Na+ and K+ are maintained by a
23
Na+=K+ pump, an active transport that moves three Na+ ions out of and two
K+ ions into the cell with the consumption of an APT(Adenosine Triphosphate
- a type of chemical energy within cells) ([39](chapter 15). The Ca+ pump helps
maintain the concentration gradient of Ca+. There are two types of Ca+ pump
involved. One type transports Ca+ out of the cell. The other type transports
intercellular Ca+ into the sarcoplasmic reticulum(SR), an intercellular structure
that stores Ca+ ions. The release of Ca+ ions stored in a SR are triggered by
raised intracellular calcium resulting from transmembrane calcium in ux [40].
2.3 Mathematical Description of Cardiac Elec-
trophysiology
The development of experimental technologies necessitates the detailed under-
standing of the electrical behavior of a cellular membrane. The experiments in-
clude the measurement of the voltages across membranes in difierent spatial posi-
tions, current  ows, ion concentrations and opening states of single ion channels
[5] (pp.157). The commonly used measurement method is the voltage clamp tech-
nique, in which two electrodes are inserted into a cell to \clamp" voltage to a flxed
value and the current needed to maintain a desired voltage is measured. The quan-
tities obtained by the experiments are partly used to create mathematical models
of difierent levels of abstraction. The models allow the numerical simulation of the
electrophysiological behavior of cells, and assist the reconstruction of the measured
data and the further discovery of unknown phenomena.
In this section, We flrst introduce the voltage clam technique, and present the
mathematical formulation of the membrane models. The Hodgkin-Huxley model
24
Command Voltage +
-
Voltage clamp amplifier I
Axon
+
-
V-controlledcurrent
I
Vc+-
+Vm
Vext -
Figure 2.5: Measurement of ionic current using voltage clamp technique.
for squid axon and the ventricular myocardium model developed by Beeler-Reuter
are described next.
2.3.1 Voltage Clamp
The voltage clamp allows the measurement of ionic currents under the in uence of
the membrane potential [41]. As shown in Figure 2.5, the voltage clamp technique
involves two electrodes that are placed inside a cell. One electrode, used to measure
the intercellular potential Vm, is connected to an amplifler. The other side of the
amplifler is connected to the exterior of the cell to take the extracellular potential
Vext. The measured potential feeds into a voltage clamp amplifler to compare with
a voltage to be maintained by the membrane potential, called a command voltage
Vc. The output of the voltage clamp amplifler, i.e. the amplifled value of the
difierence between Vm and Vc, controls a feedback current  owing into the cell to
make these two potentials the same. The feedback current I, which is injected into
the membrane through the other electrode, is then the mirror image of the current
generated by the cell membrane at the command voltage, that is, when Vm = Vc.
In order to measure currents of individual ions, other experimental procedures can
25
be also involved, such as ion substitution, channel blockers, and speciflc clamp
protocols [41][42].
A wealth of new knowledge concerning ion channels resulted from the invention
of the patch clamp method, a reflned version of the voltage clamp. In this tech-
nique, a glass pipette with a tip diameter of about 1 or 2 ?m [2](pp. 42) is used
to form an exterior contact with a tiny area of a membrane. The contact is made
very tight, by applying a small suction, such that all the ions that  ow through the
opening channels in the sealed patch of membrane  ow into the pipette. If there is
only one ion channel in the patch, the resulting electrical current, measured with
an ultrasensitive electronic amplifler, is due to the opening and closing of the sin-
gle channel, and the measured value transits randomly between zero and a certain
non-zero value. The patch clamp method allows the study of a single ion channel
and gains further insight into the electrophysiology of membranes.
2.3.2 Ionic Membrane Models
Cell Membrane as Capacitor Circuit
The electrical behavior of a cell membrane can be approximated by a circuit with
a capacitor connected in parallel to current sources as illustrated in Figure 2.6(a)
[7]. From the equivalent circuit, the membrane current can be expressed mathe-
matically as the following equation:
Iext = CmdVmdt +Iion; (2.2)
Iion =
NX
i=1
Ii; (2.3)
where Iext is the membrane current (inwards current positive); Iion is the total
ionic current (outwards current positive); Vm is the membrane potential; Cm is the
26
?Cm
I1 I2 IN
-
+
Vm
Iext
IA EA
RA -
+
Vm
(a) (b)
Extracellular
Intracellular
+
-
Figure 2.6: Equivalent circuits for Hodgkin-Huxley type of membrane models.
(a)Membrane is approximated by a circuit with a capacitor and ionic currents.
(b)Each ionic current represents a type of ionic  ow and can be modeled as a
nonlinear resistor connected with a Nernst potential. RA = 1=(gAy).
membrane capacitance; and t is time. Iion characterizes the total transmembrane
currents caused by difierent ion  ows. The component Iext represents the external
in uence applied on the cell. If all the ionic currents  ow into the membrane
capacitor, i.e. Iext = 0, equation (2.2) and (2.3) simulate the electrical activity of
an isolated cell. The permeability of the membrane for each type of ion varies along
the changes of the membrane potential, and reversely, the varying ionic currents
are attributed to the changing of the membrane potential. Next, we will discuss
the mathematical model for transmembrane ionic currents.
Mathematical Models of Ionic Channels
Almost all the ionic models of biological cells are inspired by the Hodgkin-Huxley
equations, which formulate in detail the ionic  ows between the intracellular and
extracellular media. As suggested by Hodgkin and Huxley [7][42], the ionic perme-
27
ability of the membrane can be expressed in terms of ionic conductances, whose val-
ues are determined by the number of opening channels, and thus can be modulated
by the membrane potentials. The  ow of ion is dependent on the conductances of
the ion channels, and is also in uenced by the ion concentration gradient and the
difierence of the potential across the membrane. Therefore, an ionic current can
be expressed in terms of Ohm?s law:
IA = gA ?y ?(Vm ?EA); (2.4)
where IA is the current associated to the ion A (mA=cm2); gA is the maximal
conductance of the corresponding ion channel (mS=cm2); y is a dimensionless
variable which represents the proportion of the ion channel in an open state (0 <
y < 1); Vm is the membrane potential (mV); and EA is the equilibrium potential of
ion A (mV), expressed by equation (2.1). Like in the Ohm?s law, an ionic current
is expressed as a voltage multiplied by a conductance. This can be illustrated by
Figure 2.6(b), in which RA is the equivalent conductance of the ion channel and
is described as the inverse of the product of gA and y. The voltage across RA
is Vm ? EA. The sign of the calculated result from equation (2.4) indicates the
direction of the ionic current. A positive sign represents an outward  ow, and a
negative sign represents an inward  ow.
Ion channels can be governed by more than one gate [43]. A channel is only
fully open when all the gating variables reach the maximum values. Therefore, y
can be modeled as a product of serial gates:
y = y1y2:::yM; (2.5)
where yi (i = 1 ? M) is the open proportion of the ith gate. The behavior of the
gates can present difierent properties, with some having more opening probability
28
at higher Vm, called activation gates, and the others having more opening proba-
bility at lower Vm, called inactivation gates. A current  ow is totally inhibited if
any one of the gates is fully closed.
The gating variable is controlled by a voltage-dependent gating mechanism for
most channels. The open or closed state of the gate is determined by the membrane
potential, but most of the gates do not respond instantaneously to the voltage. The
time-dependence of the gating variables is conveniently modeled by a flrst order
difierential equation, as shown in the following, and its rate of change depends on
two coe?cients, namely a gate opening rate and a gate closing rate:
dyi
dt = fii(1?yi)?fliyi; (2.6)
where fii is the opening rate coe?cient of the ith gate and fli is the closing rate
coe?cient. For the voltage-dependent gating channels, these rates are controlled
by the membrane potential. In order to obtain the functions connecting fii and fli
with the membrane potential, experimental measurements have been carried out
and data processing is performed, as presented in the following.
When a cell is in the resting state, yi has a resting value given by:
yi0 = fii0(fi
i0 +fli0)
: (2.7)
The \0" in the variable names in the above equation denotes that these are the
values associated with the resting state. When the membrane potential is changed
suddenly, say, to V 0m, fii and fli instantly take up the values related to the new
potential. Therefore, the changes of yi along the time can be expressed by the
analytical solution of equation (2.6), given in the following:
yi = yi1 ?(yi1 ?yi0)e(? t?i); (2.8)
29
where the time constant is:
?i = 1=(fi0i +fl0i); (2.9)
and the flnal value is:
yi1 = fi0i=(fi0i +fl0i): (2.10)
fi0i and fl0i are the opening and closing rates associated with the new membrane
potential V 0m.
The time course of yi can be measured by recording the membrane currents
using the voltage clamp technique [42], and setting the command voltage to a step
function jumping from the resting membrane potential to V 0m. The values of ?i
and yi1 can be obtained by drawing a smooth curve from the equation (2.8) to flt
the experimental data points. fi0i and fl0i can then be calculated from the following
relations which are derived from equation (2.9) and (2.10):
fi0i = yi1?
i
; (2.11)
fl0i = (1?yi1)?
i
: (2.12)
By changing the value of V 0m, we can flnally collect enough fi0i and fl0i data points at
difierent V 0m to obtain the formulas which express fii and fli in terms of Vm. As we
will see later in the Hodgkin-Huxley and the Beeter-Reuler models, the formulas
are difierent for difierent types of gates.
The derivative of fii and fli using the above experimental method relies on the
assumption that V 0m stays flxed when a gating variable is changed from its original
value to the value related to V 0m. As in the time course of an action potential, the
membrane potential is always changing, violating this assumption in which case
the expressions of fii and fli for the dynamic process may not be exactly consistent
with the experimental results.
30
Equation Set of Ionic Membrane Model
Putting together all the information introduced above, the ionic membrane model
is summarized as the following. The membrane potential is described as a capac-
itor charged by the ionic currents and an external current (a rewritten version of
equation (2.2)):
dVm
dt = ?
1
Cm
NX
i=1
Ii + IextC
m
; (2.13)
and each of Ii is formulated by (refer to equation (2.1) and (2.4))):
Ii = gi
MiY
j=1
yij(Vm ?Ei) (i = 1 ? N); (2.14)
Ei = RTzF ln
?[C
i]e
[Ci]i
!
: (2.15)
The above equations assume there are N types of ionic currents, I1 ? IN, and
there are Mi number of gates, yi1 ? yiMi controlling the ion channel of ion type i.
When the variation of ionic concentration is considered, additional equations are
required to describe the extracellular and intracellular concentration of ion i, [Ci]e
and [Ci]i [18][44][45]. The usage of the Nernst equation for the cases in which the
concentrations are changing is an approximation. The Nernst equation is satisfled
only for ions at electrochemical equilibrium, i.e., no net movement of the ion occurs
[26](pp. 22).
The opening portion of individual gates can by expressed by (refer to equation
(2.6)):
dyij
dt = fiij(1?yij)?flij ?yij (j = 1 ? Mi); (2.16)
fiij = Fij(Vm); (2.17)
flij = Gij(Vm); (2.18)
where Fij and Gij are functions of Vm, obtained from experiment as introduced
31
Cm
INa IK Il
-
+
Vm
Iext
RNa
Extracellular
Intracellular
ENa
-
EK
+
-
El
+
-
RK Rl
Figure 2.7: Equivalent circuits for the squid axon model of Hodgkin & Huxley.
RNa = 1=(gNam3h); RK = 1=(gKn4); and Rl = 1=gl.
previously.
The above set of equations describe the basic electrophysiological behavior of
a single cell.
Hodgkin-Huxley Model
The Hodgkin-Huxley model describes the dynamic electrophysiology of a giant
squid axon membrane from measurements of its electrical behavior using voltage
clamp [7]. The membrane is modeled as a circuit consisting of a capacitor, resistors,
and voltage sources, as illustrated in Figure 2.7. The membrane potential is still
deflned with a flrst-order difierential equation (refer to equation (2.13)):
dVm
dt = ?
1
CmIion +
Iext
Cm; (2.19)
with all the notations following the same meanings of those in equation (2.2). The
model splits the ionic current into three terms:
Iion = INa +IK +Il; (2.20)
32
where INa is the voltage- and time-dependent sodium current, IK is the voltage-
and time-dependent potassium current, and Il is a voltage-dependent leakage cur-
rent made up by chloride and other ions. The currents are determined by the
conductances gNa, gK, and gl, respectively, as well as the difierence between the
membrane potential and their equilibrium potentials:
INa = gNa(Vm ?ENa) = gNam3h(Vm ?ENa); (2.21)
IK = gK(Vm ?EK) = gKn4(Vm ?EK); (2.22)
Il = gl(Vm ?El) = gl(Vm ?El); (2.23)
where gNa, gK, and gl are conductance constants, m, h, and n are gating variables,
and ENa, EK, and El are the Nernst potentials.
gNa and gK represent the maximal conductances for sodium and potassium
ions. The conductance for the leakage current Il is assumed to be a constant gl.
The values of these conductances are given by:
gNa = 120mS=cm2; gK = 36mS=cm2; gl = 0:3mS=cm2: (2.24)
The ionic concentrations are supposed to be invariant in the model, and this leads
to non-varying equilibrium potentials:
ENa = 55mV; EK = ?72mV; El = ?49:387mV: (2.25)
The gating variables m, h, and n are controlled by their opening and closing
rate coe?cients, and described by the following equations:
dm
dt = fim(1?m)?flm ?m; (2.26)
dh
dt = fih(1?h)?flh ?h; (2.27)
dn
dt = fin(1?n)?fln ?n; (2.28)
33
Table 2.2: Units used in the Hodgkin-Huxley equations
Description Unit Description Unit
Currents ?A=cm2 Potentials mV
Conductance mS=cm2 Time t ms
where the voltage-dependent rate coe?cients, extracted by curve fltting the exper-
imental data [7], are given by:
fim = 0:1(Vm +35)(1?e?Vm+35
10 )
; (2.29)
flm = 4e?Vm+6018 ; (2.30)
fih = 0:07e?Vm+6020 ; (2.31)
flh = 1e?Vm+30
10 +1
; (2.32)
fin = 0:01(Vm +50)(1?e?Vm+50
10 )
; (2.33)
fln = 0:125e?Vm+6080 : (2.34)
The units used in the Hodgkin-Huxley equations are listed in Table 2.2.
The Hodgkin-Huxley formulation assumes that the conductances of ion chan-
nels are continuous and deterministic. The assumption arises from the limitation
of the experimental method, voltage clamp, which is at the cellular level and only
capable of measuring the macroscopic ionic currents resulting from a population
of channels. The model can not be used to describe the discrete ion channels
and their random opening and closing behaviors, which are later formulated by
Markov processes after the patch clamp technique was invented and studies of
single channels became possible[46].
The Hodgkin and Huxley model was the flrst to describe the way in which
individual ionic currents vary with membrane potential and time. Its structure
34
forms a basis for almost all models of excitable membrane behavior.
2.3.3 Beeler-Reuter Model
Next, we focus on presenting the mathematical model used in our VLSI design
work, namely the Beeler-Reuter model. The employment of Beeler-Reuter?s model
in this work for simulating the action potential relies on its less complex formu-
lation, while still keeping the description of the basics of membrane ionic  ows
between the intra- and extracellular media.
The Beeler-Reuter model was constructed to describe the electrophysiology of
a mammalian ventricular myocardium in 1977 [15], and incorporated the majority
of the experimental evidence achieved at that time by using the voltage-clamp
techniques. The model represented a numerical simulation of the ventricular action
potential, and described a system containing four types of ionic current and six
gating variables. The ion  ux of the model includes: a time-independent outward
potassium current, IK1; a time-activated potassium outward current Ix1; a fast
voltage- and time-dependent inward current carried primarily by sodium, INa; and
a slow voltage- and time-dependent inward current, Is, carried mainly by calcium
ions. The formulation of the time and voltage dependence of the gating variables
in the Beeler-Reuter model follows the Hodgkin-Huxley equations.
The Beeler-Reuter model reproduces a typical ventricular myocardial action
potential as shown in Figure 2.8. The threshold of excitation is -60 mV (refer to
Figure 2.3(b)) [15]. During the rapid depolarization phase, the upstroke velocity
is 115 V=sec. The peak of the upstroke is about 30 mV. In the plateau phase
the potential reaches maximum 17 mV, and the positive potential lasts for 153
msec. The maximum rate of repolarization is about 11 V=sec. The duration of
35
+40 
-40 
-80 
Time [ms]
0 
Vm 
 [mV
]
Threshold
Figure 2.8: Simulated ventricular action potential using the Beeler-Reuter model.
the action potential, measured at the point where re-polarization is 90% complete,
is 285 msec. The resting potential is -84 mV.
Ionic currents in Beeler-Reuter Model
The formulation of the potassium current IK1 and Ix1 is based on the available ex-
perimental results in ventricular myocardium and follows the approach presented
by McAllister et al. in [14], which mathematically described the electrophysi-
ology of a Purkinje flber. The experimental evidence indicates the presence of
two outward currents: a background current and a single time-activated current.
The background current, modeled as IK1, is determined only by the membrane
potential, and thus can be simply described by a function in terms of Vm. The
time-activated outward current, modeled as Ix1, is observed to show a non-linear,
rectifying, current-voltage relationship when fully activated. In order to charac-
terize this current, the Beeter-Reuler model adopted the formulation proposed by
McAllister et al., which is analogous to the linear case, as shown in the following:
Ix1 = Ix1 ?x1; (2.35)
36
where Ix1 is the maximum current Ix1 can reach, and x1 is a gating variable, formu-
lated by a flrst-order ODE like equation (2.6). Both IK1 and Ix1 are outward due
to the dominating outward chemical force caused by the much larger intracellular
potassium concentration compared to the extracelllar media.
The inward sodium current INa can be divided into two components [5](pp.
174), as described as follows:
INa = INaC +INaV : (2.36)
INaC represents the time independent component in the sodium current, and is to
reproduce the measured steady sodium leakage current. It is expressed as the fol-
lowing equation, which is very similar to equation (2.4), except that the membrane
conductance is not controlled by any gating variables in this case:
INaC = gNaC(Vm ?ENa): (2.37)
The second component INaV is time- and voltage-dependent, and is responsible for
the fast upstroke of the action voltage. This sodium current is governed by three
types of gates, described as follows:
INaV = gNaV m3hj(Vm ?ENa); (2.38)
where gNaV are the maximal conductance for the time-varying sodium current,
and m, h, and j are the three gating variables. This formula adopts the sodium
gating variable m determined for squid axon by Hodgkin & Huxley, to simulate the
rapid depolarization phase (refer to equation (2.21)). The re-polarization process,
in which the sodium current diminishes due to closing of the inactivation gates,
cannot be simulated with a single simple parameter for the ventricular action
potential. Therefore, two parameters, h and j, are introduced to describe the
37
process. All these three gating variables are formulated using equation (2.8), with
difierent fii and fli. The Beeler-Reuter model assumes that the ionic concentration
of sodium keeps unchanged during the action potential, and hence the sodium
equilibrium potential ENa is a constant.
The slow inward calcium current in the Beeler-Reuter model is expressed by:
Is = gsdf(Vm ?Es): (2.39)
When all the gates are fully open, Is is determined by the linear current-voltage
relation: gs(Vm ? Es) [47][48], in which, gs deflnes the maximal conductance of
calcium, and Es is the equilibrium potential for calcium. In reality, only a fraction
of the ionic channels are open, and this is described by using the gating variables
d and f. Their time derivatives are functions of the membrane potential. Here we
adopt the original notations from [15], and d and f correspond to y1 and y2 in equa-
tion (2.5). Particular attention was given by Beeler-Reuter to model the calcium
equilibrium voltage, Es, by considering the change of calcium concentration inside
a cell. The formulation of the calcium concentration is a flrst approximation to
the experimental results obtained by Bassingthwaighte & Reuter in [49], in which
the interrelationship of clamped membrane potential, observed calcium current,
and the equilibrium potential were studied. The in ux of calcium is modeled by
treating the calcium current as though it  ows into a small distribution volume
within the cell, from which the calcium concentration is reduced exponentially by
an uptake mechanism. The equilibrium potential Es is then calculated from the
calcium concentration with the Nernst equation.
38
Equations and Values that Deflne Beeler-Reuter Model
The time derivative of the membrane potential is deflned as before:
dVm
dt = ?
1
CmIion +
Iext
Cm; (2.40)
but now,
Iion = IK1 +Ix1 +INa +Is; (2.41)
where IK1, Ix1, INa, and Is sum up to be the total ionic current Iion, and Iext is
an external current that can be used as a stimulus to trigger the action potential.
The membrane capacity Cm is set at 1 ?F=cm2 in the model, which is a generally
accepted value for the capacity of biological membranes [15]. The unit of the
membrane potential, as well as other potentials used in the model, is set to mV.
The time t has a unit of ms.
Some of the equations of the ionic currents have been given previously. Here
we put these current equations together and provide the values of the constants
and the function expressions of some variables (the gating variables and Es are
given in equations (2.46)?(2.60)).
IK1 = 0:35
" 4(e0:04(V
m+85) ?1)
e0:08(Vm+53) +e0:04(Vm+53) +
0:2(Vm +23)
1?e?0:04(Vm+23)
#
; (2.42)
Ix1 = Ix1x1 = 0:8(e
0:04(Vm+77) ?1)
e0:04(Vm+35) x1; (2.43)
INa = gNaC(Vm ?ENa)+gNaV m3hj(Vm ?ENa)
= (0:003+4m3hj)(Vm ?50); (2.44)
Is = gsdf(Vm ?Es) = 0:09df(Vm ?Es): (2.45)
Since the experimental results were measured on a space-clamped patch of mem-
brane of about one square centimeter, all ionic current in the above equations are
actually current densities, and have a unit ?A=cm2.
39
The six gating variables x1, m, h, j, d, and f are deflned by a set of equations
similar to equation (2.6):
dz
dt = fiz(1?z)?flz ?z; (2.46)
where z is x1, m, h, j, d, or f, and their opening and closing rate coe?cients are
described by the following equations:
fix1 = 0:0005e
0:083(Vm+50)
e0:057(Vm+50) +1 ; (2.47)
flx1 = 0:0013e
?0:06(Vm+20)
e?0:04(Vm+20) +1 ; (2.48)
fim = ?(Vm +47)e?0:1(V
m+47) ?1
; (2.49)
flm = 40e?0:056(Vm+72); (2.50)
fih = 0:126e?0:25(Vm+77); (2.51)
flh = 1:7e?0:082(V
m+22:5) +1
; (2.52)
fij = 0:055e
?0:25(Vm+78)
e?0:2(Vm+78) +1 ; (2.53)
flj = 0:3e?0:1(V
m+32) +1
; (2.54)
fid = 0:095e
?0:01(Vm?5)
e?0:072(Vm?5) +1; (2.55)
fld = 0:07e
?0:017(Vm+44)
e0:05(Vm+44) +1 ; (2.56)
fif = 0:012e
?0:008(Vm+28)
e0:15(Vm+28) +1 ; (2.57)
flf = 0:0065e
?0:02(Vm+30)
e?0:2(Vm+30) +1 : (2.58)
The equations for calcium intercellular concentration [Ca]i and the equilibrium
potential Es are given here:
Es = ?82:3?13:0287ln[Ca]i; (2.59)
d[Ca]i
dt = ?10
?7Is +0:07(10?7 ?[Ca]i): (2.60)
40
Table 2.3: Initial conditions for the Beeler-Reuter Model
Vm x1 m h j d f [Ca]i
(mV) (mol=l)
?84 0:0088 0:01979 0:9464 0:937 0:003763 1 10?7
The unit of the calcium concentration is mol=l.
When a cardiac cell is resting, the difierential system is at an equilibrium state,
and all the variables remain constant until the system is disturbed by a stimulus.
The values of the derivative variables at a resting membrane are taken as the initial
conditions for the difierential equations. The initial values are obtained by setting
all the time derivative terms to be equal to zero, that is the variables are not
changing along the time. The initial conditions are listed in Table 2.3.
All the above equations are incorporated into the Beeler-Reuter model, which
numerically allows simulation of the action potential of mammalian ventricles. The
variable dependence in the equation set is illustrated in Figure 2.2. The model
formulates four ionic currents, IK1, Ix1, INa, and Is, and describes a degree-eight
system of flrst order difierential equations, which are for the membrane potential,
Vm, six gating variables, x1, m, h, j, d, and f, and a varying calcium concentration,
[Ca]i.
41
Chapter 3
VLSI Design of the Beeler-Reuter
Model
3.1 Overview
To achieve solutions to hard mathematical problems, numerical analysis is usu-
ally carried out, in which continuous problems are discretized and functions are
represented by a flnite amount of data. For a given problem, the numerical anal-
ysis includes flnding an iterative method that leads successive approximations to
converge to solutions, discretization of continuous domains with a flnite number of
points, and the study of errors and numerical stability [50]. Software programming
is usually involved in implementing the algorithms of numerical computation. In
this chapter, we will present a very difierent technique for obtaining numerical
solutions of mathematical problems, that is using analog VLSI circuits.
The primitives of the analog VLSI implemented computation arise from the
physics of the devices, and their functioning greatly depends on the I ?V char-
acteristics of the devices used, especially the transistors. Computing with VLSI
42
EquationTransformation IdealComponent
Circuit
ParameterScaling CircuitFunctional
Units
Initial ValueShifting VLSIImplementation
OriginalModel
Circuit DesignModel Reformulation
Reformulated Equations
Section 3.2.2
Section 3.2.1
Section 3.2.3
Section 3.2
Section 3.4
Figure 3.1: Design  ow of VLSI Realization.
takes advantage of the integration properties of capacitors, and thus the calcula-
tion of complex difierential equations can be greatly sped up. The methodology
of the VLSI design process is illustrated in Figure 3.1. The design entry starts
from the original mathematical model, i.e. the Beeler-Reuter model in our case,
and splits the design  ow into two paths. One is, shown on the left side of the
flgure, the steps of reformulating the mathematical descriptions, and the other,
shown on the right side, is more directly related to the circuit designs, in which
the shaded oval is our flnal goal. Equation transformation (the left top oval) is the
process that changes equations to the formulas that are feasible for the realization
using circuits. The transformed equations next go through parameter scaling, a
step that scales the values of variables and parameters into proper ranges for the
circuit implementation. The scaling process needs the overall consideration of all
equations, and usually requires the guidance of the numerical information, which
can be achieved from the simulation of ideal component circuits (the right top
oval), and demands the acknowledgment of the restrictions of working ranges of
43
the functional circuits (the right middle oval), which is reversely decided by the
scaling step for the requirement of the function types to be realized. The scaled
equations then undergo initial value shifting (the left bottom oval), in which ini-
tial conditions of the difierential equations are altered to all zeros for alleviating
their complexity in the circuits. The reformulated equations are put into the flnal
implementing step, which takes the functional circuit units as building blocks and
the ideal component circuit as a rough blueprint.
In the rest of the chapter, we will flrst introduce the equation reformulation
techniques, namely, equation transformation, parameter scaling, and initial value
shifting, and then provide the resultant equations of the Beeler-Reuter model after
applying the reformulation process. The functional circuits that realize a set of
mathematical expressions are described next. The ideal component circuits and
the VLSI design of the Beeler-Reuter model are presented flnally.
3.2 Model Reformulation
3.2.1 Transformation of Mathematical Description
It is not surprising that circuits are capable of doing mathematical calculations.
One simple example is we can perform an addition a+b using a circuit in which,
two current sources with values a and b are wired together with a third wire. The
result of a+b is then the current in the third wire. Beneflted from Kirchhofi?s law,
linear operations, such as addition, subtraction, magnifying or diminishing with a
constant factor, can be easily realized by using currents as signal representatives.
Non-linear operations like multiplication and division are also possible using cir-
cuits. CMOS current multipliers have been published in [51]?[54], and multipliers
44
taking voltage as signals are also reported in [55]?[57], as well as current/voltage
dividers [58][59]. Other operations such as exponentiation are also realizable with
circuits [60][61]. However, the abundance of mathematical function types is far be-
yond the availability of the implementation circuits, and thus most mathematical
operations do not have corresponding circuit realization. In addition, all sorts of
operations can be combined and nested, which make the circuit implementation
very di?cult or even impossible. The complexities of employing circuits in func-
tion computation prompts the need for transforming the equations to other forms
that are feasible to be mapped into circuits.
Fortunately, for most equations in the Beeler-Reuter model, methods have been
found for circuit realization by straightforward means, which results in relatively
simple circuit topologies. There are only two equations which really need to be
transformed in the model, equations (2.42) and (2.49), as repeated in the following:
IK1 = 0:35
" 4(e0:04(V
m+85) ?1)
e0:08(Vm+53) +e0:04(Vm+53) +
0:2(Vm +23)
1?e?0:04(Vm+23)
#
; (3.1)
fim = ?(Vm +47)e?0:1(V
m+47) ?1
: (3.2)
Since the numerical range of Vm is about -90 mV to 40 mV [27], the above equa-
tions present zeros in the denominators in the working range of Vm; one is at
Vm = ?23 mV, and the other is at Vm = ?47 mV, and let us call them the
zero-denominator points. Note that the zeros in the denominators do not cause
poles because the corresponding numerators also have zeros at the same Vm. By
expanding the denominators using Tylor series around Vm = ?23 mV for IK1 and
Vm = ?47 mV for fim, the terms containing Vm can be canceled in the denomina-
tors and the numerators, and this makes the functions of IK1 and fim continuous
45
at the zero-denominator points, i.e. we have:
limV
m)?23
IK1 = 2:82; (3.3)
limV
m)?47
fim = 10: (3.4)
However, the VLSL design of the equations of IK1 and fim are very challenging,
because the circuits have to handle the zero-divided-by-zero case, and there is no
circuit that can directly represent the equations as a whole. Hence, we intentionally
remove the zero-denominator points with the equation transformation process, in
which an equation is represented by a difierent equation that gives values very
close to those calculated from the original equations in a required working range,
and this will be explained the following.
The transformation is performed using the Matlab Curve Fitting Toolbox. Due
to the exponential tendency in the original data, we decide to let the flt equations
still contain exponential terms, and specify the type of flt to customize the equa-
tions that favor the implementation of using emitter-coupled pairs (refer to section
3.3.2). The flt curves are also shown in Figure 3.2 and 3.3, which indicate the flts
are quite successful because we nearly can not distinguish the fltted curves from
the original data. The resulting transformed equations are:
IK1 = 2:742+ 3:632e?0:05966(V
m?32:77) +1
? 5:879e0:1177(V
m+88:23) +1
(3.5)
for equation (3.1), and
fim = 109:7+Vm ? 62:92e?0:05687(V
m+76:2089) +1
(3.6)
for equation (3.2).
The goodness of flt is evaluated with Sum of Squares Due to Error (SSE), R-
square, adjusted R-square, and Root Mean Squared Error (RMSE) [62], and Table
46
Vm  [mV]
IK1 [a0A/cm2]
Figure 3.2: Plot of transformed Ik1 equation and its original equation.
Vm  [mV]
?m
Figure 3.3: Plot of transformed fim equation and its original equation.
Table 3.1: Statistics of goodness of flt for equations (3.5) and (3.6)
Sum of Squares R-square Adjusted R-square Root Mean
Due to Error Squared Error
IK1 0.09225 0.9995 0.9995 0.02634
fim 5.621 1 1 0.2033
47
3.1 lists the flt statistics for equation (3.5) and (3.6). SSE is deflned as the sum
of squared residuals, which are the difierences between the response values (the
original data) and the predicted values (fltted curves). A SSE closer to 0 indicates
the flt is more useful for prediction. R-square is the square of the correlation
between the response values and the predicted values, and adjusted R-square is
based on R-square and adjusted with the residual degree of freedom, which is
deflned as the number of response values minus the number of fltted coe?cients
estimated from the response values. Both R-square and adjusted R-square take
values between 0 and 1, and a value closer to 1 means a better flt for adjusted R-
square, but not always for R-square. RMSE, also known as the flt standard error,
is an estimate of the standard deviation of the random component in the data.
Like SSE, a RMSE closer to 0 indicates a flt that is more useful for prediction.
As shown in Table 3.1, the goodness of the flt for Ik1 and fim is quite satisfactory,
with adjusted R-square being 0.9995 for Ik1, and 1 for fim.
3.2.2 Parameter Scaling for VLSI Design
Parameter scaling is needed to linearly convert the values of variables and some
constants to be in the ranges that are feasible to be represented with electrical
signals in circuits. As introduced in chapter 2, the units utilized in the Beeler-
Reuter model are \small", for instance, the unit for the membrane potential is mV,
and the unit for currents is ?A=cm2. Consequently, the values of variables in the
equations are normally big, for example, the values of membrane potential can be
up to about a hundred (here we care only the value, i.e. \a hundred", not its unit
mV, because VLSI implementation models mathematical equations and ignores
their physics meanings, and this decides that VLSI can be also used to simulate
48
non-electrical systems). Therefore, in order to describe such big numbers in the
VLSI circuits, parameter scaling is necessary to shrink the values. In addition, we
care only about the equations and the values of variables for circuit realization,
and ignore their physical meanings in the model, and thus the circuit design of the
cardiac cell model does not have to preserve the original electrical meanings of the
variables. For example, the calcium equilibrium potential Es does not have to be
voltage in the circuit, and can be represented with a VLSI current and scaled into
the normal magnitude of a current. Moreover, there are some variables, such as
ionic concentrations and gating variables, that are not electrical signals. Hence, in
order to represent these variables, parameter scaling is usually required to convert
their values.
As introduced in chapter 2, the cardiac cell model consists of a set of equations
that are inter-dependent, and no variables can be decided by a single formula.
Hence, the in uence of parameter scaling can be divided into two parts, one is the
efiect on the equation which deflnes the scaled variable, the other is the efiect to
the equation(s) that takes the variable as an input. For the flrst case, variables can
appear as totally independent terms in equations, like equation (2.42)?(2.45) that
deflne ionic currents, equation (2.47)?(2.59) that describe gate opening/closing
rates, and equation (2.60) that expresses calcium equilibrium potential; or, the
variables can be deflned in time derivative equations, like equation (2.40) which
deflnes Vm, equation (2.46) that formulates gating variables, and equation (2.60)
that describes the calcium concentration.
Let us abstract an equation that formulates a variable x to be:
x = F(z1;z2;::zN); (3.7)
where zi (i = 1;2:::N) are the variables that are determining x through the function
49
F. Suppose the equation that uses x (assuming there is only one equation that
takes x as an input) can be abstracted to be:
y = G(x;z1;z2;::zN); (3.8)
where y is another variable, and G is a function that relies on variable x and zi.
After performing scaling x0 = kx, where k is a constant, equation (3.7) becomes:
x0 = k ?F(z1;z2;::zN) = F0(z1;z2;::zN); (3.9)
and the expression for y is changed to the following as a result of scaling x:
y = G(x
0
k ;z1;z2;::zN) = G
0(x0;z1;z2;::zN): (3.10)
After scaling, reorganization, such as combining coe?cients, is usually needed to
construct the equations to be friendly for circuit implementing.
The above scaling technique may sound quite simple, however, the parame-
ter scaling process for circuit design requires comprehensive consideration of all
variables, constant terms, and coe?cients in the equations, and demands the ac-
knowledgment of the existence and limitations of available mathematical functions
in the form of circuits and apply the acknowledgment into scaling to flt the working
ranges of the circuits. To scale an equation set, like the Beeler-Reuter model, in
which more than twenty variables and equations are involved, a simulation of its
non-scaling ideal components with PSpice (or with Matlab SimuLink alternatively)
is usually necessary to carry out in order to achieve useful numerical ranges before
the actual scaling can be done.
The representatives of variables in circuits (either voltages or currents) are also
decided during the scaling process, because VLSI currents are normally small nu-
merically, for example in the magnitude of several 10?6 with unit A, and the values
50
of voltages are big, for instance in the magnitude of a few decimals with unit V.
The electrical representatives of some variables may be pre-determined by their
appearances in equations, relying on the convenience of their VLSI implementa-
tion. For example, we select calcium concentration [Ca]i to be a current, because
it appears in a logarithms in equation (2.59) and logarithms can be realized using
a bipolar transistor whose base-emitter voltage is decided by the logarithm of its
collector current. Another example, the equations of the opening/closing gate ra-
tios (equations (2.47)?(2.58)) have many exponential terms whose powers contain
Vm, and this leads to the desirability to let Vm to be a voltage in circuits due to
the circuit realization of an exponential being, again, via bipolar transistors.
The above discussion of parameter scaling does not include the scaling for
time. Time scaling can be treated absolutely independent of the scaling of other
parameters, since it only afiects the width of a signal on the time axis, and can
be easily adjusted by multiplying the capacitance values of all the capacitors that
realize time derivatives by a constant.
3.2.3 Initial Value Shifting
As introduced in chapter 2, the variable values associated with the resting state
of a membrane are taken as the initial conditions of the difierential system, and
their values are listed in Table 2.3. This means that a reset mechanism is required
to set the initial states of the devices which implement the time derivative. In the
case that capacitors are used to realize the time derivative in difierential equations,
the initial voltages across the capacitors need to be set. This requirement brings
complexity for circuits. It is necessary for the reset circuits to work as batteries to
provide desired voltages and currents for charging the capacitors at a reset stage,
51
and also to be able to be switched ofi (or not in uence) the capacitors when the
circuits are operating as an activated cell. Here we propose an initial value shifting
method to avoid the reset circuits. We shift the resting state of the difierential
system to the origin and make the initial conditions all zeros while keeping the
same waveform of Vm created by the model.
In the Beeler-Reuter model, there are eight variables that are time dependent,
these being the membrane potential Vm (equation (2.40)), six gating variables, x1,
m, h, j, d, and f (equation (2.46)), and a varying calcium concentration, [Ca]i
(equation (2.60)). Their equations are special cases of:
8>
><
>>:
b0dxdt = b1 +b2x;
x(t = 0) = x0;
(3.11)
where b0 is a constant, b1 and b2 can be constants or variables that do not explicitly
depend on time. In the case of Vm, the second term on the right side b2x does not
exist, i.e. b2 = 0. The initial value shifting takes place by replacing x with bx+x0,
that is, setting x = bx+x0, equation (3.11) then becomes:
8
>><
>>:
b0dbxdt = b1 +b2(bx+x0);
bx(t = 0) = 0:
(3.12)
Equation (3.11) and (3.12) can be both mapped into the circuit depicted in
Figure 3.4, by setting difierent parameters. For simplicity of discussion, here we
ignore the scaling issue, and assume the parameters in the equation are suitable
for VLSI realization. In Figure 3.4 the capacitor C, with a capacitance of b0, is
charged by a current source b1 and a controlled output current from G1, whose
transconductance is equal to b2. The output current from G1 is b2(x+Va), where
x is the voltage across C and Va is a bias voltage. The linear voltage to voltage
converter E, serving as the output stage of the whole circuit, takes the input x and
52
C=b0
b1
G1 
a
x
b2(x+Va)
-Va
Gain= b2
E 
Gain= 1
x+Va
Figure 3.4: Initial value shifting for difierential equations.
?Va and subtracts them to create the output. When Va is connected to ground,
i.e. Va = 0, the diagram represents equation (3.11) and the initial voltage of C is
required to be set to x0. The output of E is then x?0 = x. Following equation
(3.12), the initial condition shifting process is performed by setting Va to x0 in the
circuit. The initial voltage of C then can be left to 0, with the output of E being
unchanged: bx+Va = (x?x0)+x0 = x.
The initial value shifting method takes the advantage of the circuit realization
of the linear transconductance and linear voltage converter, which have a pair
of input voltages whose difierence linearly determines the output. It is worth
mentioning that the presented method is not limited to shifting the initial value
for the mathematical description with the form of equation (3.11). It can be
extended to more general situations as formulated with:
b0dxdt = F(x) (x(t = 0) = x0): (3.13)
The equation after the initial value shifting is then written as:
b0dbxdt = F(bx+x0) (bx(t = 0) = 0): (3.14)
53
3.2.4 Reformulated Beeler-Reuter Equations for VLSI De-
sign
After applying equation transformation, parameter scaling, and capacitance initial
value shifting on the Beeler-Reuter model, the resulting equations can be mapped
directly into circuits. The new set of equations have one-to-one correspondence to
equations (2.40)?(2.60), with (2.42) and (2.49) being replaced by (3.5) and (3.6).
To avoid redundancy, we list in Table 3.2 only the equations that relate the new
variables, denoted with a single quotation mark, with the original variables. These
equations result from the parameter scaling and initial value shifting, and can be
taken into the original model equations (except replacing (2.42) and (2.49) with
(3.5) and (3.6)) to create a new set of equations, which are the ones we use for
VLSI design. The new equations represent a difierential system that has all-zero
initial conditions. In Table 3.2, the second column shows the variable conversions,
and the third column lists the electrical representative of the variables in our VLSI
circuits.
3.3 Circuit Blocks of Function Units
In this section we will present the sub-circuits that implement the mathematical
functions needed by the Beeler-Reuter cardiac cell model. The determination
of the necessary functions is dependent on the mathematical description of the
original model formulations and also the parameter scaling process, in which the
electrical representatives of the signals are decided. The circuits of the function
units are summarized in Table 3.3, and most of them are based on existing circuit
architectures, with some being modifled in order to obtain required functions or
54
Table 3.2: Relations of original variables and the scaled and shifted variables, and
their electrical representative in circuits.
Variable Conversion Equation Circuit Representative
Vm = 100V 0m Voltage
Iion (includes = 106I0ion Current
IK1, Ix1, INa, and ICa)
z (includes x1, m, = 2?105z0 Current
h, j, d, and f)
fix1 = 103fi0x1 Current
flx1 = 103fl0x1 Current
fim = 107fi0m Current
flm = 107fl0m Current
fih = 2?105fi0h Current
flh = 2?105fl0h Current
fij = 2?104fi0j Current
flj = 2?104fl0j Current
fid = 104fi0d Current
fld = 104fl0d Current
fif = 103fi0f Current
flf = 103fl0f Current
Es = 108E0s Current
[Ca]i =[Ca]0i Current
55
Table 3.3: Summary of function circuits. Ci (i = 1;2:::N) are constant.
Function Description Circuit
Iout = C1Vin Linear voltage to current converter Figure 3.5(a)
Vout = Vin1 ?Vin2 Voltage adder/subtractor Figure 3.6
Vout = Vin Voltage bufier Figure 3.7
Iout = C1eC2(Vin+C3)+1 Sigmoid function Figure 3.12
Iout = C1e(C2Vin+C3) Exponential function Figure 3.13
Iout = C1eC2xeC3(Vin+C4)+1 Exponentialized sigmoid function Figure 3.12, 3.13
Vout = C1 ln(C2Iin) Logarithm function Figure 3.14
Iout = Iin1Iin2C1 Two inputs current multiplier Figure 3.15
Iout = Iin1Iin2Iin3:::Iin;N+1C1C2:::C
N
Multiple inputs current multiplier Figure 3.17
achieve higher accuracy in the circuit realization of the mathematical expressions.
The working ranges of the circuits are discussed in terms of the accuracy the
circuit implementation of the functions can reach. A scheme for implementing big
capacitors with small capacitors using a NPN-based circuit is introduced at the end
of this section. For the readers that are familiar with VLSI circuits, this section
can be skipped to continue the reading from section 3.4, in which the circuits for
realizing the cardiac cell model are presented.
The circuits presented here are not limited by a particular fabrication process
and can be implemented by any BiCMOS technologies. For the technologies with
smaller feature sizes, VLSI circuits can be more densely integrated and are more
power e?cient. For our heart model application, we give the flrst priority to
the accuracy of the circuit realization of mathematical functions, and hence we
prefer to using technologies with larger feature sizes, which have less high order
56
efiects introduced by short-channel MOS devices and thus are more feasible to
implement equations with their I-V characteristics. As an example, we select AMI
Semiconductor 1.5 ?m ABN technology to realize the cardiac cell model. AMI 1.5
?m ABN process is a n-well CMOS process with two metal layers and two poly
layers. It provides an NPN option and also can be used to make capacitors with
PiP (poly2 over insulator over poly).
The presented VLSI circuit of the cardiac cell model is essentially MOSFET.
Many circuit blocks in our VLSI design use the voltage-controlled-current proper-
ties of MOS devices, and take advantage of the thin insulating layers under the
gates and the resulting good isolation MOS devices provide to the preceding and
the succeeding circuits. We select current to represent signals in most situations
as shown in Table 3.2, because current allows performing addition and subtrac-
tion simply, and can represent positive and negative numbers with large ranges
easily. This is opposite to voltages, whose magnitudes are limited by the dynamic
range of circuits, which is highly dependent on the supply voltages. Therefore,
most arithmetic circuit blocks presented take currents as input and output. Cur-
rents also can be consumed by capacitors for integration. NPN devices are used
in our application to realize exponential functions. It is worth mentioning that
the proposed circuits with NPN transistors do not depend on the reverse satura-
tion currents of the devices, and hence the circuits can be transmitted to other
fabrication technologies with minimal modiflcation of the bipolar devices.
In the rest of this section, we will introduce the basic arithmetic function blocks
used in the VLSI design of the Beeler-Reuter model. The function circuits can be
organized into three categories: linear function circuits, exponential circuits, and
multipliers. Table 3.3 summaries the functions of the presented circuits.
57
The following presentation will frequently use the equations that describe the
I ? V characteristics of MOS and NPN transistors. To avoid verbosity, some
commonly used symbols and notational conventions are listed alphabetically below:
? Cox is the gate oxide capacitance per unit area.
? Ici (or Ic) is the collector current of NPN transistor Qi (i=1,2...).
? IMi is the drain-to-source current of transistor Mi (i=1,2...).
? Is0 is the reverse saturation current of a NPN transistor.
? Kni (or Kn, Kpi, Kp, K) is the fabrication-dependent parameter of NMOS
device Mi (i=1,2...): Kni = ?nCox2 WiLi .
? ?n (or ?p) is the mobility of electrons (holes).
? Li (or L) is the gate length of transistor Mi (i=1,2...).
? ?ni (or ?n, ?pi, ?p, ?) is the channel-length modulation parameter of NMOS
(or PMOS) transistor Mi (i=1,2...).
? Vgsi is the gate-source voltage of transistor Mi (i=1,2...).
? Vbei is the base-emitter voltage of NPN transistor Qi (i=1,2...).
? Vtni (or Vtpi, Vtn, Vtp, Vt)is the threshold voltage of NMOS (or PMOS) tran-
sistor Mi (i=1,2...).
? VT is the thermal voltage.
? Wi (or W) is the gate width of transistor Mi (i=1,2...).
58
3.3.1 Linear Function Circuits
With currents being signal representatives, the linear operations, such as addition,
subtraction, and multiplication with a constant, can be easily realized using current
mirrors and wire connections, governed by the Kirchhofi?s law. In the following, we
will introduce the circuits that perform linear operations based on voltage inputs
or (and) outputs.
Linear Voltage to Current Converter (VCC)
One of the techniques to realize a linear voltage to current converter (VCC) is
using difierential pairs as shown in Figure 3.5(a) [63]. M1 and M2 are the NMOS
input stage, whose tail current I0 is provided by the current sink transistor M9.
The current distribution in the two input transistors M1 and M2 is controlled by
their gate voltages, and these two currents are transmitted to the output stage by
current mirrors that are composed of M3 ? M8. The output Iout current is the
difierence between the two currents passing through transistors M6 and M8. By
varying the sink current I0, and the ampliflcation factors of the current mirrors,
the transconductance of the difierential pair can be adjusted.
Assume that M1 and M2 are identical and both working in the saturation
region, neglecting the channel length modulation, the source currents of M1 and
M2 are expressed by:
IM1 = K(Vgs1 ?Vt)2 = K(?V +Vgs2 ?Vt)2; (3.15)
IM2 = K(Vgs2 ?Vt)2; (3.16)
where ?V is the difierence between the two voltage inputs:
?V = V+ ?V? = Vgs1 ?Vgs2: (3.17)
59
(a) (b)
M9 Vbb
Vss
Vdd
M1
M3 M4
M2V+ V-
I0
M5 M6
M7 M8
Iout
M9 Vbb
Vss
Vdd
M1
M3 M4
M2V+ V-
M5 M6
M7 M8
I- I+
M10
M11
M12
M13
Figure 3.5: CMOS difierential pair works as a linear voltage to current converter.
(a)Difierential pair circuit. (b)Altered difierential pair that works as a linear re-
sistor.
From equations (3.15) and (3.16), we can derive the following equations by using
I0 = IM1 +IM2:
IM1 = 12
0
@I0 +K?V
s
2I0
K ??V
2
1
A; (3.18)
IM2 = 12
0
@I0 ?K?V
s
2I0
K ??V
2
1
A: (3.19)
Therefore, the transconductance of the circuit is expressed by:
gm = IM8 ?IM6?V = AIM1 ?IM2?V = AK
s
2I0
K ??V
2; (3.20)
where A is the ampliflcation factor of the current mirrors composed of M3?M8,
and M4 and M6. From equation (3.20), we can derive the range of ?V for a given
linearity requirement on the transconductance by a Tylor series expansion of the
square root:
?V 2 < ?4I0K ; (3.21)
60
where ? is the maximum allowable transconductance difierence in percentage, i.e.
? = max((gm0?gm)=gm0)?100% where gm0 is the transconductance at ?V = 0,and
gm is an arbitrary transconductance.
Ignoring the slight changes in I0 caused by the variation of Vds9, the non-
linearity of the transconductance relies on the term ?V 2 in equation (3.20). Hence,
according to the equation, in order to achieve a good linearity of the transconduc-
tance, the magnitude of the difierential input ?V needs to be kept much smaller
than 2I0K , or at the same ?V, I0 can be enlarged, under the restriction that all
transistors are still in the saturation region.
The circuit of the difierential pair can be used to implement linear resistors, as
illustrated in Figure 3.5(b). The transistors and wires drawn with dotted lines are
newly added upon the circuit shown in Figure 3.5(a). The added transistors make
the circuit topology of the left side and the right side perfectly symmetrical, and
thus I+ = ?I?. When V+ > V?, the current  owing into the V+ port is provided
by M6 and M8, and the current  owing out of the V? port is sourced from M10
and M11. The working condition of the circuit in Figure 3.5(b) is that the inputs
V+ and V? are within the range which makes the output transistors M6, M8, M10
and M11 working in the saturation region.
Linear Current to Voltage Converter (CVC)
The circuit of a linear current to voltage converter, shown in Figure 3.6, is proposed
in [64]. It is composed of two NMOS and two PMOS transistors. M1 and M2 are
identical in size. M3 and M4 are also identical in size, and are connected like a
current mirror. Va is a bias voltage, which determines the current in M1. When
the input current Iin is zero, the currents in all the transistors are the same, and
61
Vsgp
Vdd
Vout
Iin
I3
I1
I4
I2Va
Vb
M1 M2
M3 M4
Figure 3.6: Linear current to voltage converter.
thus Vout = Vb when there is no load drawing current from the output. Let us use
Vout0 to denote the output voltage associated with Iin = 0, and IQ to represent the
drain current  owing in M1 when Iin = 0.
If there is a small negative input current, Vb can be taken as nearly unchanged
due to the small current Iin. Since the current I4 is mirrored from I3, and I2 is
mirrored from I1, when the output does not take any current, the output voltage
Vout needs to increase signiflcantly to cancel out the difierence between I2 and I4
using the channel length modulation efiect. As we will see in the following, the
change in the output voltage is almost linearly decided by the input current.
For the transistors M1 and M3, the following equations describe their currents:
I3 + Iin = I1; (3.22)
I1 = Kn ?(Va ?Vtn)2(1+?nVb); (3.23)
I3 = Kp ?(Vdd ?Vb ?Vtp)2(1+?p(Vdd ?Vb)): (3.24)
Similarly the currents in M2 and M4 are expressed by the following equations:
I2 = I4; (3.25)
62
I2 = Kn ?(Va ?Vtn)2(1+?nVout); (3.26)
I4 = Kp ?(Vdd ?Vb ?Vtp)2(1+?p(Vdd ?Vout)): (3.27)
From equation (3.22)?(3.27), we can derive the relation of Vout vs. Iin as follows:
Iin = (Vb ?Vout)Kn(Va ?Vtn)
2(?n +?p +?n?pVdd)
1+?p(Vdd ?Vout) : (3.28)
If we neglect the change of Vb caused by the small input Iin, that is, take Vb = Vout0,
and treat ?p(Vdd?Vout) as much less than 1 and, thus, ignore it in the denominator
of equation (3.28), the change of the transimpedance, i.e. the relation of the change
of the output and the input, then becomes purely linear, which is expressed by:
?Vout
Iin = ?
1
Kn(Va ?Vtn)2(?n +?p +?n?pVdd) (3.29)
? ? 1I
Q(?n +?p)
; (3.30)
where ?Vout is the change in Vout, caused by Iin and deflned as: ?Vout = Vout?Vout0.
The nonlinearity of the ?Vout to Iin ratio comes from the terms ?p(Vdd ?Vout)
and Vb in equation (3.28). ?p(Vdd ?Vout) can cause a variance of the ?Vout-to-Iin
ratio on the order of 1% assuming ?p ? 0:01. Vb can be derived from equation
(3.22)?(3.24) and is determined by the following equation:
Kp?(Vdd?Vb?Vtp)2(1+?p(Vdd?Vb))+Iin = Kn?(Va?Vtn)2(1+?nVb): (3.31)
Take ?n = ?p ? 0 in equation (3.31), Vb is expressed by:
Vb = (Vdd ?Vtp)? 1qK
p
q
Kn(Va ?Vtn)2 ?Iin: (3.32)
According to equation (3.32), Iin needs to be small compared to Kn(Va ?Vtn)2 in
order to diminish the dependence of Vb on Iin. Kn(Va?Vtn)2 can be approximated
as the current in M1 when Iin = 0, and hence the condition for a good linearity of
63
the ?Vout-to-Iin ratio is that Iin needs to be kept small compared to IQ, i.e. the
current of M1 when Iin = 0.
Vout0 usually varies with the selection of difierent sizes of the transistors for a
current to voltage converter (CVC), and the transistor sizes also decide the values
of the ?Vout-to -Iin ratio. It is not easy to adjust the circuit to achieve both the
specifled ?Vout-to-Iin ratio and Vout0 with a restricted value. Therefore, in the
heart VLSI circuits, CVCs are employed in pairs, and this will be discussed in
section 3.3.2.
A working condition of the CVC circuit is that the output does not take any
current, due to the great sensitivity of the output voltage to the current changes
in M2 and M4. In order to provide a current to a succeeding circuit that utilizes
the linearly changed output voltage from a CVC, it is necessary to add a voltage
bufier to copy the output voltage and also to deliver a certain among of current.
The voltage bufier is introduced in the following.
Voltage Bufier
An ordinary CMOS voltage bufier [65], shown in Figure 3.7(a), is composed of
complementary source followers. On the input side, the current of M1 and M3 is
decided by the bias current In0 and Ip0, and In0 = Ip0. Because M2 is set identical
to M1, and M4 is to M3, the currents In0 and Ip0 are fully mirrored on the output
side if no output current is drawn, and hence Vout = Vin. The circuit has very low
voltage ofiset, with the voltage ofiset mainly due to transistor mismatch [66].
In the VLSI design of the heart model, the voltage bufier services as a stage that
isolates the next circuit block, which takes voltage as an input but also consumes
some amount of current depending on the input voltage, from the previous stage.
64
(a) (b)
Vin Vout
Vdd
Vss
Ip0
In0
M3
M1 M2
M4
Vdd
Vss
M1 M2
M3
M4
M6
M5
Vin Vout
(c)
Vdd
Vss
M1 M2
M3
M4
M6
M5
Vin Vout
Figure 3.7: Voltage bufier. (a)Traditional voltage bufier. (b)Proposed voltage
bufier for the heart implementation. (c)Another version of proposed voltage bufier,
modifled for reducing drain-source voltage difierence of M1 and M2.
For a circuit like the CVC presented in the previous section, its output voltage can
be easily distorted by its output current. Therefore, the bufier is required to take
no current from the previous stage, or equivalently the two bias currents In0 and Ip0
need to be matched very well in order to prohibit current  owing through the input
terminal. Since it is unavoidable for the non-ideal bias currents to vary a little bit
with the changing input voltage, the requirement of no input current is obtained
by letting the changes of In0 and Ip0 be in the same direction and also by the same
amount which brings much complexity to the transistor circuits. This makes the
ordinary voltage bufier presented in [65] not applicable to our implementation. In
order to achieve high enough input resistance, the gate of a NMOS device, which
does not draw any current, is selected as the input stage of the proposed bufier.
The topology of the circuit is depicted in flgure 3.7(b). Transistors M3 ? M6,
working in the saturation region, compose a cascode current mirror, which keeps
the currents going through M1 and M2 the same. The output transistor M2
65
is identical to M1. When the efiect of the drain-source voltage on the current
is ignored, the gate-source voltages are the same for M1 and M2 due to their
identical currents, and this results in Vout = Vin. Suppose the current required by
the next stage is much less than the current provided by M6, the in uence of the
output current to the output voltage is then very small and can be neglected. This
assumption can be satisfled by selecting the parameters of M1 to let the current
in M6 much bigger than the required current output.
The I ?V characteristics of M1 and M2 can be described as follows, with the
efiect of their drain-source voltages taken into account:
IM1 = K(Vin ?Vss ?Vt)2(1+?Vds1) (3.33)
= K(Vin ?Vss ?Vt)2(1+?Vds2 +??Vds);
IM2 = K(Vout ?Vss ?Vt)2(1+?Vds2) (3.34)
= K(Vin +?Vout ?Vss ?Vt)2(1+?Vds2);
where ?Vout is the output error caused by the channel modulation efiect, ?Vout =
Vout ?Vin, and ?Vds is the difierence of the drain-source voltages of M1 and M2
is ?Vds = Vds1 ?Vds2. Ignoring the current difierence of IM1 and IM2 introduced
by the non-idealism of the current mirror, after combining equations (3.33) and
(3.34), and replacing Vds2 (i.e. Vgs2) with Vin ?Vss, we express ?Vout in terms of
?Vds:
?Vout = (Vin ?Vss ?Vt)(
s
1+ ??Vds1+?(V
in ?Vss)
?1); (3.35)
which can be expanded with Taylor approximation into:
?Vout ? 12(Vin ?Vss ?Vt)? ??Vds1+?(V
in ?Vss)
(3.36)
? 12(Vin ?Vss ?Vt)???Vds:
66
Vdd
V1
V2
Vout
M1
M2
Vdd
CurrentMirror
CurrentMirror
Vout
M2
V1
V2
Vss
Vbias
M1
M4
M5
M3
Va
Vb
(a) (b)
Vc
Figure 3.8: Voltage subtractor. (a)Maundy voltage subtractor. (b)Improved volt-
age subtractor.
Equation (3.36) expresses the voltage ofiset in the output in terms of the dif-
ference in the drain-source voltages of M1 and M2. The reduction of the average
of ?Vds can e?ciently minimize the voltage ofiset. Figure 3.7(c) shows another
version of the proposed voltage bufier, which e?ciently reduces the magnitude of
?Vds. The currents in the NMOS transistors M3 and M6 are the same as the
currents in M4 and M5, and, thus, the gate-source voltages of M3 and M6 are
equal. Therefore, the drain-source voltages of M1 and M2 are equal. Ibias is used
to constrain the currents in M1 and M2.
Voltage Subtractor and Adder
Figure 3.8(a) illustrates the idea of a voltage subtractor presented by Maundy in
[55], where the circuit is introduced to make low-voltage multipliers. The sizes M1
and M2 are equalized, and both transistors work in the saturation region. When
the channel length modulation efiect is not considered, the following equations are
67
satisfled:
8
>>>>
>><
>>>
>>>:
IM1 = Kp(Vgs1 ?Vtp)2 = K(Vdd ?V1 ?Vtp)2
IM2 = Kp(Vgs2 ?Vtp)2 = K(Vout ?V2 ?Vtp)2
IM1 = IM2 = I
(3.37)
) Vout = Vdd ?V1 +V2: (3.38)
Hence this circuit operates like a subtractor. Multiples of this circuit can be cas-
caded to construct mathematical functions with mixed additions and subtractions
[55].
The performance of Maundy?s voltage subtractor is degraded due to the channel
length modulation efiect. If the channel length modulation efiect is taken into
account, the equations for M1 and M2 become the following:
8
>><
>>:
IM1 = Kp(Vdd ?V1 ?Vtp)2(1+?p(Vdd ?Vout));
IM2 = Kp(Vout ?V2 ?Vtp)2(1+?pVout):
(3.39)
Assume there is a drop of amount ?V1 on input V1, and input V2 stays unchanged.
With IM1 = IM2, according to equation (3.39), Vout does not increase the same
amount as V1 reduces, and the difierence between the increase and that of the
ideal case is on the same order as ?p or bigger, decided by the values of V1, V2, and
Vdd. We propose two methods to improve Maundy?s subtraction circuit, based on
considering the two transistors M1 and M2 separately. M1 is the key transistor
that decides the current of the circuit. In the non-ideal case, the output voltage
can in uence the current. If M1 can be replaced by a perfect current source, and
the value of the current is independent of the output, the inaccuracy introduced
by the drain-source voltage of M1 can be eliminated. For M2, if a feedback circuit
can be added to control the drain voltage of M2 to make it follow Vout, the channel
length modulation efiect caused by the changing Vds2 can be diminished.
68
Figure 3.8(b) depicts the circuit of the improved voltage subtractor. In the
circuit, M1 and M2 serve the same functions as in Figure 3.8(a). Two current
mirrors are added to isolate the input V1 from the output. The lower current
mirror is required to employ large size transistors to minimize the change of Vc
caused by the change of V1. Cascode current mirrors are also recommended to
reduce the channel modulation efiect on the current. The transistors M3 ? M5
consist of a feedback loop that makes Va follow Vout. M4 and M5, the same in
size, incorporate a subtraction circuit just like the one shown in Figure 3.8(a),
except NMOS transistors are used this time. Hence the following equation is
approximately satisfled following equation (3.38):
Vout ?Vb ? Vbias ?Vss = constant: (3.40)
According to the equation, when Vout changes, Vb changes in the same direction
with about an identical amount. If M3 is big in size, and the changed current
caused by V1 results in little change on Vgs3, then Va can be taken as following Vb.
As a result, Va follows Vout, and the drain-source voltage efiect of M2 is mitigated
by the feedback circuit comprised of M3?M5.
Another implementation of a voltage subtractor (or adder) is from [67], called
the pool circuit, which, as depicted in Figure 3.9, consists of two difierential pairs.
The currents of M1 and M2 , derived with the same method for equation (3.23),
are described be the following equations:
IM1 = 12
?
IB1 +Kn?V
s2I
B1
Kn ??V
2
!
; (3.41)
IM2 = 12
?
IB1 ?Kn?V
s2I
B1
Kn ??V
2
!
; (3.42)
where ?V is the difierence of the gate voltages of M1 and M2:
?V = V1 ?Vout: (3.43)
69
Vdd
M3 M4
Vss
M1 M2
Vdd
M7 M8
Vss
M5 M6
IB1 IB2
I1 I2
V1 V3 V2
Vout
Figure 3.9: Pool circuit works as voltage adder/subtractor.
Since M3 and M4 is connected as a current mirror, the current in M4 is forced to
equal that in M1. Therefore, the output current from the left difierential pair is
the difierence of the currents in M1 and M2:
I1 = IM4 ?IM2 = IM1 ?IM2 (3.44)
= Kn(V1 ?Vout)
s2I
B1
Kn ?(V1 ?Vout)
2:
Similarly, I2 can be derived from the difierential pair on the right side:
I2 = IM7 ?IM5 = IM6 ?IM5 (3.45)
= Kn(V2 ?V3)
s2I
B2
Kn ?(V2 ?V3)
2:
If IB1 = IB2, with the use of the relation I1 = ?I2, combining equations (3.44)
and (3.45) yields the expression for the output voltage in terms of V1, V2 and V3:
Vout = V1 +V2 ?V3: (3.46)
Therefore, the pool circuit can be used as a subtractor with V1 or V2 being flxed,
or a adder with V3 being flxed. Compared to Maundy?s voltage subtractor/adder,
the pool circuit is more  exible to the inputs and can work alone as an operational
70
Vdd
M3 M4
Vss
M1 M2
Vdd
M7 M8
Vss
M5 M6
IB1 IB2
V1 V3 V2
VoutVa Vb Vc
Figure 3.10: Improved pool circuit.
unit, whereas Maundy?s circuit needs to be cascaded to cancel Vdd in equation
(3.38), and is limited by the saturation region of M1 and M2 in Figure 3.8 and
thus requires V1 > V2 +2Vtp, derived from:
8>
><
>>:
V1 > Vout +Vtp;
Vout ?V2 > Vtp:
(3.47)
The linearity of the pool circuit is limited by the channel length modulation
efiect, and the efiect in uences both the left side and the right side of the circuit.
We propose a new version of the pool circuit, shown in Figure 3.10, to alleviate the
channel length modulation efiect. The wire connections among M3, M4, M7, and
M8 are changed to make the currents of M3 and M4 mirrored from the currents of
M7 and M8 separately. Hence, IM1 ?IM2 = IM5 ?IM6 is satisfled, and this leads
to the same expression of the output voltage as equation (3.46). The advantage
of the new circuit is that the right part of the circuit is totally independent of
the output. M3, M4, M7, and M8 can be applied with large size transistors to
minimize their source-gate voltage changes caused by the inputs, and make Vb and
Vc roughly equal.
71
3.3.2 Exponential Functions
We select NPN transistors to implement the exponential functions needed in the
heart model. The exponential property of a NPN transistor comes from its I ?V
characteristics in the forward-active region that can be described by:
Ic = Is0e
Vbe
VT ; (3.48)
where Is0 is the reverse saturation current, and VT is the thermal voltage. In order
to obtain an exponential function, a voltage is taken as the input and connected
to the base of a NPN transistor. The current of the collector is then exponentially
decided by the base voltage. NPN transistors can also be used to calculate loga-
rithms, if the base-emitter voltage is taken as output when the collector current is
treated as the input.
Difierent from the gates of MOS transistors, NPN devices demand currents on
the bases. When the base serves as an input terminal, we can no long consider
the NPN well isolated from the preceding circuit and neglect its in uence. In
addition, the base voltage needs to be controlled to be less than about 0.7 V to
avoid large current and high power consumption. These constraints add new design
requirements absent from MOS circuit designs. In the following, we will discuss
circuits based on NPN transistors. It is worth to mention that in our circuits, the
equations of the input-output relations of the function circuits are independent of
Is0. We use a single size for all NPN transistors and add compensation circuits to
cancel the Is0 terms in the equations. The advantage is that the function circuits
can be more easily exported to difierent technologies which may have variant device
parameters, and still keep the same input-output relationship.
72
M1 Vbias
Vss
Vdd
Q1 Q2
M2 M3
V+ V-
I0
Va
M4 M5
I+ I-
Figure 3.11: Emitter-coupled circuit.
Exponential Circuit Using Difierential Pairs
The most frequently used exponential function in the Beeler-Reuter model has the
form:
F(x) = C1eC
2(x+C3) +1
; (3.49)
where C1; C2; C3 are constants. This equation is very similar to the current
expressions for the emitter-coupled pair circuit, shown in Figure 3.11. M1 provides
a current I0 whose value is decided by a bias voltage Vbias, and M2 and M3 serve
as active loads of Q1 and Q2. M4 and M5 copy the currents in Q1 and Q2 and
send them to the outputs. The transistors are matched on the left side and the
right side. The currents of Q1 and Q2 are described by the following equations:
Ic1 = Is0e
V+?Va
VT ; (3.50)
Ic2 = Is0e
V??Va
VT ; (3.51)
I0 ?fiF = Ic1 +Ic2; (3.52)
where fiF is deflned by the transistor forward current gain flF and is approxi-
mately equal to 1: fiF = flFfl
F+1
? 1. Equation (3.50)?(3.52) leads to the following
73
Emitter Coupled Pair
x Vin     Vout
V Buffer1G1 
V-
IoutV+
-C3 Vin     Vout
V Buffer2G2 
Iout
Ia Vb
Io
H1 
H2 
Vc
Vd
Ve
Figure 3.12: Circuit that implements equation (3.49).
expressions:
Ic1 = I0
e(??VVT ) +1
; (3.53)
Ic2 = I0
e(?VVT ) +1
; (3.54)
where we take fiF = 1, and ?V is the difierence of the two inputs ?V = V+?V?.
The great similarity between equations (3.53) (3.54) and equation (3.49) in-
dicates that the emitter-coupled pair circuit can be used to provide the function
described in equation (3.49). The main task to adopt it arises from the difierences
of the exponential powers in the equations, which require linear transformations
between x and ?V. To solve the problem, linear modules introduced in section
3.3.1 are applied to perform the conversion, as illustrated in Figure 3.12. G1 and
G2 are identical voltage to current converters (VCC), whose transconductance is
g1. They can be implemented with the difierential circuit introduced previously in
Figure 3.5. H1 and H2 are identical current to voltage converters (CVC), whose
transimpedance is h1. They can be realized with the circuit presented in section
3.3.1, Figure 3.6. The voltage Vb is expressed by the following equation:
Vb = Iah1 +bH1 = h1(x?g1 +bG1)+bH1 (3.55)
= x?g1h1 +bG1h1 +bH1;
74
where bG1 and bH1 are the output ofisets of G1 and H1 when their inputs are zero.
Similarly, Vc can be expressed by:
Vc = ?C3g1h1 +bG1h1 +bH1: (3.56)
Two voltage bufiers replicate Vb and Vc to be Vd and Ve and provide currents to the
bases of the NPN transistors in the emitter-coupled circuit. Following equations
(3.53) and (3.54), Iout can be expressed by:
Iout = I0
e?
g1h1(x+C3)
VT +1
; (3.57)
or
Iout = I0
e
g1h1(x+C3)
VT +1
; (3.58)
depending on which side of the emitter-coupled pair the current comes from. Let
I0 equal to C1, and g1h1V
T
equal to C2, the equation of Iout then becomes the same
as equation (3.49).
The advantages of using double VCCs and CVCs are the ofisets in their outputs
can be canceled. From the other point of view, if only a single VCC and CVC are
used, the ofiset terms bG1h1+bH1 in equation (3.55) will be required to implement
C3. This makes the design very challenging, because C2 and C3 in equation (3.49)
are, though independent of each other, now implemented to be determined by
the parameters of H1, and h1 and bH1 can not be adjusted independently in the
introduced CVC circuit.
Single Exponential Term
In this section, we introduce how we realize the circuit to calculate a single expo-
nential term, which can be described with the equation:
F(x) = C1e(C2x+C3); (3.59)
75
Vin
G1 
G2 
Ia
Pool 
Circuit 
V+1
V+2
V-
Vout
Iin
Vin     Vout
V Buffer I
out
Q1 Q2
Vb
VcVd VeH1 
H2 
Vbias 
Figure 3.13: Circuit implementation of equation (3.60).
where C1; C2; C3 are constants. The above equation is equivalent to:
F(x) = C4e
C5(x+C6)
VT ; (3.60)
where C4; C5; C6 are constants and restricted by:
8
>><
>>:
C5
VT = C2;
C4e
C5C6
VT = C1eC3:
(3.61)
The similarity between equation (3.60) and equation (3.48) allows the circuit imple-
mentation of equation (3.60) using the exponential properties of NPN transistors.
Looking at the circuit shown in Figure 3.13, G1 is a linear transconductance,
and H1 is a linear transimpedance. G1 and H1 convert the input Vin to Vb, which
can be expressed by:
Vb = Iah1 +bH1 = h1(Ving1 +bG1)+bH1 (3.62)
= Ving1h1 +bG1h1 +bH1;
where g1 and h1 are the transconductance and transimpedance of G1 and H1, and
bG1 and bH1 are their ofisets when the inputs of G1 and H1 are zero. Suppose G2
and H2 are exactly the same as G1 and H1, taking Vin = Vbias into equation (3.62)
we then obtain the equation for Vc using the equation for Vb:
Vc = Vbiasg1h1 +bG1h1 +bH1: (3.63)
76
The second pair of transconductor and transresistor G2 and H2 are used to cancel
the ofiset bG1h1 +bH1 as we will see soon.
Q1 and Iin compose a logarithm function, and the voltage Vd is decided by
the current Iin. Assuming the forward current gain flF of Q1 is large, so that the
current  owing into the base can be ignored and Iin ? Ic;Q1, Vd is then described
by the following formula, derived from equation (3.48):
Vd = VT ln IinI
s0
: (3.64)
A pool circuit, introduced in section 3.3.1, puts Vb, Vc, and Vd together following
equation (3.46), and yields an output voltage expressed by:
Ve = Vb ?Vc +Vd: (3.65)
Replacing Vb, Vc, and Vd with equation (3.62), (3.63), and (3.64) generates:
Ve = (Vin ?Vbias)g1h1 +VT ln IinI
s0
: (3.66)
Note that the ofiset in G1 and H1 are canceled by that of G2 and H2.
Ve is sent to the output transistor Q2 through a voltage bufier to control the
output current Iout. The voltage bufier, presented in section 3.3.1, reproduces the
voltage Ve in its output and provides current to the base of Q2. Hence the collector
current of Q2 can be described in terms of Ve using the I?V characteristics of an
NPN transistors:
Iout = Is0eVeVT (3.67)
= Is0e
(Vin?Vbias)g1h1+VT ln IinIs0
VT
= Iine
g1h1(Vin?Vbias)
VT :
Substitute Iin with C4, and substitute g1h1 with C5, Vbias with ?C6, and Vin with
x in equation (3.67), we obtain the same expression as equation (3.60). Hence,
77
the circuit in Figure 3.13 is a realization of equation (3.60). The circuit does not
require C5 to be positive, because Vin and ?C6 can be connected to the negative
input terminal of G1, and thus in equation (3.67) there is an additional \-" sign
in front of g1h1 for this case. Parameter scaling is involved when C4 can not be
directly represented by Iin in a circuit.
The circuit is restricted by the base voltage of Q2, or equivalently Ve, with tak-
ing account of the power on Q2. Since the collector currents increase dramatically
with increasing base-emitter voltage, the current can be amazingly high when Ve
reaches above a certain value. To avoid this situation, Ve is normally constrained to
be under 0.7 V. The restriction can be expressed mathematically by the following
equation, using equation (3.66):
(Vin +C6)g1h1 +Vbe1 < 0:7; (3.68)
) (Vin +C6)g1h1 < 0:7?Vbe1: (3.69)
The restriction can be satisfled by selecting proper C4 and C6 under the constraint
described by equation (3.61). Note the input voltage of the \V bufier" Ve is not
required to be higher than 0. When Ve < 0, the collector current of Q2 is nearly
(if not absolutely) 0.
The circuit of Figure 3.13 can be used to replace the I0 in Figure 3.12 to realize
functions of the form:
F(x) = D1e
D2x
eD3(x+D4) +1; (3.70)
where D1 ?D4 are constants.
Logarithm Function
The circuit that is used here to realize logarithm functions is depicted in Figure
3.14. The NPN transistors Q1 and Q2, connected like diodes, take Ix and Ic as
78
Ix
Q2Q1
Ic G 
Iout
Vb
Va
Figure 3.14: Circuit implementation of logarithm function.
their collector currents, with Ic, a constant, working as a bias, and Ix being taken
as a variable. Va and Vb are determined by the base voltages and described by:
Va = VT ln IxI
s0
; (3.71)
Vb = VT ln IcI
s0
: (3.72)
The output current of the circuit Iout is controlled by Va and Vb through a linear
voltage to current converter G1 (refer to section 3.3.1), and is expressed with:
Iout = g1(Va ?Vb) = g1VT lnIx ?g1VT lnIc = g1VT ln IxI
c
; (3.73)
where g1 is the transimpedance of G1. Therefore the circuit in Figure 3.14 is the
implementation of functions of the form:
F(x) = C1 lnC2x; (3.74)
where C1 and C2 are constants, if we let g1 = C1=VT and let Ic = 1=C2. Current
mirrors and parameter scaling can be used to handle the difierence between C1
and g1VT.
3.3.3 Multipliers
Analogue multipliers are important building blocks in analog signal processing.
The literature [56] gives a review on the multiplier circuits that are based on the
79
Gilbert cell. Those circuits take voltages as inputs, and output currents. In the
VLSI design of the heart model, because in most cases signals are represented
by currents (refer to Table 3.2) due to the easy manipulation of current signals
compared to voltage signals, we focus on current multipliers here. Various CMOS
current multiplier circuits have been published, with transistors operating in either
the subthreshold region or the strong inversion region [51]?[54]. Many of them,
like the circuits introduced in literature [51] and [52], are composed of both NMOS
and PMOS transistors and require their parameters to be equal, i.e. Kn = Kp
or (and) Vtn = Vtp. This restriction is usually not easy to satisfled because of
the limitation of fabrication technologies. In this section, we present a current
multiplier proposed by Tanno. After that, we will introduce a multiplier based on
bipolar transistors for implementing cascaded multiplications. A scheme of making
big capacitors based on the bipolar multiplier will be presented last.
Tanno Multiplier
The Tanno multiplier is developed based on the quarter square technique, which
is deflned by:
Io = (Ix +Iy)2 ?(Ix ?Iy)2 = 4IxIy: (3.75)
The technique uses addition and subtraction operations, as well as squaring oper-
ations. The requirement decides that current-squaring circuits are needed in the
multiplier circuits.
The squaring circuit used by the Tanno multiplier is depicted in Figure 3.15(a)
[54]. The transistors work in the strong inversion region. M1 and M2 are identical,
and they incorporate the bias circuit, in which Vb is decided by the bias current
IB. M3 and M4 are also identical in size. The sources and gates of M1?M4 are
80
(a) (b)
IB
M1
M2
M3
M4
Io
Iin
Bias Block
Vb
Va
IB
M1
M2
Bias Block
M3
M4
M5
M6
M7
M8
M9
M10
Ix Ix IyIy
Ia Ic IdIb
Iout= Ia+ Ib- Ic- Id
Current Subtractor Iout
Vb
Va
Vc
Figure 3.15: Tanno current multiplier. (a)Basic block of Tanno multiplier.
(b)Tanno multiplier.
connected into a loop, and hence their gate-source voltages follow:
Vgs1 +Vgs2 = Vgs3 +Vgs4: (3.76)
Their gate-source voltages can be expressed in terms of their currents:
Vgs1 = Vgs2 =
s I
B
KM1 +Vt; (3.77)
Vgsi =
s I
i
KM3 +Vt (i = 3;4); (3.78)
where we assume the body-source voltages of the transistors are the same, so that
the threshold voltages are all equal to Vt. When KM1 = KM3, we have:
2
q
IB =
q
IM3 +
q
IM4: (3.79)
From Figure 3.15(a), the current in M3 is:
IM3 = Io = IM4 ?Iin: (3.80)
81
Combining equation (3.79) and (3.80) yields:
IM3 +IM4 = 2IB + I
2
in
8IB; (3.81)
and thus:
Io = (Iin ?4IB)
2
16IB : (3.82)
The sizes of M1 and M2 are not necessary to be the same as for M3 and M4.
Assume KM1 = mKM3, then Io of Figure 3.15(a) is described as:
Io = (mIin ?4IB)
2
16mIB : (3.83)
The Tanno multiplier is composed of one bias block, four output blocks, and a
current adder/subtractor, as depicted in Figure 3.15(b). The current Ia, Ib, Ic, and
Id can be expressed by equation (3.83) with the replacement of Iin with 0, Ix +Iy,
Ix and Iy respectively. The output of the multiplier circuit is then described by
the following:
Iout = Ia +Ib ?Ic ?Id = m8 ? IxIyI
B
: (3.84)
The current subtractor in Figure 3.15(b) can be realized with the circuit shown in
Figure 3.16, in which Iin1 and Iin2 are copied with current mirrors and changed
directions so that Iout is Iin1 ?Iin2.
A limitation of the Tanno multiplier relies on the assumption that the threshold
voltages of M1?M4 are identical in Figure 3.15(a) so that the threshold voltage
Vt can be canceled out to obtain equation (3.79). This assumption may not be true
if the transistor bodies are not connected like in Figure 3.15(a), because Vc changes
along Iin, and this makes the source voltage of M3 vary and, thus, the threshold
voltage of M3 can be difierent from that of M1 due to the body efiect. For the
AMI 1.5 ?m ABN process we select for the VLSI design, the difierentiation of Vt
82
Vdd
Iin1 Iin2 Iout
Figure 3.16: Current subtractor.
results from the connection of the bodies of all NMOS transistors to Vss, and hence
the threshold voltages of M1, M2, M3 and M4 are no long equalized. Therefore,
we use PMOS devices to build the Tanno multiplier in our circuit to tie the bodies
and the sources of PMOS transistor to n-wells for diminishing the body efiect.
Compared to other current multipliers, the Tanno multiplier has relatively sim-
ple circuit topology. The Tanno multiplier has another advantage that its two
inputs are completely equivalent, i.e. its output current is identical to the result
obtained after switching Ix and Iy, ignoring the errors caused by transistor mis-
match. This is opposite to some multipliers, for example the Wiegerink Multiplier
[53], in which the product results difier after switching the input multiplicands
due to the non-symmetrical inputs caused by unavoidable high order efiects of the
transistors.
The working ranges of the Tanno multiplier are restricted by the prerequisite
of the circuit that all transistors work in the saturation region. Taking the circuit
branch of M7 and M8 as an example, when Ix ? 0, we have:
IM7 +Ix = IM8 > Ix > 0: (3.85)
83
The upper bound of the magnitude of the input current Ix is limited by the gate
voltage of M8, which needs to be lower than Vb to make M7 and M8 work in the
saturation region. Vb is given by Vb = 2Va = 2(
q I
B
mKM8 + Vt). Therefore IM8 can
be bounded by:
IM8 = KM8(Vgs;M8 ?Vt)2
< KM8(Vb ?Vt)2 (3.86)
= KM8[2(
s I
B
mKM8 +Vt)?Vt]
2
< KM8(2
s I
B
mKM8)
2
= 4IBm: (3.87)
Combining equation (3.87) with (3.85), we have:
4IBm > IM8 > Ix > 0
) Ix < 4mIB (Ix ? 0): (3.88)
Similarly, when Ix < 0, we can derive that:
?Ix < 4mIB (Ix < 0): (3.89)
Putting the equations for Ix ? 0 and Ix < 0 together, and applying them to the
circuit branches of inputs Ix+Iy and Iy, we derive the working range of the inputs
of the Tanno multiplier as described below:
8
>>>>
>><
>>>
>>>:
jIx +Iyj < 4mIB;
jIxj < 4mIB;
jIyj < 4mIB:
(3.90)
84
Q1
I1
Q2
I2
Q3
I3
Q4
Ic1
Q5
Ic2
Pool Circuit1
V+1V
+2V
-
Vout
Vin  VoutV BufferVb
Iout
Pool Circuit2 
V+1V
+2V
-
Vout
Va
Q6
Figure 3.17: Multiplier based on NPN transistors.
Multipliers Based on Bipolar Devices
The Tanno multiplier perform the multiplication of two inputs. In the case that
the product of several variables is to be calculated, the multiplier can be cascaded.
However a cascaded multiplier has a problem of accumulating errors, and due to
the peripheral circuits necessary to duplicate currents in the inputs and sometimes
to inverse the current directions in the input or output, the extra circuits may
introduce further inaccuracies. An operation like the one presented in the Beeler-
Reuter model, m3hj (refer to the equations (2.44)), requires the multiplication to
be carried out four times, which results in a big complicated Tanno type of multi-
plier circuit, and this motivates us to flnd a circuit for easy cascaded multiplication.
We observe that the input variables of m3hj are all positive from their physical
meanings, and this inspires us to take advantage of the exponential property of
NPN transistors for realizing the cascaded multiplication, and the implemented
multiplier is the flrst quadrant multiplier.
In the example presented below there are three multiplicands: I1, I2 and I3,
which are all positive. The circuit implementation is illustrated in Figure 3.17.
85
Q1?Q5 are connected as diodes and their base voltages are decided by I1 ?I3,
Ic1 and Ic2 in the logarithm relation described by:
Vbe;z = VT ln IzI
s0
(z is 1; 2; 3; c1 or c2; Iz > 0): (3.91)
Ic1 and Ic2 are constant bias currents. The pool circuits (refer to section 3.3.1
Figure 3.9) add I1, I2 and I3 together, and subtract Ic1 and Ic2 from the sum, i.e.
Vb is:
Vb = Va +Vbe3 ?Vbe;c2 = (Vbe1 +Vbe2 ?Vbe;c1)+Vbe3 ?Vbe;c2: (3.92)
Note that the inputs of the pool circuits are the gates of MOS transistors, and,
thus, the efiect of the pool circuit on Q1?Q5 can be treated as none. Therefore,
when Ic1 = Ic2, Q4 or Q5 can be saved with one NPN transistor driving two inputs
of the pool circuit due to the ignorable loading efiect from the pool circuit.
The voltage bufier (refer to section 3.3.1) copies Vb to the output, and provides
current to the base of the output transistor Q6. The current of Q6 collector is
exponentially decided by the base-emitter voltage, and can be expressed in terms
of the input currents from equation (3.92):
Iout = Is0e
Vb
VT = I1I2I3
Ic1Ic2 : (3.93)
The coe?cient difierence between Iout and the real product of I1 ? I2 ? I3 can be
easily flxed with parameter scaling and current mirrors.
For a general situation in which N +1 multiplicands are involved, the circuit in
Figure 3.17, specifying the case N = 2, can be extended by adding NPN transistors
andthepoolcircuits. ForanoperationwithN+1multiplicands, N+1NPNdevices
are needed to port the input currents I1 ? IN+1, N NPN transistors are used to
host the constant currents Ic1 ?IcN for canceling the Is0 terms in the equation of
86
the product result, and N pool circuits are needed to perform the addition and
subtraction. If there are identical currents among Ic1?IcN, some bias NPN devices
can saved. For example, in Figure 3.17, if Ic1 = Ic2, Q5 can be saved by connecting
the base of Q4 to the second pool circuit. The circuit represents a mathematical
equation expressed as follows:
Iout = I1I2I3:::IN+1I
c1Ic2:::IcN
: (3.94)
Like other cascaded circuits, this type of multiplier circuit also has the issue of
accumulated errors. The magnitude of the error largely depends on the accuracy
of the linear addition/subtraction operation of the pool circuits. Fortunately, for
designing circuits of the Beeler-Reuter model, the errors can be kept well under
control. In the pool circuit shown in Figure 3.10, the error mainly comes from the
channel length modulation efiect of M3 and M4. For the case that the swing of
the inputs is small, it is possible to reduce the difierence of Vds3 and Vds4 to a very
small range. Due to the restriction of the physical meanings in the heart cell model,
the input currents of the multiplier can never reach the value of absolute zero and
usually have a relative range (i.e. the maximal value divided by their minimal
value) within 105. Hence, the input voltages of the pool circuits are normally
between 0:4V and 0:7V, corresponding to a range of current changing between ?1
and ?105. Therefore, the  uctuation of the inputs are within 0:3V, which can be
treated small for reducing the channel length modulation efiect.
Shrinking Current with Bipolar-Based Circuit
Current mirrors are probably the most popular current amplifying/shrinking cir-
cuits. They work based on the linear dependence of the drain-source currents of
MOS transistors on the W:L ratio, or the linear dependence of the collector cur-
87
Qn1I
in
Qn2
In1
Qn3
In2
Pool Circuit1
V+1V
+2V
-
Vout Vin  Vout
V Buffer1
Qn4
Ip1 Ip2
Pool Circuit2
V+1V
+2V
-
Vout Vin  Vout
V Buffer2
Qp1 Qp2
Qp3
Qp4
Current Mirror1
Current Mirror2
Iout
Va
Io1
Io2
Figure 3.18: Current amplifler/shrinker based on bipolar transistors.
rents of NPN transistors on the size of their emitters. However, current mirrors
may take big layout areas. In the case that a current needs to be shrunk by 104
times, it is impractical to build a current mirror circuit composed of transistors
with one having a W:L ratio of 104 times of another. Then, cascaded current
mirrors are needed for this case. In the following, we propose a new approach to
implement current ampliflers/shrinkers, and it is based on bipolar transistors.
The circuit is depicted in Figure 3.18, in which, the dotted area is a multiplier
circuit like the one shown in Figure 3.17, except in this circuit only two multipli-
cands Iin and In1 are illustrated, and In1 is flxed. The lower half of the circuit,
based on PNP transistors, works just the same way as the upper half except work-
ing compensatively for the case Iin < 0. When Iin > 0, Iin goes through Qn1, which
decides the voltage Va, and Qp1 does not take any current. Hence the voltages on
the bases of the output transistors Qn4 and Qp4 are:
8>
><
>>:
Vbe;n4 = Va +Vbe;n2 ?Vbe;n3;
Vbe;p4 = Va +Vbe;p2 ?Vbe;p3:
(3.95)
If the selection of the currents In1, In2, Ip1 and Ip1 makes jVbe;n2 ? Vbe;n3j < jVaj
and jVbe;p2 ?Vbe;p3j < jVaj, the voltages on the bases of the output transistors Qn4
88
and Qp4 are then both positive in equation (3.95). Therefore Qn4 creates a non-
zero collector current, while Qp4 is in the cut-ofi region and conducts no current.
Similarly, when Iin < 0, Iin goes through Qp1 alternatively, and the output current
Iout is delivered by Qp4 while Qn4 is cut ofi. Therefore the collector currents in Qn4
and Qp4 are expressed by:
Ic;n4 =
8>
><
>>:
IinIn1
In2 (Iin ? 0);
0 (Iin < 0):
(3.96)
Ic;p4 =
8>
><
>>:
0 (Iin ? 0);
IinIp1
Ip2 (Iin < 0):
(3.97)
The current mirrors on the right reverse the direction of Ic;n4 and Ic;p4, and
send the currents to the output:
Iout = kIin; (3.98)
where the ampliflcation factor k is decided by k = In1In2 = Ip1Ip2.
The proposed current amplifler/shrinker circuit can be used to magnify ca-
pacitances in our VLSI heart circuits. Capacitors serve as integration devices in
our VLSI design, and their voltages are the outputs that drive the inputs of the
succeeding circuits. Since the inputs of the succeeding circuits are gates of MOS
transistors and take no currents, we can build the integration circuits with smaller
capacitors using the current shrinker circuits and still keep the waveforms of the
voltages across the capacitors. This will be introduced in section 3.4.7.
3.4 VLSI Design of Beeler-Reuter Model
In this section, we are going to present the VLSI circuit for implementing the
Beeler-Reuter model. The mathematical description of the model for VLSI realiza-
89
Iout   
IK1
Vin  
Iout   
Ix1
Vin  
Iout   
INa
Vin  
Iout   
Is
Vin  
Cm IextV
m?E 
Vm
-Vm0
Figure 3.19: Top-level circuit diagram of the Beeler-Reuter model.
tion was summarized in section 3.2.4. The circuit building blocks of mathematical
functions need for the cardiac cell model have all been introduced in section 3.3.
A top level block diagram for realizing the cardiac cell model is shown in Figure
3.19, and it is mapped from the equations that describe the membrane potential
Vm with the initial condition shifting process (refer to section 3.2.3), i.e. equations
(2.40) and (2.41), which are combined and rewritten as follows:
CmdVmdt = ?(IK1 +Ix1 +INa +Is)+Iext: (3.99)
In the diagram in Figure 3.19, the membrane capacitor Cm is charged by flve
current sources, an external current Iext which serves as the stimulus to activate
the cardiac cell model, and four ionic channel currents, IK1, Ix1, INa, and Is. E
is used to shift the resting initial voltage of the capacitor to zero, and Vm0 is
the original initial condition of Vm. Note that the positive input terminal of E,
connected with one side of the capacitor, does not take any current and, hence,
E does not contribute any current to charge the capacitor. The ionic current
modules are nonlinearly controlled by the membrane potential Vm and do not
draw any currents at their inputs.
Figure 3.20 provides more details for the four ionic currents. IK1, having the
simplest formulation among the four, is determined only by Vm, and thus can be
described by a voltage-controlled current source, as shown in Figure 3.20(a). Ix1
90
Vm Ix1Vm
IK1
G 
m3
Vm
gNaV
ENa
INa
gNaC
(b)
(a)
(c)
a
b
c
?
?
?
Vm
gs
Is
Is
Is
aG 
Es
F 
(d)
+
?
Vin      IoutI
K1
Vin      Ioutm
Vin      IoutI
x1
Iin      IoutE
s
Vin      Iouth
Vin      Ioutm
Vin      Ioutx1
Vin      Ioutd
Vin      Ioutf
Figure 3.20: Diagrams of ionic current modules: (a)IK1, (b)Ix1, (c)INa, (d)Is.
can be separated into two parts as stated in equation (2.43): the Vm dependent Ix1,
and a potential- and time-dependent gating variable x1. Hence, Ix1 can be modeled
as the product of Ix1 and x1, as illustrated in Figure 3.20(b). The diagram for INa
is shown in Figure 3.20(c), in which gating variable m is cubed and multiplied by
the gating variables h and j, and the product, amplifled by a factor of gNaV and
then added with the constant sodium conductance gNaC, is multiplied with the
output of the linear transconductance G whose output current is linearly decided
by Vm ? ENa; this implements equation (2.44). The diagram for Is is shown in
Figure 3.20(d). As described in equation (2.45), Is is the product of a constant
transmembrane conductance gs, two gating variables d and f, and the difierence
of Vm and the calcium equilibrium potential Es. The resulting Is is duplicated into
two copies in Figure 3.20(d), with one being the output to charge the membrane
capacitor Cm, and the other feeding back to control the generation of Es.
91
Vm? IK1?
I0
F 
Gain=1
Vdd
Vm? IK1?Vin      Iout
IK1
Vin       Iout
Sigmoid 1
Vin       Iout
Sigmoid 2
Figure 3.21: IK1 module is composed of a constant current source and two sigmoid-
function circuits. Refer to Appendix A, schematics page Sch-2 to Sch-4 for tran-
sistor circuits.
Next we will introduce how the ionic currents are implemented with transistor
circuits using the function blocks presented in section 3.3, and explain the structure
of each ionic current module based on Figure 3.20. In the rest of the presentation,
we use the notations with a single quotation mark to represent the variables as-
sociated with reformulated equations. The simulation results of the VLSI circuits
are also provided in the following sections.
3.4.1 Time-Independent Potassium Current IK1
The description of IK1 is given in equation (3.5), and is rewritten as follows after
the variables are scaled:
I0K1 = 2:742?10?6 + 3:632?10
?6
e?5:966(V0m?0:3277) +1 ?
5:879?10?6
e11:77(V0m+0:8823) +1: (3.100)
Following equation (3.100), the overall diagram for modeling IK1 can be depicted
in Figure 3.21, where the constant current source I0 represents the flrst term on
the right side of the equation, and the two sigmoid sub-circuits implement the two
terms containing exponentials in the equation. The sigmoid functions are realized
92
Vm [mV]
IK1
(Vm
)
Figure 3.22: Ideal and simulated curves: IK1 vs. Vm.
using the coupled-emitter circuit shown in Figure 3.12. The current to current
converter F is a current mirror and is used to switch the direction of the output
current of the top sigmoid block in Figure 3.21. The complete transistor circuits
for ionic current IK1 are provided in Appendix A, schematics page Sch-2 to Sch-4.
The simulation result of implementing equation (3.100) using transistor circuits
is shown in Figure 3.22. The x-axis is the membrane potential Vm, whose value
ranges from -90 mV to 50 mV during an action potential. The solid line depicts the
ideal IK1 vs. Vm curve, i.e. calculated directly from equation (2.42), and the dashed
line is the simulated IK1 which is scaled back to the original magnitude (the scaling
equation is provided in Table 3.2, IK1 = 106I0K1) in order to be compared with the
ideal values. We adopt the average signal-to-error ratio (SER) and R-square as the
metrics to evaluate how close the circuit realization is to the original equations.
The deflnition of the signal-to-error ratio follows the signal-to-noise ratio, which
describes the power ratio of a given transmitted signal to the background noise,
93
and it is formulated as the following:
SER = 1N
NX
i=1
10log E
2
i
(Ei ?Oi)2 (Unit: dB); (3.101)
where Ei are the expected values, i.e. the values calculated from the mathematical
model, and Oi are the observed value, i.e. the circuit simulation results. R-square
is a statistic that measures how similar the waveforms of data sets are to each other
[68], and is also called the square of the multiple correlation coe?cient. R-square
is deflned as:
R-square = 1?
PN
i=1(Ei ?Oi)
2
PN
i=1(Ei ? ?O)2
; (3.102)
where Ei and Oi have the same meaning of those in equation (3.101), and ?O is the
mean over the N observed values. R-square takes any value between 0 and 1, and
a value closer to 1 indicates higher similarity between two data sets. The average
SER for realizing the function of IK1 with VLSI circuits is 36.0 dB. The R-square
between the ideal data and the simulated data for IK1 is 0.9966.
3.4.2 Gating Variables
Gating variables are described by flrst-order difierential equations, which are given
in equation (2.46) and are rewritten in the following:
dz
dt = fiz ?(fiz +flz)z; (3.103)
where z is one of the gating variables, x1, m, h, j, d, or f. The six gating
variables of the Beeler-Reuter model can be implemented with a circuit of the
same structure, as illustrated in Figure 3.23, with difierent parameters. In the
diagram shown in Figure 3.23, fi Function and fl Function are two black boxes
that model the opening and closing rates of the gating variable, and their VLSI
94
Vin      Iout
a0Function
Vin      Iout
a1Function
Vm? Cz
G 
a
b
c
z?m0
m0
d0
--+
I1
I2
+
- Vz
Iout2
F ?
?
a2z?
a3z?
Gain=g
Gain=1
z?
z?
Vbias
Vin     Iout
z
Vm? z?
Figure 3.23: Detailed structure of gating variable circuits x1, m, h, j, d, and f.
circuits are all difierent as they depend on the equations of the opening and closing
rates (provided in equations (2.47)-(2.58)). The fi Function and the fl Function
are controlled by voltage V 0m, i.e. the scaled membrane potential Vm, and create
output currents fi0z and fl0z. The fi Function has two identical outputs, with one,
joined with fl0z at point a, going to the upper multiplier, and the other alone feeding
into the lower multiplier. The products of the multipliers, in the form of currents,
are subtracted and charge the capacitor Cz. The voltage of the capacitor Vz is
connected to a linear transconductor G, whose negative terminal is connected to a
constant voltage Vbias, which is used to shift the initial voltage of the capacitor Cz
as introduced in 3.2.3. The output current of G is replicated into two copies, and
one copy is the output which is sent to the succeeding circuits in the cardiac cell
system, and the other copy serves as another multiplicand of the upper multiplier.
Here we use z0, fi0z and fl0z to indicate the variables after being scaled, to distinguish
them from the original ones denoted with z, fiz and flz.
95
In Figure 3.23 the voltage across the capacitors is described by:
CzdVzdt = I2 ?I1: (3.104)
Using the relation of Vz and z0 and the description of the input/output of the
multipliers (refer to section 3.3.3 equation (3.84)) provided as follows:
z0 = (Vz ?Vbias)?g; (3.105)
I1 = (fi0z +fl0z)z0m0; (3.106)
I2 = fi0zd0m0; (3.107)
where g is the transconductance of G, d0 is a constant current, and m0 is a constant
factor of the multipliers. We therefore obtain the following equations:
Cz
g
dz0
dt = fi
0
zd0m0 ?(fi
0
z +fl
0
z)z
0m0;
) dz
0
dt = fi
0
z
d0m0g
Cz ?(fi
0
z +fl
0
z)
z0m0g
Cz : (3.108)
Suppose the scaling factor for fiz and flz is Afiz, and the scaling factor for z is
Az, we have:
8>
>>>>
><
>>>>
>>:
z = Azz0;
fiz = Afizfi0z;
flz = Afizfl0z:
(3.109)
Taking the above scaling equations into equation (3.103), we derive the following:
Azdz
0
dt = Afizfi
0
z ?Afiz(fi
0
z +fl
0
z)Azz
0;
) dz
0
dt =
Afiz
Az fi
0
z ?Afiz(fi
0
z +fl
0
z)z
0: (3.110)
96
Table 3.4: Vbias for initial conditions of gating variable capacitors (unit: V).
Gating variable x1 m h j d f
Vbias -0.0044 -0.00989 -0.4732 -0.468 -0.00188 -0.5
Comparing equation (3.108) with (3.110) yields the expression for the relation
between the scaling factors and the constant parameters of the circuit:
Az = 1d
0
; (3.111)
Afiz = m0gC
z
: (3.112)
The selection of the scaling factors and the circuit parameters are restricted by
the realization circuit of G (refer to section 3.3.1). The range of variation of Vz is
required to be small enough in order to achieve high linearity in the transconductor
G. In our VLSI design, the variation of the difierential input of G is constrained
to be within 0.5 V in the circuit. The constant current source d0 is selected to
be 5 ?A, and, hence, z is scaled down by a factor of 2 ? 105 (refer to Table 3.2
for the scaling of z ). m0 and g are set to 105 and 10?5 respectively. The scaling
factor of fiz and flz is determined by their original maximum values, which, after
being scaled, are restricted to be within the working ranges of the multipliers. The
maximum magnitude of fiz and flz are difierent for for difierent gating variables,
ranging from 10?3 to 102 , and this decides difierent values of the capacitance Cz
in the circuit. The scaling factors of fiz (and flz) are given in Table 3.2 for the
six gating variables in the Beeler-Reuter model, and the corresponding Cz will be
provided in section 3.4.7, discussed together with the time scale.
In Figure 3.23, the negative input terminal of G is connected to a non-zero
voltage Vbias to shift the initial voltage of Cz to be zero, as introduced in section
3.2.3. The magnitude of Vbias is decided by the original initial condition of the
97
gating variable, provided in Table 2.3, and the ratio of the capacitor voltage Vz
to the original gating variable z. Table 3.4 lists the values of Vbias for the gating
variables in the cardiac cell model, which make Vz(t = 0) = 0.
The transistor circuits of realizing the gating variables of the Beeler-Reuter
model are provided in Appendix A, and will be mentioned in the following sections
when presenting corresponding ionic currents.
3.4.3 Time-Dependent Potassium Current Ix1
Opening/Closing Rate of x1
The circuit for realizing gating variable x1 has been introduced in the previous
section and given in Figure 3.23. Here we only present the black boxes fi Function
and fl Function in the x1 circuit. The mathematical descriptions of fix1 and flx1
are provided in equations (2.47) and (2.48). They are rewritten as follows:
fi0x1 = 2:012?10
?3e8:3(V0m?0:5)
e5:7(V0m+0:5) +1 ; (3.113)
fl0x1 = 8:67?10
?5e?6(V0m+0:9)
e?4(V0m+0:2) +1 ; (3.114)
where V 0m, fi0x1 and fl0x1 are the membrane potential Vm, the opening rate fix1 and
the closing rate flx1 after being scaled. Here we reorganize the constant coe?cients
in the equations for easier circuit design. For instance, in the original equation
of fi0x1 (2.47), the coe?cient 0:0005 and the exponential with a constant power
e0:083?50 is re-organized as follows:
0:0005?e0:083?50 = 0:0005?e(+0:083?100)e(0:083?50?0:083?100) = 2:012?e(?0:083?50):
(3.115)
It is the term 2:012?10?3e(?8:3?0:5) in equation (3.113) after being scaled. Equa-
tions (3.113) and (3.114), very similar to exponential-sigmoid equation (3.70), can
98
Vm [mV]
????x1
(Vm
)
Vm [mV]
(a) (b)
????x1
(Vm
)
Figure 3.24: Ideal and simulated opening and closing rate of x1. (a)fix1, (b)flx1.
be implemented with the emitter-coupled pairs combined with a single exponential
circuit block introduced in section 3.3.2. The block diagram of the exponential-
sigmoid circuit is shown in Figure 3.12, with I0 being replaced with the output
current created by the circuit depicted in Figure 3.13. The transistor circuits for
fi0x1 and fl0x1 are given in Appendix A, schematics page Sch-6 and Sch-7. The
simulation results of the transistor circuits are shown in Figure 3.24(a), for fix1,
and 3.24(b), for flx1. In both flgures, the solid curves depict the ideal relations of
the opening/closing rates vs. Vm, calculated from equations (2.47) and (2.48), and
the dashed curves are the simulated results of the implementing circuits (schemat-
ics page Sch-6 and Sch-7 in Appendix A) after being scaled back to the original
magnitude for comparison (refer to Table 3.2 for the scaling of fix1 and flx1). The
average SER (signal-to-error ratio) for the circuit realization of fix1 is 15 dB, and
is 20.5 dB for flx1, and the corresponding R-square is 0.9907 for fix1 and 0.9978 for
flx1.
99
Vm [mV]
Ix1
(Vm
)
Figure 3.25: Ideal and simulated curves: Ix1 vs. Vm.
Maximum Current of Time-dependent Potassium Channel: Ix1
As depicted in Figure 3.20(b), Ix1 is the product of x1 and Ix1. The expression
for Ix1, given in equation (2.43), is transformed with the Matlab Curve Fitting
Toolbox for easier circuit realization, and is rewritten as follows after scaling the
variables:
Ix10 = 4:264?10?6 ? 50?10
?6
e4:345(V0m+1:31) +1; (3.116)
where we select the scale factor of Ix1 to be 106, i.e. Ix1 = 106 ? Ix10. Following
equation (3.116), Ix10 is realized in a VLSI circuit with a constant current source
and the sigmoid-function circuit illustrated in Figure 3.12. The transistor circuit
for Ix10 is in Appendix A, schematics page Sch-8. The PSpice simulated Ix10
vs. V 0m, after being scaled back, is shown in Figure 3.25 with the dashed line,
compared with the ideal Ix1 vs. Vm curve of equation (2.43) depicted with the
solid line in the same Figure. The average SER (signal-to-error ratio) for the
circuit implementation of Ix1 is 27.7 dB, and the R-square between the ideal and
the simulated data is 0.9886.
100
Mb1
Mb2
Mb3
M01
M02
M03
Mxy1
Mxy2
Mxy3
Mx1
Mx2
Mx3
My1
My2
My3
Vss
Vdd
Msub1
Msub2
Msub3
Msub4
Msub5
Msub7
Msub6
Msub8
Msub9
Msub10
Msub11
Msub12
Msub13
Msub14
IoutIB
Ix+ Iy Ix Iy
Current 
Subtractor
Figure 3.26: PMOS current multiplier based on Tanno multiplier.
Multiplier for Ix1 ?x1
Ix1 is multiplied with the gating variable x1 to generate the flnal potassium current
Ix1, as illustrated in the diagram in Figure 3.20(b), and described by equation
(2.43). After the scaling and initial value shifting processes (provided in Table
3.2), equation (2.43) is changed to the following:
I0x1 = 2?105Ix10x10: (3.117)
The multiplier is based on the Tanno circuit, as introduced in section 3.3.3, with
some modiflcation to adapt it to our application. The transistor circuit of the
multiplier is depicted in Figure 3.26, which is composed of a quarter squarer,
shown as the unshaded area in Figure 3.26, and a current subtractor, shown in the
shaded area in the flgure. Difierent upper-rail voltages are applied for these two
parts of the circuit because of the constraints from the usage of the inputs and
the output of the multiplier module. These constraints follow from the fact that
the input currents, coming from current mirrors of the previous stage, restrict the
101
Ix [?A]
Iou
t [
????A
]
Iy= 5?A
Iy=2.5?A
Iy=0?A
Iy=-2.5?A
Iy=-5?A
Figure 3.27: Simulation results of modifled Tanno multiplier.
voltages on the inputs to make the output transistors of the current mirrors work
in the saturation region. Since the voltages at the inputs are decided by the input
transistors Mxy1, Mx1, and My1, the sourcing PMOS transistors of the current
mirrors have very small working ranges in the saturation region if the previous
stage and the quarter squarer are powered by the same upper-rail supply voltages.
Hence, ground is employed instead of Vdd to serve as the upper rail voltage for
the quarter squarer circuit. For the output of the multiplier, the output voltage,
determined by the succeeding circuit, is limited to around 0 for our application, and
thus the current subtractor remains powered by Vdd. Difierent from the original
Tanno multiplier, PMOS transistors are used alternatively, and cascode structures
are employed in the bias block Mb2 and Mb3, and in the current subtractor to
mitigate the channel modulation efiect. The complete circuit of this multiplier
with transistor parameters is provided in Appendix A, schematics page Sch-9.
The PSpice simulation result of the modifled Tanno current multiplier is shown in
Figure 3.27.
102
Table 3.5: Equations of the opening and closing rate variables for INa.
Variable Reformulated Equation Sig-to-Err R-Square
Ratio (dB)
fi0m 1:097?10?5 +V 0m ?10?5 ? 6:292?10?6e?5:687(V0m+0:762089)+1 25.9 0.9989
fl0m 1:096?10?5e?0:05690(V0m+0:9) 17.0 0.9943
fi0h 1:625?10?5e?25(V0m+0:9) -54.5 0.8776
fl0h 8:5?10?6e?0:082(V0m+0:225)+1 39.8 0.9999
fi0j 5:53?10?5e?25(V0m+0:78)e?2(V0m+0:78)+1 -52.0 0.9972
fl0j 1:5?10?5e?10(V0m+0:32)+1 28.2 0.9995
3.4.4 Time-Dependent Inward Sodium Current INa
Opening/Closing Rates of Sodium Gating Variables
There are three gating variables for the fast sodium current, m, h, and j. Their
realization circuits, illustrated in Figure 3.23, have been presented in section 3.4.2.
Hence, we will only introduce the circuit implementing the black boxes for the
opening and closing rate variables next.
The original equations for the opening and closing rates are provided in equa-
tions (2.49)?(2.54), which are rewritten in Table 3.5 to take into account the refor-
mulation processes as described in section 3.2. The third and the fourth columns
of Table 3.5 list the SER (signal-to-error ratio) and R-square of the PSpice sim-
ulated results for the transistor circuit realization of the opening and closing rate
variables vs. the values given from the equations, and this will be discussed after
the simulation results are presented.
fi0m is realized following its expression shown in Table 3.5 by a constant current
source for the flrst term in the expression, a linear voltage to current converter for
103
the second term, and a sigmoid-function circuit for the third term. The difierential
circuit introduced in section 3.3.1, Figure 3.5, is adopted to realize the linear
voltage to current converter. The sigmoid-function circuit is shown in Figure 3.12
and presented in section 3.3.2. The detailed transistor circuit for implementing fim
is provided in Appendix A, schematics page Sch-11. The PSpice simulation results
of the circuit for fim compared with the ideal values calculated from the equation
are shown in Figure 3.28(a). fl0m is implemented with the single exponential circuit
depicted in Figure 3.13 according to the expression of fl0m given in Table 3.5, and
its transistor circuit is presented in Appendix A, schematics page Sch-12. The
circuit simulation results for flm are shown in Figure 3.28(b).
The expression for fi0h in Table 3.5 contains a single exponential term, and thus
can be modeled with the exponential circuit illustrated in Figure 3.13. The circuit
for realizing the sigmoid function for fl0h is depicted in Figure 3.12. The transis-
tor circuits for implementing fih and flh are provided in Appendix A, schematics
page Sch-13 and Sch-14, and their PSpice simulation results are shown in Figure
3.28(c) and (d) respectively with the dashed curves, compared with the ideal values
depicted with the solid curves.
The circuit for realizing fi0j uses the sigmoid circuit in Figure 3.12, but with
replacement of I0 with the output current of the single exponential circuit shown
in Figure 3.13, whose value is exponentially decided by the input voltage. The
equation of fl0j is realized, again, with the sigmoid-function circuit illustrated in
Figure 3.12. The transistor circuits for implementing fij and flj are provided
in Appendix A, schematics page Sch-15 and Sch-16. The curves of the circuit
simulated and ideal values of fij and flj vs. Vm are shown in Figure 3.28(e) and
(f).
104
Vm [mV]
????m
 (
Vm
)
????m
 (
Vm
)
Vm [mV]
(a) (b)
Vm [mV]
????h 
(Vm
)
Vm [mV]
(c) (d)
????h 
(Vm
)
Vm [mV]
????j 
(Vm
)
Vm [mV]
(e) (f)
????j 
(Vm
)
Figure 3.28: Ideal and simulated opening and closing rates of gating variables for
INa. (a)fim. (b)flm. (c)fih. (d)flh. (e)fij. (f)flj.
105
All circuits for the sodium opening and closing rates show high R-square values
as listed in Table 3.5, which indicate good correlation between the ideal and the
simulated results. Most circuits show good average SER, except fih and fij. The
low SER of fih and fij is because that, in a large range of the working region of
Vm, fih and fij are nearly zero, and this can be seen in Figure 3.28(c) and 3.28(e).
According to the deflnition of SER, i.e. equation (3.101), when the expected value
Ei tends to be zero and the observed value Oi is not as close to zero as Ei, the
SER then becomes negative and indicates a big difierence between the expected
results and the observed value. The ideal value of fih drops fast from around 3 at
Vm = ?90 mV to about 10?5 at Vm = ?50 mV and reaches the order of 10?15 with
increasing Vm. Similarly fij drops from around 0.09 at Vm = ?90 mV to the order
of 10?16 with increasing Vm. We observe that in the equation of gating variable
h (or j), that is equation (3.103), when fih (or fij) is close to zero, flh (or flj) has
high enough value to keep the term flz?z dominating the right side of the equation
(if z is not zero). Therefore we assume the close-to-zero values of fih and fij are
not as important as their larger values. Thus we ignore the big SERs of fih and fij
for the close-to-zero values. Recalculating the average SERs for fih in the range
Vm = [?90;?60] mV and for fij in the range Vm = [?90;?50] mV, where their
values are not so close to zero, the circuit implementation for fih and fij presents
a SER of 10.7 dB and 14.3 dB for each.
Sodium Current INa
The equation for implementing sodium current INa is from equation (2.44) and is
given as follows:
I0Na = 1:28?1023m03h0j0(V 0m ?0:5)+3?10?7(V 0m ?0:5); (3.118)
106
Vm?
Q1 Q
2 Q3
Q4
Ic1
Q5
Ic2 Vdd
m?3h?j?
Q6
F 
G1 G2 
?
Ena? INa?
Gain=g1 Gain=g2
Gain=f
NPN Based 
Multiplier
Ena? a
Vin  Iout
m?
Vin  Iout
h?
Vin  Iout
j?
Vin Vout
V BufferMulti-
Input
Pool 
Circuit
Vm?Vh?
Vj?
Vout
Vc1V
c2
Figure 3.29: Circuit diagram of INa. Refer to Appendix A, schematics page Sch-10
for transistor circuits.
where I0Na is the scaled sodium current, V 0m is the scaled membrane potential,
and m0, h0, and j0 are the scaled gating variables. The scaling factors are listed
in Table 3.2. Following this equation, Figure 3.29 provides more details of the
implementation circuit for INa. The three modules on the left top corner model
the gating variables m0, h0, and j0 and can be realized with the circuit structure
introduced in Figure 3.23. G1 and G2 are linear transconductors, with gains
of g1 = 5? and g2 = 0:3?, respectively, and realize V 0m ? E0Na with difierent
ampliflcation factors, where E0Na is, according to equation (3.118), 0.5 V. The
shaded area is a NPN-based multi-input multiplier that implements m03h0j0. The
output of the NPN-based multiplier is then multiplied with the output current
created by G1, and the product is added to the output of G2 to achieve the flnal
I0Na.
In the NPN-based multiplier depicted in Figure 3.29, Q1, Q2, Q3 take the
output currents from the gating variable modules and convert them into voltages.
107
Vb1
Vb2
Vdd
Vss
Vout
Vm?
Vm?
Vm?
Vj?
Vh?Vc1
Vc1
Vc1
Vc2
Diff 1 Diff 2
I1 I2
I3 Diff 3
I4 Diff 4
I5 Diff 5
Multi-
Input
Pool 
Circuit
Vm?Vh?
Vj?
Vout
Vc1V
c2
Figure 3.30: Circuit of a multi-input adder/subtractor for INa module. Refer to
Appendix A, schematics page Sch-17 for transistor circuits.
Q4 and Q5 provide biasing voltages. V-Bufier and Q6 transfer voltage back to
current. F serves to change the direction of the output current from the multiplier
and has a gain f = 1. Difierent from Figure 3.17, which uses cascaded pool circuits
as the adder/subtractor (refer to section 3.3.1), this NPN-based multiplier for INa
is based on the pool circuit shown in Figure 3.9 which is expanded for multiple
inputs as illustrated in Figure 3.30. Difi 1 corresponds to the difierential circuit
on the left side of the pool circuit in Figure 3.9. Difi 2 to Difi 5, with each being
equivalent to the difierential circuit on the right side of the pool circuit in Figure
3.9, are identical and are all connected to the output point Vout. Vout is connected
to the input of V-Bufier as shown in Figure 3.29, which does not draw any current
(refer to Figure 3.7(b) for V-Bufier). Hence, we have I1 = ?(I2 + I3 + I4 + I5).
The expression of I1 follows equation (3.44), and Ii (i = 2;3;4;5) follow equation
(3.45). If the non-linear terms (V1 ?Vout)2 and (V2 ?V3)2 in equations (3.44) and
108
(3.45) are ignored, the currents can be linearly added to yield the output expression
for Figure 3.30:
8>
>>>>
>>>
>><
>>>>
>>>
>>>:
I1 = Kn(V1 ?Vout)
q2I
B1
Kn ;
I2;D2 +I2;D3 +I2;D4 +I2;D5 = Kn(V2;D2 ?V3;D2 +V2;D3
?V3;D3 +V2;D4 ?V2;D4 +V2;D5 ?V3;D5)
q2I
B2
Kn ;
I1 = ?(I2;D2 +I2;D3 +I2;D4 +I2;D5);
(3.119)
) Vout = V1 +(V2;D2 ?V3;D2 +V2;D3 ?V3;D3 +V2;D4
? V2;D4 +V2;D5 ?V3;D5); (3.120)
where, in order to make the symbols agree with those in equation (3.45), we use
I2;Di (i = 2;3;4;5) to represent current Ii from Difi i, and use V2;Di and V3;Di
(i = 2;3;4;5) to represent the voltages of the positive and negative terminals of
Difi i. Following equation (3.120), the output Vout in Figure 3.30 is described by:
Vout = 3V 0m +V 0h +V 0j ?(3Vc1 +Vc2); (3.121)
where Vm0, Vh0, Vj0, Vc1, and Vc2 are the emitter-base voltages of Q1?Q5 in Figure
3.29. The output of the NPN-based multiplier is then expressed by:
IQ6 = m
03h0j0
(Ic1)3 ?Ic2: (3.122)
The detailed transistor circuits of the INa diagram in Figure 3.29 is provided
in Appendix A, schematics page Sch-10.
3.4.5 Time-Dependent Slow Inward Calcium Current Is
Opening/Closing Rates of Calcium Gating Variables
There are two gating variables for the calcium current Is: d and f. They are
realized by the gating variable circuit provided in Figure 3.23. The equations for
109
the opening and closing rates of gating variable d are rewritten from equations
(2.54) and (2.55) as follows:
fi0d = 2:456?10
?5e?(V0m+0:9)
e?7:2(V0m?0:05) +1 ; (3.123)
fl0d = 1:48?10
?5e?1:7(V0m+0:9)
e5(V0m+0:44) +1 ; (3.124)
where, again, we use single quotations to distinguish the scaled variables from
the original ones where the scaling factors are provided in Table 3.2. Here we
re-organize the constant coe?cients of the equations for easier circuit design. For
example, 0:095e[?0:01?(?5)] in the original equation (2.54) is converted as follows:
0:095e[?0:01?(?5)] = 0:095?e+(0:01?95) ?e(0:01?5?0:01?95) = 2:456?10?5e(?0:01?90);
(3.125)
anditisrewrittenas2:456?10?5e(?0:9) inequation(3.123). Similarly, theequations
of fif and flf are rewritten by re-organizing the constant coe?cients as follows from
equation (2.56) and (2.57):
fi0f = 1:97?10
?5e?0:8(V0m+0:9)
e15(V0m+0:28) +1 ; (3.126)
fl0f = 2:16?10
?5e?2(V0m+0:9)
e?20(V0m+0:3) +1 : (3.127)
Equations (3.123) to (3.127) are all (3.70)-type of equations, and, thus, can
be implemented with the sigmoid-function circuit shown in Figure 3.12, with I0
being substituted by the exponential circuit provided in Figure 3.13. The tran-
sistor circuits for realizing equations (3.123)-(3.127)are provided in Appendix A,
schematics page Sch-20 to Sch-23. Figure 3.31(a)(b) shows the simulation results
of the implementation circuits for fid and fld and also the curves of the ideal cases,
calculated directly from equations (2.54) and (2.55). The SER (signal-to-error ra-
tio) of the circuit realization for fid is 37 dB, and 26.8 dB for fld, and the R-square
110
Vm [mV]
????d 
(Vm
)
????d 
(Vm
)
Vm [mV]
(a) (b)
Vm [mV]
????f 
(Vm
)
Vm [mV]
(c) (d)
????f 
(Vm
)
Figure 3.31: Ideal and simulated opening and closing rate variables for Is. (a)fid.
(b)fld. (c)fif. (d)flf.
111
G4 
F1 
Gain=f1
F2 
Gain=f1
G1 
G2 
Ic1
Q2
G3 
NPN-basedMultiplier
d?f ?
Vm?-Es?
Iout
Q1
Vm?
Is?
a Is?
Vc
Gain=g1
Gain=g2
Vb
Es Block
b
Gain=g4
Gain=g3
Es?
Ic2
[Ca]?
C
Vin  Ioutd?
Vin  Ioutf ?
Log
Figure 3.32: Circuit diagram of Is. Refer to Appendix A, schematics page Sch-19
for transistor circuits.
for fid and fld are 1.0000 and 0.9933, respectively. The simulated results and the
ideal values for fif and flf are shown in Figure 3.31(c)(d). The SER of fif is 28.7
dB, and the SER of flf is 22.4 dB. The R-square for fif and flf are 0.9998 and
0.9990, respectively.
Calcium Current Is
Equations (3.128) to (3.130) provide the reformulated equations for the calcium
current Is after the scaling and initial condition shifting processes (refer to section
3.2.2 and 3.2.3, and Table 3.2) are performed on the original equations of the
Beeler-Reuter model (2.45), (2.59), and (2.60):
I0s = 3:6?105d0f0(V 0m ?E0s ?104); (3.128)
E0s = ?8:23?10?7 ?1:30287?10?7 ln[Ca]0i; (3.129)
d[Ca]0i
dt = ?0:1I
0
s +0:07(10
?7 ?[Ca]0
i): (3.130)
The diagram for implementing Is is illustrated in Figure 3.32. The d0 and f0
modules on the left consist of the gating variable circuit depicted in Figure 3.23.
112
The voltage V 0m is linearly converted by the same scale factor as E0s by G4, subtracts
E0s at node b, and then feeds the current into an NPN-based multiplier. The NPN-
based multiplier follows the same design strategy presented in section 3.3.3, and a
similar multiplier has been introduced as when we discussed the multiplier circuit
for INa. Therefore, we do not repeat the details for the NPN-based multiplier
for Is here; its transistor circuit can be found in Appendix A, schematics page
Sch-26. The calcium equilibrium potential E0s is created at node b by the circuit
illustrated in the shaded area of Figure 3.32, i.e. the Es Block. I0s is generated by
a NPN-based multiplier, which takes the scaled d0, f0 and V 0m ?E0s as the inputs,
and delivers its output current, i.e. I0s, to F1 and F2. F1 and F2 duplicate I0s and
sends one copy via F2 to charge the membrane capacitor in the succeeding circuit
and the other via F1 as a feedback to the Es block.
In the Es Block in Figure 3.32, a capacitor C and a linear transconductor G1
with a transconductance g1 construct the derivative function for equation (3.130).
The voltage across the capacitor Vc, is described by the following equation:
CdVcdt = ?I0s ?Vc ?g1; (3.131)
where Vc has the relation with [Ca0] that can be expressed by the following equation
using g2, the gain of G2:
Vc = [Ca]
0
i
g2 +Vb; (3.132)
where Vb is the bias voltage on the negative terminal of G2. Taking equation
(3.132) into (3.131) creates:
C
g2
d[Ca]0i
dt = ?I
0
s +g1 ?(?Vb ?
[Ca]0i
g2 ): (3.133)
Equation (3.133) can be equivalent with equation (3.130) by selecting proper values
of the circuit parameters. In our VLSI design, g1 is set to 7?, g2 is set to 10?, Vb is
113
-0.01 V, and C is 0.1 ?F. The transistor circuit for implementing equation (3.130)
is shown in Appendix A, schematics page Sch-24. [Ca]0i, generated by G2, is sent
to a logarithm circuit, shown in the oval in Figure 3.32, which implements equation
(3.129). The same as the logarithm circuit introduced in Figure 3.14, the logarithm
circuit for Es is composed of two NPN transistors Q1 and Q2, a constant current
sources Ic1, and a linear transconductor G3 with a gain g3 = 5?. The logarithm
circuit together with Ic2 realize equation (3.129). The two bias currents Ic1 and
Ic2 are set to 10?A and 0:674?A, respectively. The detailed transistor circuit for
creating Es is also provided in page Sch-24 in Appendix A.
3.4.6 Action Potential Vm
Puttingthecircuitsoftheioniccurrentspresentedintheprevioussectionstogether,
the VLSI circuit that imitates the behavior of a membrane potential is established
as shown in the block diagram in Figure 3.19 and in the schematics page Sch-
1 in appendix A. More than three thousand transistors are used to build the
full circuit for modeling the electrical activity of a single cardiac cell. When the
external stimulating source sends a current impulse of 20 ?A for a duration of
about 0.003 second, the circuit creates an action potential across the capacitor Cm
and the waveform is shown with the dashed line in Figure 3.33. The simulation
results, after being scaled back to the original magnitude, indicate that the circuit-
generated Vm has a resting potential of about -84.7 mV, a peak of the upstroke
of 47.7 mV, and a re-polarization duration (measured at 90% of re-polarization
completion) of 300 ms. The results are compared with the ideal Vm which is
calculated from the original Beeler-Reuter equations and depicted with the solid
line in Figure 3.33. At the same amount of stimulus (compared with the de-scaled
114
Time [ms]
Vm
[mV
]
20 ?A
3 ms
Stimulus Impulse
Figure 3.33: Ideal and simulated action potential.
stimulus), the ideal action potential has a resting potential of -83.0 mV, the peak
of the upstroke of 41.7 mV, and a re-polarization duration of 287 ms.
The waveforms of the ionic currents IK1, Ix1, INa, and Is during an action
potential are provided in Figure 3.34, with solid curves representing the ideal
ionic currents and the dashed curves depicting the simulated values of the VLSI
circuits. IK1 and Ix1 are outward currents and responsible for discharging the
membrane capacitance. Their maximum values are about 5.0 ?A=cm2 and 1.5
?A=cm2 respectively. INa andIs areinwardcurrents(beingnegative), whichcharge
the membrane capacitance. The magnitude of Is is on the same order as IK1 and
Ix1, which is about 4.3 ?A=cm2. INa, responsible for creating the steep upstroke at
the start of an action potential, has a very big magnitude compared to the other
ionic currents, which is around 140 ?A=cm2. INa is almost zero most of the time,
except in the spike in the upstroke phase of an action potential. This is because
most of the time, either the activation gating variable m or the inactivation gating
variables h and j are closed, and only in the upstroke phase when m transits
to open and h and j start to close, there is a narrow time period that all the
115
Time [ms]
(a)
Time [ms]
(b)
Time [ms]
(c)
Time [ms]
(d)
IK1
 [
????A
/cm
2 ]
Ix1 
[????A
/cm
2 ]
INa 
[????A
/cm
2 ]
Is 
[????A
/cm
2 ]
Figure 3.34: Ideal and simulated membrane ionic currents. (a)IK1. (b)Ix1. (c)INa.
(d)Is.
116
gating variables are open to allow the passage of the sodium ions. This property
of sodium channel decides the un-excitability of a cardiac cell for another action
potential when it is still in the flrst action potential and the gating variables m, h
and j have not  ipped back to the state of a resting cell. This will be discussed
again in chapter 4, when reentry and refractory period are introduced.
3.4.7 Discussions
Experimental Data of Beeler-Reuter Model
We have evaluated how close the VLSI realization is to the original equations in the
Beeler-Reuter model using the average signal-to-error ratio (SER) and R-square for
Vm-dependent equations. In this section, we will discuss the correlation between
the experimental data from which the cardiac cell model is extracted and the
equations of the model, and compare the correlation to that between the equation
and the PSpice simulated results, as illustrated in Figure 3.35. The comparison of
the data and the equations is di?cult and sometimes not possible due to missing
details in literature. The mathematical equations of the model are mostly based
on multiple results of experiments carried out by difierent researchers, and can
be adjusted from the original data with re-scaling or shifting along axis by Beeler
and Reuter to simulate the ventricle, whereas the information of how the difierent
data is combined to obtain the formulas and how the equations are re-scaled and
shifted are not mentioned in [15]. In addition, most experimental data is provided
by graphs rather than numbers, and this brings more di?culties to make numerical
comparison of the data and the models. Hence, we do not provide full comparison
of the data and the Beeler-Reuter model, and only present two examples to roughly
show how close the resulted equations in the model are to the original experimental
117
Experimental
Data
Extracted
Equations
Circuit 
simulation
results
How Close How Close
Compare
Equations are extracted
from experimental data
Circuits are made to 
calculate equations
Figure 3.35: Comparison among experimental data, extracted equations, and cir-
cuit simulation results.
data, and also how close the equations are to simulation of the presented VLSI
design.
In the flrst example, we discuss the formulation of IK1. The formula of IK1 in
the Beeler-Reuter model is converted from the equation used by McAllister etc. in
[14], whose parameters were selected to give a current vs. voltage diagram similar
to the experimental data published in [69]. The data points, totally 15, are shown
graphically in [69], and we measured them with a ruler for their values and re-draw
the data in Figure 3.36. The flnal equation of IK1 in the Beeler-Reuter model has
been re-scaled and shifted from McAllister?s equation and this information is not
published, and, thus, we scale and shift the equation of IK1 (2.42) by the following
formula to let the resulted equation give a best flt to the experimental data:
gIK1 = a1FK1(a3Vm +a4)+a2; (3.134)
where gIK1 is the potassium current after being re-scaled and shifted for matching
the data points, function FK1(Vm) is the right side of equation (2.42) for IK1, and
a1, a2, a3, and a4 are constants for scaling and shifting the X and Y axis. a1,
a2, a3, and a4 are decided to let equation (3.134) best flt the experimental data,
and the resulting gIK1 curve is illustrated in Figure 3.36. Compared to Figure 3.22,
118
Vm [mV]
IK1 [10-7A]~
Figure 3.36: Experimental data [14] (dots) and its fltting curve of steady-state
current-voltage relation of outward current in purkinje flbers.
which shows the IK1 curve calculated from the equation and the PSpice simulated
IK1 curve, the data points in Figure 3.36 are more scattered around the equation
curve than the simulated results. The average signal-to-error ratio (SER) of the
experimental data and the curve of equation (3.134) in Figure 3.36 is 23.1 dB, and
the R-square is 0.9490, showing worse correlation compared to the SER and the
R-square of the Beeler-Reuter equation calculated curve and the PSpice simulated
curve of Figure 3.22, 36.0 dB and 0.9966.
In the second example, we discuss the experimental data for obtaining the
equations of the gate opening/closing rates (fid and fld) of the gating variable d
(equations (2.55) and (2.56)). The equations of fid and fld are decided by the
measurement results of the time constant ?d and the steady value d1 at difierent
command voltage steps (refer to chapter 2 section 2.3.2). The relation between fid
119
Vm[mV]
Rat
e co
nst
ant
s o
f de
cay
 1/?
d(
sec
-1 )
Vm[mV]
Ste
ady
-sta
te a
tiva
tion
d  8
(a)
(b)
Figure 3.37: [48] (a) Steady-state activation d1 of calcium current Is; data ob-
tained from flve cow ventricular trabeculae superfused with solution containing
0.2(?), 0.45(?) or 1.8 mM CaCl2 (? ? 24). (b) Rate constant for decay 1=?d,
measured with solution containing 0.2(?), 0.45(?) or 1.8 mM CaCl2 (?).
and fld, and ?d and d1 is formulated as, rewritten from equations (2.9) and (2.10):
?d = 1=(fid +fld); (3.135)
d1 = fid=(fid +fld): (3.136)
The ?d and d1 curves used for deflning the equations of fid and fld in the Beeler-
Reuter model represent a collection of experimental results published in [47], [48],
and [70]-[73]. Here, in Figure 3.37 [48], we show one set of the data measured
by Reuter and Scholz. Figure 3.37 depicts ??1d in (a) and normalized d1 in (b)
in a semi-logarithm scale. Difierent symbols represent difierent samples in three
solutions with difierent calcium concentration. We get the equation data of ?d
and d1 from equations (3.135) and (3.136) from the mathematical description of
fid and fld in the Beeler-Reuter model, and then compare the calculated values
to the experimental data. Again, the values of the experimental data is obtained
120
by measuring the data points in the graphs with a ruler. The SER of the experi-
mental data and calculated values is 13.5 dB for d1 and 13.2 dB for ?d, and the
corresponding R-square is 0.8800 and 0.7115. The equation calculated ?d and d1
are also compared to those obtained from the PSpice simulation. Again equations
(3.135) and (3.136) are used to achieve simulated data of ?d and d1 from the sim-
ulated fid and fld values (simulation curves are shown in Figure 3.31(a)(b)). The
SER of the simulated data and the calculated values is 40.4 dB for d1 and 38.2
dB for ?d, and the corresponding R-square is 0.9999 and 0.9986.
In the above two examples, we showed the experimental data for IK1, ?d and
d1, and numerically estimated how close the data is to the extracted mathematical
model by computing the average signal-to-error ratio and the R-square coe?cient.
The correlation metrics were compared with those derived from evaluating the
similarity between the model and the PSpice simulated results. The comparison
shows that the transistor circuits? simulated data has higher correlation to the
model than the experimental data. Hence, we conclude the accuracy of the VLSI
design of the Beeler-Reuter model is satisfactory.
Time Scale and Capacitors
There are eight capacitors in the VLSI cardiac cell model: one is for the membrane
capacitance, one is for the difierential equation of the calcium concentration, and
the rest are for the gating variables x1, m, h, j, d, and f. The magnitude of
the capacitance is decided by the original equations in the Beeler-Reuter model
and also variable scaling factors (when ignoring time scaling). Take the membrane
capacitance Cm as an example, it is deflned in equation (2.40) and related to Vm
and the ionic currents. Since Vm is numerically shrunk by a factor of 100 (ignoring
121
Table 3.6: Capacitances in the VLSI realization of Beeler-Reuter model (unit: F).
Schematics Vm Ix1 INa Is
Cap Name Cm Cx1 Cm Ch Cj Cd Cf C[Ca]i
Cap in F 10?7 10?6 10?10 5?10?9 5?10?8 10?6 10?7 10?7
the unit) in the VLSI realization and Iion is reduced by a factor of 106 (refer to
Table 3.2), according to equation (2.40), Cm needs to be 104 times smaller than
the original needed value 1 F in order to balance the right and the left sides of the
equation, and, thus, in the circuit, the capacitance of Cm is 0.1 mF. Note that the
units of variables employed by the cardiac cell model do not in uence the variable
scaling and are ignored, and only the values of the variables matter.
Time can be scaled easily by increasing or decreasing the capacitors with a
same factor. When the capacitors are enlarged, it takes longer time to charge
the capacitors, and, thus, the waveforms of the variables on the time axis are ex-
panded. Vise versa, reduced capacitances spend less time to be charged and make
the waveforms narrower. The Beeler-Reuter model uses ms as its time unit, and
an action potential lasts for about 300 time units (300 ms). Therefore, when the
model is mapped into the circuit, which is based on the second as the time unit, the
duration of an action potential becomes 300 s. In order to make the circuit? gen-
erated waveforms consistent with those generated from the Beeler-Reuter model,
the capacitance of the eight capacitors for implementing the Beeler-Reuter model
in VLSI are shrunk by 1000 times. All the simulation results shown previously are
with the scaled capacitors. The values of the capacitors are summarized in Table
3.6.
The VLSI circuits of the Beeler-Reuter model presented previously are able to
122
be realized on-chip, except for the capacitors. As shown in Table 3.6, the sizes of
the capacitors range from 100 pF to 1 ?F. A capacitance as big as 1 ?F is hard to
be realized on a chip. In the following we will discuss the methods for implement-
ing big capacitors with smaller capacitors to solve this issue. Capacitor multipliers
have been proposed for increasing the capacitances of compensation capacitors in
[74]-[76]. However their schemes are not applicable to our circuit design of the
heart model as next discussed. The compensation capacitors are used to introduce
poles and are small, and, thus, the capacitor multiplication factors which are less
than hundreds are su?cient to make on-chip compensation capacitors. In the VLSI
design of the heart model, as provided in Table 3.6, the capacitance can be as high
as 1 ?F, and this requires a multiplication factor of at least 104. Difierent from the
compensation capacitors, the capacitors in the heart model are used as integration
devices, and, hence, demand higher accuracy and also bi-direction. The additional
requirement of accuracy and bi-direction brings more di?culties in using the pub-
lished capacitor multiplication techniques. Capacitor multiplication schemes can
be divided into two categories: voltage-mode and current-mode. Voltage-mode,
working the same way as amplifying Miller capacitance, can not be used in our
case, because it amplifles the voltage signal with the same factor as the capacitor,
which makes a voltage signal too high to be supplied by the power sources. The
current-mode method amplifles capacitance by decreasing the capacitor currents.
For example, if a capacitor is 1 ?F and its maximum input current is 1 ?A, the
capacitor can be substituted with a capacitor 106 smaller by shrinking the current
by a factor of 106, i.e. the maximum input current becomes 1 pA, with the time
scale unchanged. A di?culty of implementing this method is the current shrinker.
Implemented with current mirrors, a current reduction factor of 106 means not only
123
C
k
1
k
Iin I
c=Iin/ k
VcC
Iin V
c Succeeding
Circuit 
Succeeding
Circuit 
Figure 3.38: Current amplifler/shrinker circuit used for amplifying capacitance.
transistors of very big sizes in a cascaded circuit topology, but also tremendous
challenge in designing circuits with accuracies to pA because high-order efiects of
transistors can easily overwhelm the tiny currents.
Here we propose a capacitor ampliflcation circuit that can be used in the VLSI
heart model. The circuit is based on current-mode, and has been mentioned in sec-
tion 3.3.3 when we introduce the bipolar-based shrinking-current circuit provided
in Figure 3.18. A diagram is given in Figure 3.38 to show how a current shrinker
can be used to magnify the capacitance in our application. On the left side of the
flgure, the capacitor C serves as an integration device and its voltage, simulating
a derivative variable in the heart model, drives a succeeding circuit which does
not draw any current from the capacitor (connecting to the gates of MOS transis-
tors). This circuit can be implemented with a capacitor k times smaller by using
a current shrinker of factor 1k (refer to Figure 3.18 for the current shrinker circuit),
illustrated on the right side of Figure 3.38, while the capacitor voltage Vc remains
the same for driving the succeeding circuit.
In the following example, we show the simulation results of the presented circuit
for magnifying the capacitance for 104 times. Since the present educational VLSI
MOSIS runs, using AMI 1.5 ?m ABN technology, do not support PNP transistors,
the PNP devices in the current shrinker shown in Figure 3.18 are replaced with
NPNs. The circuit is modifled accordingly and redrawn in Figure 3.39, in which Qpi
124
Qn1I
in
Qn2
In1
Qn3
In2
Pool Circuit1
V+1V
+2V
-
Vout Vin  Vout
V Buffer1
Qn4
Ip1 Ip2
Pool Circuit2
V+1V
+2V
-
Vout
Vin  VoutV Buffer2Qp1 Q
p2 Q
p3
Current Mirror1
Iout
Va
Io1
Io2
E 
Vbb
Vb2
Vb1
Vb
Vc
Qp4
Figure 3.39: Current amplifler/shrinker circuit using NPN transistors.
(i = 1;2;3;4) are the NPN transistors and they keep the names of PNP transistors
used in Figure 3.18. A voltage-to-voltage converter E is added to convert the
output voltage of Pool Circuit2 Vb to Vc, and Vc = Vb1 + Vb2 ?Vb, where Vb1 and
Vb2 are constant voltages. Note that in order to use NPN transistors, E needs to
be added, which introduces more errors. Hence, if PNP is available, the original
circuit shown in Figure 3.18 is preferred. The circuit in Figure 3.39 works the same
way as the circuit in Figure 3.18. When Iin > 0, Qn1 and Qn4 are activated, Qp1
and Qp4 are cut ofi, and the output current Iout is supplied by Qn4. When Iin < 0,
Qp1 and Qp4 are activated, Qn1 and Qn4 are cut ofi, and the output current Iout is
supplied by Qp4.
We set the circuit parameters to let Iout = 10?4Iin, and build the circuits
like the ones shown in Figure 3.38 without \Succeeding circuit". The original
capacitanceisC = 1?F, and with thecurrentshrinker, theintegrationfunction can
be implemented with a capacitor 104 smaller, i.e. C=k = 100pF. The waveforms
of the capacitor voltage are shown in Figure 3.40(a), with the solid line being the
voltage across capacitor C, and the dashed line being the voltage across capacitor
125
Time [s] Time [s]
Cur
ren
t [nA
]
Vol
tag
e [V
]
(a) (b)
Figure 3.40: Simulation results of using current shrinker circuit to amplify capac-
itance. (a)Voltages of original capacitor and amplifled capacitor. (b)Currents of
original capacitor and amplifled capacitor.
C=k. The currents going into the capacitors are shown in Figure 3.40(b), where
the solid curve is the current of capacitor C (i.e. Iin) minimized by a factor of
104 in order to be compared with the current in C=k, and the dashed curve is the
current of capacitor C=k. Here we adopt the current running through the capacitor
of the gating variable x1 during an action potential as the input current Iin. The
difierence between the ideal Iout (i.e. Iin=104) and simulated Iout is caused by
error sources such as the Early efiect of the NPN transistors, and the non-idealism
of the voltage adders, i.e. the pool circuits, voltage bufiers, and the voltage-to-
voltage converter E. More accurate voltage adders can be used to reduce the error
and improve the linearity between Iin and Iout. Figure 3.41 shows a scheme to
mitigate the Early efiect by using the cascode structure. The flgure depicts only
the upper half of the circuit in Figure 3.39 for improvement, and the lower half
can be modifled the same way. The shaded areas are newly added. Qn5 and Qn6,
together with Qn1 and Qn4 incorporate the cascode structure, which reduces the
126
Qn1Iin
Qn2
In1
Qn3
In2
Pool 
Circuit1
V+1V
+2V
-
Vout Vin  Vout
V Buffer1
Qn4
Current 
Mirror1
Iout
Io1
Pool 
Circuit3
V+1V
+2V
-
Vout Vin  Vout
V Buffer3
Qn6
Qn5
Figure 3.41: Cascode structure can be used to minimize the Early efiect on Qn4
and improve the accuracy of Iout.
impact of the Early efiect on the output current.
The presented capacitor ampliflcation circuit is limited by the leakage current
of the bipolar transistors. When a current is shrunk to be too small, the leakage
current dominates the overall output current of the current shrinker and breaks
the linear relationship between Iin and Iout. Hence, bipolar devices with low leak-
age currents are desirable for implementing the proposed capacitor ampliflcation
method.
Note it is possible to implement the circuit in Figure 3.18 with all NPN and
PNP transistors being replaced by NMOS and PMOS transistors that work in the
subthreshold region. The current in a subthreshold MOS transistor is given by
[77]:
IDS = I0e(kVGB=VT)(e(?VSB=VT) ?e(?VDB=VT)) (3.137)
= I0e(kVGB=VT)[e(?VSB=VT)(1?e(?VDS=VT))]; (3.138)
where I0 is the subthreshold current-scaling parameter, k is the subthreshold ex-
ponential coe?cient, and VT is the thermal voltage. When VDS > 5VT, the term
127
e(?VDS=VT) can be ignored, and the above equation can be rewritten to:
IDS = I0e(?VSB=VT)e(kVGB=VT) (3.139)
= I0e(kVSB?VSB)=VTe(kVGS=VT): (3.140)
When VSB is a constant, the current IDS is only exponentially dependent on the
gate-source voltage VGS, which makes the I-V properties of subthreshold MOS
devices very similar to bipolar transistors. Therefore, the capacitor ampliflcation
circuit in Figure 3.18 can work the same way when its bipolar devices are replaced
by MOS transistors by connecting the gate, source and drain of MOS devices
where the base, emitter and collector of the NPN devices originally are connected
respectively, and connecting the body to Vss. With MOS transistors, the currents
of the input and output terminals, and the bias currents In1, In2, Ip1, and Ip2 in
the circuit in Figure 3.18 have to be restricted within proper ranges to necessitate
the subthreshold conditions of the replacing MOS transistors.
In the VLSI design of the cardiac cell model, the integration capacitors can be
further minimized by shrinking the time scale. The current VLSI circuit is able to
function properly as a cardiac cell with a time scale shrunk by up to 1000 times
of the current time scale, that is, the generated action potential has a duration of
about 300 ?s after the time is scaled, and the capacitances are 1000 times smaller
of the values listed in Table 3.6. For further shrinking the time scale, analysis of
the high frequency operation of the VLSI circuit needs to be done, and this is one
of the future works.
Impact of Parameter Variations on Design
In this section, we analyze the threshold voltage variation sensitivity of the follow-
ing circuits: linear voltage to current converters, linear current to voltage convert-
128
ers, voltage bufiers, and the pool circuits. These circuits are the building blocks of
sigmoid functions and exponential functions, which, due to their exponential prop-
erties, are considered to be very sensitive to the  uctuation of circuit parameters.
Due to the imperfect control of fabrication processes, chips often display varia-
tions in transistor characteristics, such as the threshold voltage and drain currents
[78][79]. Variabilities are roughly classifled into two types [80]. One is chip-to-chip
variability, which results from factors such as processing temperature and equip-
ment properties. The other is within-chip  uctuations, which is caused by random
factors, such as the spatial distribution of dopant and channel length variation. In
the following discussion, we focus on the impact of within-chip threshold voltage
variation on our VLSI design for the heart model. We consider that the thresh-
old voltage is one of the most critical parameters that in uences the accuracies
of circuit realization of mathematical equations, because many output vs. input
equations introduced in section 3.3 are dependent on the threshold voltage, and,
current mirrors, one of the most often used circuit structures, require the proper
match of the threshold voltages of the transistors at the input and output stages.
The Monte Carlo statistical tool provided by PSpice is used to investigate the
threshold voltage variation sensitivity of the heart model circuits. The Monte Carlo
analysis varies device parameters randomly and independently within a specifled
tolerance and performs multiple runs of analysis (DC, AC, or transient). In our
simulation, the tolerance of the threshold voltage is determined from 62 sets of
test data (MOSIS run numbers T18K to T1CL) [81] from difierent fabrication
runs for the AMI Semiconductor 1.5 ?m ABN technology. Figure 3.42 shows the
histogram of the NMOS threshold voltage and the Gaussian fltting curve for this
experimentally determined data. The mean of the Gaussian curve is 0.53788 V,
129
Threshold Voltage [V]
Count
s
Figure 3.42: Histogram and Gaussian flt for the NMOS threshold voltages mea-
sured from 62 sets of data (AMI 1.5 ?m technology).
and the standard deviation  is 0.017 V. Therefore, the variation for our purpose
is taken to be 0:017=0:53788 ? 3:2%. Note that this is a chip-to-chip variation, so
we need to modify it to get the within-chip variation. The within-chip variation is
estimated from the chip-to-chip variation in a way presented in the following.
It was reported in [80] that in a 0.35 ?m process, the variation of the saturation
current within a wafer is 3%, and its variation within 58 lots is 15%. Since a wafer
holds many chips, we assume that the within-chip variation does not exceed the
variation over a single wafer, and the lot-to-lot variation is equivalent to the chip-
to-chip variation. Then the ratio of the chip-to-chip variation to the within-chip
variation for the saturation current is higher by a factor of %15:%3 = 5. Since
the saturation current is dependent on the threshold voltage through a square law,
thus, the saturation current variation is worse than the threshold voltage. Also,
by considering that the within-chip component of the variability becomes larger
when the process scales down [80], we estimate that the within-chip variation of
the threshold voltage for AMI 1.5 ?m technology is lower than one flfth of the
130
Table 3.7: Monte Carlo analysis of threshold voltage variation on circuits
Circuit Function Figure Number  Variation
Voltage-Current Converter Figure 3.5(a) 0.12 ?A 2.1 %
Current-Voltage Converter #1 Figure 3.6 0.265 V 140%
Current-Voltage Converter #2 Figure 3.44(a) 1.3 mV 0.1 %
Voltage Bufier Figure 3.7 1.5 mV 0.2%
Pool Circuit Figure 3.9 3.9 mV 1.0%
Iin [?A]
Vou
t [
V]
Figure 3.43: Monte Carlo analysis of Current-Voltage Converter #1 (Figure 3.6).
131
chip-to-chip variation, that is 3:2%=5 ? 0:6%. Therefore, 0:6% is taken as the
tolerance of the threshold voltage in our simulation, and we use the Gaussian as
the distribution type.
The results of the Monte Carlo analysis of the threshold voltage variation on
the heart model circuits are provided in Table 3.7. The flrst column lists the
circuit functions. All the circuits have been introduced in section 3.3, except
the Current-Voltage Converter #2, which will be introduced shortly. The second
column shows the flgure numbers of the circuits. The third and the fourth columns
are average standard deviations and average variations in percentage of circuit
outputs caused by the  uctuation of the threshold voltage. The results given in
Table 3.7 correspond to the Monte Carlo analysis of 20 PSpice runs.
As we can see from Table 3.7, most circuits we analyzed presented small varia-
tions, except Current-Voltage Converter #1 introduced in section 3.3.1. As shown
in Figure 3.6, the linear Vout vs. Iin relationship relies on the channel modulation
efiect. That is, when there is a small input current, the current difierence caused
in M2 and M4 are canceled by the shifting of the output voltage. Due to the
high sensitivity of Vout to current changes at the output stage, the implemented
transimpedance is high, and the proper functioning of the linear current to voltage
converter demands that M1 and M2, and M3 and M4 match very well. The mis-
match introduced by the threshold voltage  uctuation severely shifts the output
voltage as shown in Figure 3.43, hence, this presented circuit can not tolerate the
variance of the threshold voltage we specify.
Figure 3.44(a) shows a proposed alternative linear current to voltage converter.
It is the same as the input stage circuit of a bi-directional current mirror. By using
equations of the saturation currents on both transistors as the following (under the
132
Vdd
VoutIin
I1
I2
Vss
M1
M2
(a) (b)
Iin[?A]
Vou
t [V
]
Figure 3.44: Proposed linear current to voltage converter. (a) Circuit of linear
current to voltage converter. (b) Monte Carlo analysis: 20 runs.
no load DC conditions, for which it is used):
8>
>>>
>><
>>>>
>>:
I1 = Kp ?(Vdd ?Vout ?jVtpj)2;
I2 = Kn ?(Vout ?Vss ?Vtn)2;
I1 +Iin = I2;
(3.141)
we derive a second order polynomial equation that contains the output voltage and
the input current:
A?V 2out +B ?Vout +C ?Iin = 0; (3.142)
where 8
>>>
>>><
>>>>
>>:
A = Kn ?Kp;
B = ?2?(Kp(Vdd ?jVtp)j+Kn(Vss +Vtn));
C = Kp(Vdd ?jVtpj)2 +Kn(Vss +Vtn)2:
(3.143)
When Kn is close to Kp, Vout can be approximated by the linear term of a Taylor
series expansion to give a linear function of Iin. Then the change of Vout to Iin
ratio is expressed by:
?Vout
Iin = ?
1
2?(Kp(Vdd ?jVtpj)+Kn(Vss +Vtn)); (3.144)
133
which is constant in terms of Iin. The ofiset can be obtained from equation (3.142)
by setting Iin = 0:
VoutjIin=0 = Kp(Vdd ?jVtpj)
2 +Kn(Vss +Vtn)2
2?(Kp(Vdd ?jVtpj)+Kn(Vss +Vtn)): (3.145)
Figure 3.44(b) shows the simulation results of Vout vs. Iin for the proposed linear
current to voltage converter with the transistor parameters being set to WM1:LM1=
16?:8? and WM2:LM2=6.4?:32?. Its numerical results of the Monte Carlo analysis
are provided in Table 3.7 in the row for Current-Voltage Converter #2.
By using the Monte Carlo analysis results shown in Table 3.7, we are able to
estimate the deviation of the outputs of sigmoid function circuits and NPN expo-
nential function circuits. The circuit structure of sigmoid functions is illustrated
in Figure 3.12. Similar to how we obtained the equations (3.55) and (3.56) for Vb
and Vc (i.e. Vd and Ve when the threshold voltage variation is not considered), we
re-derive the equations of Vd and Ve by taking into account the variation of the
threshold voltage as follows:
Vd = Vb +?Vbuf1 = (Iah1 +bH1 +?VH1)+?Vbuf1 (3.146)
= (x?g1 +bG1 +?IG1)h1 +bH1 +?VH1 +?Vbuf1
= Vd0 +h1?IG1 +?VH1 +?Vbuf1;
where g1, h1, bG1, and bH1, have the same meanings as for equation (3.55), and,
thus, are the transconductance and transimpedance of G1 and H1, and the output
ofisets of G1 and H1, Vd0 is the original Vd (equal to Vb) expressed by equation
(3.55), and ?IG1, ?VH1, and ?Vbuf1 are the random difierences of the outputs of
G1, H1, and V Bufier1 caused by the threshold voltage  uctuation. The equation
for Ve is re-derived in the same way:
Ve = Ve0 +h1?IG2 +?VH2 +?Vbuf2: (3.147)
134
?Iout
?Vrandom
0 VT 2VT 3VT-VT-2VT-3VT
0.8xI0
0.6xI0
0.4xI0
0.2xI0
0
1.0xI0
Iout
?V
Figure 3.45: Difierence of Iout in sigmoid function circuit caused by variation of
the threshold voltage. ?Vrandom=VT = 0:3.
where Ve0 is the original Ve (equal to Vc) expressed by equation (3.56), and ?IG2,
?VH2, and ?Vbuf2 are the random difierences of the outputs of G2, H2, and V
Bufier2.
Due totherandom variationattheinputs ofthe Emitter Coupled Pair inFigure
3.12, Vd and Ve, the output vs. the input curve has a horizontal shift ?Vrandom,
which causes a non-deterministic difierence of the output value ?Iout, as illustrated
in Figure 3.45. From equations (3.146) and (3.147), ?Vrandom is expressed by:
?Vrandom = (h1?IG1 +?VH1 +?Vbuf1)?(h1?IG2 +?VH2 +?Vbuf2): (3.148)
Bytaking?Vrandom intotheequationoftheemittercoupledpair, whichisrewritten
here from equations (3.53) and (3.54)
Iout = I0
e(??VVT ) +1
; (3.149)
where ?V is Vd?Ve in this case, we can estimate the upper bound of the percentage
difierence j?Iout=Ioutj in terms of ?Vrandom=VT. The maximum j?Iout=Ioutj occurs,
as seen from Figure 3.45, when the deterministic parts of the emitter coupled pair
135
inputs are equal, i.e. Vd0 = Ve0. Therefore, we have:
max(j?IoutI
out
j) = j( I0e0 +1 ? I0
e(0+
??Vrandom
VT ) +1
)? I0e0 +1j
= j1? 2
e(
??Vrandom
VT ) +1
j?j?Vrandom2V
T
j; (3.150)
wherewereplacee?(?Vrandom=VT) byitslinearTaylorapproximation1??Vrandom=VT,
assuming j?Vrandom=VTj < 1. By using the data of the standard deviation in Table
3.7 (here we use Current-Voltage converter#2) and the slope of the line in Figure
3.44(b) for h1, we can estimate the magnitude of ?Vrandom from equation (3.148):
j?Vrandomj? 2?(8500?0:12?+1:3m+1:5m) ? 7:64mV. Hence, the upper bound
of the percentage difierence max(j?IoutIout j) is about 7:64mV=2VT = 14:7%. Figure
3.45 shows a normalized sigmoid curve in the case when ?Vrandom=VT = ?0:3. The
corresponding output difierence is 15%. We still use R-square and signal-to-noise
ratio (SER) to evaluate how close the shifted curve is to the original curve. The
R-square and SER for the curves in Figure 3.45 are 0.9903 and 22.8 dB respectively.
Using the same analysis methods as presented above, we can derive the equation
for the percentage difierence of the output of our NPN exponential circuit shown in
Figure 3.13. The equation for the base-emitter voltage of Q2 is, by taking account
of the threshold voltage variation, expressed as follows:
Vbe;Q2 = Vb ?Vc +Vd +?Vbuf +?Vpool
= (Vb0 +h1?IG1 +?VH1)?(Vc0 +h1?IG2 +?VH2)+Vd +?Vbuf +?Vpool
= Vbe;Q2;0 +(h1?IG1 +?VH1)?(h1?IG2 +?VH2)+?Vbuf +?Vpool
= Vbe;Q2;0 +?Vrandom; (3.151)
where ?IG1, ?IG2, ?VH1, ?VH2, ?Vbuf, and ?Vpool are the random variation
of the outputs of G1, G2, H1, H2, V Bufier, and Pool Circuit caused by the
136
 uctuation of the threshold voltage, and Vbe;Q2;0 is the base-emitter voltage of
Q2 when the threshold voltage  uctuation is 0. By using the data in Table 3.7
and Figure 3.44(b), we estimate the magnitude of ?Vrandom to be: j?Vrandomj ?
2?(8500?0:12?+1:3m)+3:9m+1:5m = 10:0mV.
The output percentage difierence caused by ?Vrandom is derived as follows:
j?IoutI
out
j = j(Is0e
Vbe;Q2;0
VT ?Is0e
Vbe;Q2
VT )?Is0e
Vbe;Q2;0
VT j
= j1? Is0e
Vbe;Q2
VT
Is0e
Vbe;Q2;0
VT
j
= j1?e?Vrandomj?j?VrandomV
T
j; (3.152)
wherewealsouseTaylorexpansiontoapproximatetheexponentialterme?Vrandom ?
1+?Vrandom. By using the ?Vrandom value estimated previously, the percentage dif-
ference can be obtained from equation (3.152), which is about ?VrandomV
T
= 38:5%.
The corresponding R-square of the original exponential curve and the random
shifted curve is 0.8689, and SER is 9.1 dB.
As analyzed above, the output percentage difierences caused by the threshold
voltage randomness of the MOS transistors in sigmoid function circuits and NPN
exponential circuits are highly dependent on the thermal voltage VT. By increasing
VT, the percentage difierence can be reduced. Despite of the fact that VT can
not be adjusted manually, we can increase the numerators of equations (3.150)
and (3.152) by factors of integers by changing the circuit structures of the emitter
coupled pair, and the NPN transistor circuitries in the exponential function circuit.
Figure 3.46(a) shows the new emitter coupled pair, modifled for enhancing circuit
immunity to device parameter changes. This circuit, compared with the one in
Figure 3.11, uses two levels of cascode NPN devices. The equations of the output
137
M1 Vbias
Vss
Vdd
Q1a Q2a
M2 M3
V+ V-
I0
Va
M4 M5
I+ I-
Q1b Q2b
Iin
Iout
Q1a Q
2a
Q1b Q
2b
From V Buffer
(a) (b)
Vd
Figure 3.46: (a)Modifled emitter coupled pair using cascode NPN devices;
(b)Cascode NPN structure used in exponential function circuit.
currents are changed from equations (3.53) and (3.54) to:
Ic1 = I0
e(? ?V2VT ) +1
; (3.153)
Ic2 = I0
e( ?V2VT ) +1
; (3.154)
which are equivalent to the expressions for the case when VT is magnifled by
two times for the original emitter coupled circuit. The modiflcation for the NPN
exponential function circuit is shown in Figure 3.46(b). Cascode-connected Q1a
and Q1b replace the Q1 in Figure 3.13, and Q2a and Q2b replace the Q2. The
expression of Iout in Figure 3.13 then becomes:
Iout = Iine
g1h1(Vin?Vbias)
2VT :
With the same methods, we can connect more than two NPN transistors in cascode,
with the limitation being that of power supply voltages. As a result, the output
percentage difierence caused by circuit parameter variation, expressed by equations
(3.150) and (3.152), can be reduced by the number of the cascode-connected NPN
transistors.
138
Chapter 4
Propagation of Electrical Activity
in Cardiac Tissue
4.1 Introduction
The specialized excitatory and conductive cardiac muscles are capable of conduct-
ing electrical signals, i.e. action potentials, throughout the heart, and incorporate
an excitatory and transmission system which is tightly coupled with the mechanical
activity and responsible for the contraction and relaxation of the heart muscle.
Cardiac cells are typically cylindrical, 30-100 ?m long and with radius 8-20 ?m
[36][38]. They connect end to end to form elongated flbers, which are arranged in
parallel to construct flbrous bundles. The two ends of a cardiac cell are separated
from their neighbors by intercalated disks, the membranes that separate adjacent
cells in cardiac tissues. The conduction of electrical signals from one cell to the
next cell depends on gap junctions (nexi), very permeable pathways present in the
intercalated disks, allowing relatively free difiusion of ions and thus the passage of
the electrical signals. Gap junctions exist in the borders between cells connected
139
longitudinally, and they are sparse or absent in the borders between cardiac flbers
that lie side by side. Therefore, current  ows more easily in the direction along
the longer axis than the transverse directions [82].
Two difierent approaches of modeling the spread of the electrical activity have
been developed by treating cardiac tissues as either a discrete system of coupled
cells, or a continuous medium, called cellular automata and core-conductor models
respectively [5][83]. The cellular automata is a highly abstract model, in which a
discrete network representing the spatial structure of cardiac tissue is constructed,
with each node (a cardiac cell) taking one of a flnite set of states at each moment.
The state of each node evolves according to deterministic rules and is decided
by the last states of this node and its neighbors. Because of the simplicities of
the cellular automata model, it is capable of e?ciently simulating the excitation
propagation throughout the whole heart. However, the model is not related to
cellular electrophysiology and, thus, can not describe biophysical details of cel-
lular excitability [83]. The weakness of the automata model is overcome by the
core-conductor model, also known as excitable dynamics equations, which, by com-
bining the models of cardiac ionic membranes, describes the cardiac propagation
behaviors using coupled PDE-ODE (partial difierential equation - ordinary difier-
ential equation) equations. The core-conductor models incorporate the greatest
amount of present knowledge about the ionic basis for the action potential, and
thus become dominant computational approaches to study wave phenomena in
cardiac tissue for today.
In this chapter, we will present how to obtain the numerical results of the partial
difierential equation of the cardiac propagation model using VLSI circuits. The
propagation model belongs to a special category of partial difierential equation,
140
namely reaction-difiusion equations. The circuit mapping of the generalized form
of reaction-difiusion equation is introduced, and we then show the circuit setup for
simulating the electrical propagation over a small piece of ventricular muscle using
simplifled spatial conflguration of the cardiac muscle. The modeling of the cardiac
cell structure details, such as the orientation of cardiac flbers, cellular geometry,
and the anatomic structure of the heart, are beyond our efiorts. Hence, here we
do not intend to construct a circuit to simulate the electrical propagation over the
whole heart, but propose a methodology of calculating reaction-difiusion equations
with analog circuits, especially for the cardiac system.
In the rest of the chapter, we will flrst introduce the model of 1-dimensional
cable theory and its extended multi-dimensional equations which describe the prop-
agation of cardiac electrical activity using partial difierential equations. Next we
will present the circuit mapping of reaction-difiusion systems, a methodology of
computing reaction-difiusion equations with circuits. It is worth to mention that
the proposed methodology is also applicable to the reaction-difiusion systems that
are irrelevant to electrical phenomenons, such as population model of genetics
[84][85], model of morphogenesis [86], and model of the Belousov-Zhabotinsky re-
action [87][88]. The simulation of a VLSI circuit and an ideal component circuit
that implement cardiac propagation is provided flnally.
4.2 Electrical Propagation Model of the Heart
The core-conductor class of models are evolved from the cable theory [89][90],
which was developed more than a century ago by Kelvin in [91] to describe the
transmission of an electrical telegraph. The cable equations were flrst modifled
to apply in the biological sciences in 1879 by Hermann [89][92] for analysis of
141
C
Vm a
Ip
Iin Iout
Iq
m
?x
(a) (b)
A cardiac muscle fiber
A cardiac muscle fiber in cable model
Extracellular
Intracellular
Vin Vout
Figure4.1: 1-Dimensional cablemodelforcardiactissue. (a)A6-cellcardiacmuscle
flber (top) and its representative (bottom) in cable model - a continuous media
that does not consider the boundaries of cells. (b)Circuit diagram of cable model.
nerves. Since then such models have been proven successful in studies of electrical
propagation in nerves [7]-[9] and cardiac muscle [12][82][83] up to today.
4.2.1 1-Dimension Cable Model in Cardiac tissue
The 1-dimensional cable model prototypes a single flber and characterizes the
electrical propagation over a bundle of parallel muscle flbers [31]. The cable model
is a continuous model that does not consider the boundaries of cells, as illustrated
in Figure 4.1(a). In the top of Figure 4.1(a) shows a portion of a cardiac flber that
contain 6 cells, and in the bottom of the flgure is the representative of the flber in
the cable model, a continuous and uniform cylindric medium.
Consider a segment of cardiac flber in the cable model (let us call it a cardiac
cable to distinguish it from a real cardiac flber which has cell-to-cell connections),
as shown in Figure 4.1(b). There are three elements of currents in the cardiac
cable: 1) the ionic current, denoted with Iq in the flgure, that  ows between the
intrcellular space (the cable) and the extracellular space (the ground); 2) Ip, the
142
current of the membrane capacitor C; 3) Iin and Iout, the current  owing through
the resistance caused by gap junctions and cytoplasm 1. Ip and Iq can be expressed
by the membrane current per unit (cable surface) area Im and the ionic current
per unit area Iion multiplied with the cable surface area of the membrane segment:
Ip = C@Vm@t = 2?a?x?Cm@Vm@t ; (4.1)
Iq = 2?a?x?Iion; (4.2)
where a is the radius of the cable, ?x is the length of the small segment of cardiac
cable, and 2?a?x is the lateral area of the cable. Iin and Iout can be formulated
using Ohm?s law:
Iin = ?Vin=R = ?Vin=(??x?a2 ); (4.3)
Iout = ?Vout=R = ?Vout=(??x?a2 ); (4.4)
where R is the resistance of the cable segment which is determined by the resistivity
? times the length ?x and divided by the area of the cross section ?a2. Putting
equations (4.1)? (4.4) together by Kirchhofi?s current law, which is formulated as:
XI = I
in ?Iout ?Ip ?Iq = 0 ) Iin ?Iout = Ip +Iq; (4.5)
we derive:
?a2
?
?Vin
?x ?
?a2
?
?Vout
?x = 2?a?xCm
@Vm
@t +2?a?xIion: (4.6)
By letting ?x ! 0, equation (4.6) becomes a continuous difierential equation:
a
2?
@2Vm
@x2 = Cm
@Vm
@t +Iion: (4.7)
1Cytoplasm is the  uid enclosed by a cell membrane and fllled with dispersed particles and
organelles [2](pp. 10-12).
143
Replacing resistivity ? by the conductivity d = 1=?, and representing a=2 by the
inverse of the surface-to-volume ratio fl = (2?a?x=?a2?x) = 2=a, we obtain the
cable model of 1-dimensional cardiac flber, written as:
d
fl
@2Vm
@x2 = Cm
@Vm
@t +Iion; (4.8)
which is:
1
flr?(drVm) = Cm
@Vm
@t +Iion: (4.9)
Assumption of Cable Model
The basis of the cable model is that it treats a cardiac flber as 1-dimensional and
reduces the problem to a single spatial dimension [93]. 1-dimension is only in the
sense that the current and voltage potentials are uniformly distributed across any
cross section perpendicular to the longitudinal direction of the cable [94].
The resistance shown in Figure 4.1(b) corresponds to the average resistance
of the intracellular space, which in reality is composed of cytoplasmic and gap
junctional resistances [95][96]. This model, called the continuous model, therefore
considers the cardiac tissue as a medium with continuous difiusive properties and
is only capable of preserving macroscopic details [90]. The gap junctions, which
are relatively discrete due to their short length but sizable resistance compared to
cytoplasmic [31], introduce jumps in the voltage patterns among cardiac cells, and
this can be modeled by a variation of the original cable model, in which, the gap
junctions are represented by a change in resistance [97]-[99].
The presented model assumes that the resistivity of the extracellular media is
negligible and can be treated as isopotential, and only considers the propagation
of electrical signals within the cells by taking the extracellular space as grounded
as shown in Figure 4.1. This assumption relies on the experimental observation
144
that the extracellular medium has a much higher conductivity compared to the
intracellular medium (a typical factor of 8 is reported in the longitudinal direction
in [100]).
4.2.2 Multi-Dimensional Cardiac Propagation Model
The multi-dimensional cardiac propagation model is extended from equation (4.9)
and takes into account the anisotropic electrical properties of cardiac media. This
type of equation is also called an excitable dynamics equation and is formulated
with a reaction-difiusion equation as given below [5]:
r?(DirVm) = fl(Cm@Vm@t +Iion); (4.10)
where Vm, fl, Cm and Iion, the same as in equation (4.8), are the membrane poten-
tial, the surface-to-volume ratio, the membrane capacitance, and the transmem-
brane current per unit area. Di is the intracellular conductivity tensor, which
incorporates the conduction properties of the intracellular space. Di is a single
variable d in equation (4.8) for the 1-dimensional case, and for 2-dimensional and
3-dimensional cases, it is expressed by:
2
66
4
d  d ?
d? d??
3
77
5; (4.11)
and
2
66
66
66
4
d  d ? d ?
d? d?? d??
d? d?? d??
3
77
77
77
5
; (4.12)
respectively, where  , ?, and ? are three axes in the selected Cartesian coordinate
system. When the cardiac flbers are considered as parallel and straight, it is
145
mathematically convenient to adopt the Cartesian coordinate system and align it
with the principal axes of Di to make Di a diagonal matrix. When modeled cardiac
tissues are curved or rotated, Di generally cannot be diagonalized by selecting the
axes [36].
The most commonly applied boundary condition for equation (4.10) is a no-
 ux, or Neumann condition, expressed as:
(DirVm)?~n = 0; (4.13)
where ~n is the normal vector to the boundary. The no- ux condition, also known
as sealed boundary condition, deflnes no current  ows into or out of the edge of a
piece of cardiac muscle (including intracellular and extracellular space), and can
represent naturally the outer surface of the heart.
4.2.3 Summary of Cardiac Propagation Model
The present quantitative description of electrical signal propagation in cardiac tis-
sue is based on two classical models. They are the Hodgkin-Huxley model, which
describes the active membrane properties using the kinetics of the ionic currents,
and the cable equation of Kelvin, which relates the current and voltage along
a continuous one-dimensional structure [90]. For the propagation model we se-
lect for VLSI implementation, the Beeler-Reuter model, which has been designed
with VLSI circuits as introduced in chapter 3, is chosen to describe the mem-
brane excitable properties and the non-linear ionic currents, and the continuous
core-conductor model expressed by equation (4.10) is employed to formulate the
electrical activities transmitted in the intracellular space.
Put all the equations together for our propagation model, equation (4.10),
replacing equation (2.40), becomes the top level of formulation, with its boundary
146
Simplified cardiac 
muscle
Cardiac muscle as continuous medium in 
cable-theory-based model
z
y
x
Figure 4.2: Simplifled structure of cardiac tissue and its continuous representative
in cable-model.
condition being represented by equation (4.13). All the other equations in the
Beeler-Reuter model, from (2.41) to (2.60) are still valid. Note that in equation
(4.10), the propagation term r?(DirVm) substitutes Iext in the equation (2.40).
As explained in chapter 2, Iext represents the external in uence applied on a cell
and can be used as a stimulus to trigger the action potential. Therefore, in the
propagation model, it is the current transmitted along the myocardium flber that
activate a cardiac cell to create an action potential.
For the circuit implementation, we simplify the structure of cardiac flbers to
be as in the left side of Figure 4.2, where we assume all the muscle flbers are
parallel and straight by ignoring the curved shapes of the heart muscles and the
difierent orientations between difierent cardiac flbers. Therefore, the intracellular
conductivity tensor Di can be made a diagonal matrix by properly selecting the
coordinates. On the right side of Figure 4.2 is the representative of the cardiac
muscle in the cable-theory-based model, which is a continuous medium that does
not difierentiate the boundaries of individual cells. As illustrated in Figure 4.2,
we select the Cartesian coordinate system and assume that x and y, orthogonal
to each other, are both transversal to the flber axis; and z is parallel to the flber
longitudinal direction. Expanding the nabla operator on the diagonalized Di,
147
equation (4.10) becomes equation (4.14):
dx@
2Vm
@x2 +dy
@2Vm
@y2 +dz
@2Vm
@z2 = fl(Cm
@Vm
@t +Iion); (4.14)
where dx, dy and dz are the conductivity in the direction x, y and z respectively,
all assumed to be constant.
In section 4.3, we will show the circuit simulation of the electrical propaga-
tion over the cardiac tissues based on equation (4.14). We will only present the
simulation for 1-dimensional and 2-dimensional propagation, due to the lack of
the transillumination imaging techniques [101] needed for the visualization of 3-
dimensional propagation. The propagation over 3-dimensional cardiac tissues can
be realized by circuits using the same methodology, which will be described in
section 4.3.1.
4.3 Modeling Propagation of Cardiac Active Po-
tential with VLSI
Core-conductor models involve the fewest simplifying assumptions and incorpo-
rate the greatest amount of present knowledge about the ionic basis for the action
potential, and, thus, are widely used for the study of cardiology. Since the 1950s,
core-conductor models have been descretized to divide a continuous flber model
into segments to allow the studies of propagation with circuits and later with
high-speed computers [20]-[22][37][102]. However, the ionic-membrane-based car-
diac cell model described by the reaction-difiusion equation (4.10) (and equation
(4.14)), can substantially increase the complexity to obtain a numerical solution
with computers [96]. Equation (4.10) formulates a special coupled PDE-ODE (par-
tial difierential equation - ordinary difierential equation) system, with the PDE
148
formulating the propagation and the ODE presenting the membrane ionic mecha-
nism. In the heart model, the ODE system and the PDE system, though coupled
through the voltage and current, are independent of each other. At any point in
the space, the activation of Vm, or the generation of the action potential, relies
on the spatial conduction of electrical signals described by the PDE to transmit
the stimulating current, whereas, the shape of the action potential at that point
is totally decided by the ODE equations. The ODE system determines the sole
dependance of the gating variables and ionic concentrations on Vm at each spa-
tial point, and, thus, their boundary conditions are not available. Because of the
complexity of the problem, computation with digital computers is very expensive
and challenging. The rapid upstroke in an action potential of a ventricular cell
can require the time discretization on the order of 0.01 ms [103], and the steep
wave front in space may require the space discretization on the order of 0.1 mm.
This means to simulate 1 s of the heart beat on a uniformly grided 3-dimentional
tissue of size 1 cm3, there are 1011 unknown Vm. If the Beeler-Reuter model (a
degree-eight ODE system) is used to formulate the membrane ionic mechanism,
there are 8?1011 unknown values coupled together through time and space. The
tremendous complexity in computation with digital computers promotes the use
of analog circuits to perform the calculation.
Simulating propagation with circuits can be traced back to Schmitt?s circuit of
the electronic neural network implementation [104], and the active pulse transmis-
sion line constructed by Nagumo etc. [9][105]. Nagumo etc. used a capacitor, a
tunnel diode, and an inductor to simulate the ionic behavior of a segment of an
animal nerve axon; the segments are coupled with resistors to construct an active
transmission line. The purpose of segmenting the continuous propagation model
149
for circuit implementation is just the same as meshing numerical regions for flnd-
ing solutions with computers, which is to use a flnite number of data to present a
continuous space.
Despite the great similarity of the circuit implementation and the presented
cablemodel (botharenonlinear RCnetworks), theproposed methodology ofcircuit
computation is valid for calculating not only the electrical propagation models for
cardiac tissue, but also applicable to other reaction-difiusion systems in biology,
physics, chemistry, etc [84]-[88].
4.3.1 ComputationofDiscretizedReaction-DifiusionEqua-
tion with Circuits in Cartesian Coordinates
The partial difierential equation for a general reaction-difiusion system is formu-
lated as [87][106][107]:
r?(Ajruj) = bj@uj@t +fj(u1;:::;uM); (j = 1;2;:::M); (4.15)
where uj are characteristics of the difiusion medium, bj are constants, fj, rep-
resenting the reaction in the reaction-difiusion system, are functions of uj, and
Aj are difiusion tensors. The presented set of partial difierential equations (4.15)
describe a system that has M difiusion variables, which, with each being time-
and space-dependent, are coupled together to make the system cooperative by the
reaction functions fj. We use Cartesian coordinates and consider only the case
when the difiusion tensors are flxed, i.e., they do not change with the space axis,
and all Aj are isotropic. Hence, the Cartesian coordinates can be selected to make
150
C
I0
Va
a
RyR
x
Rz
Iz-
Iz+
Iy-
Iy+
Ix-
Ix+
IC
Figure 4.3: Using an RC network to implement reaction-difiusion systems.
Aj diagonal matrices which are expressed as:
2
66
66
66
4
ajx 0 0
0 ajy 0
0 0 ajz
3
77
77
77
5
; (4.16)
where ajx, ajy, and ajz, are the difiusion coe?cients for the directions x, y and z. In
the following discussion, we will show how the RC network shown in Figure 4.3 can
be used to implement the reaction-difiusion system formulated by equations (4.15)
(and (4.16)) after (4.15) is discretized. Note that the circuit shown in Figure 4.3
is to realize one of the equation set formulated by (4.15). To represent the whole
system, M of the same RC circuit with difierent parameters are needed, and all
circuits are coupled together through the current sources I0 in Figure 4.3. The
following discussion is for the jth equation in the reaction-difiusion system.
Equation (4.15) can be rewritten by expanding the nabla operation as follows:
ajx@
2uj
@x2 +ajy
@2uj
@y2 +ajz
@2uj
@z2 = bj
@uj
@t +fj: (4.17)
The partial difierential equation (4.17) can be discretized spatially by letting @i (i
151
is x, y, or z) transition to segments ?i with non-zero lengths, and the equation
then is changed to:
(ajx?ujx+?x2 ?ajx?ujx??x2 ) + (ajy?ujy+?y2 ?ajy?ujy??y2 )
+(ajz?ujz+?z2 ?ajz?ujz??z2 ) = bj@uj@t +fj; (4.18)
where ?uji+ and ?uji? (i is x, y, or z) are the changes of uj in the positive and
negative i direction. Equation (4.18) becomes equation (4.17) when the ?i (i is x,
y, or z) tend to be 0.
Let us now consider the 3-dimensional RC network shown in Figure 4.3, in
which each node is joined to ground through a capacitor C and a current source
I0, and distributed nodes are bridged with resistors. At node a, from Kirchhofi?s
current law, the sum of the currents that  ow into the point is 0, and, hence, the
following equation is fulfllled:
(Ix+ ?Ix?)+(Iy+ ?Iy?)+(Iz+ ?Iz?) = Ic +I0; (4.19)
where Ic is the capacitor current, and Ii (i is x+, x?, y+, y?, z+, or z?.) are
the currents in the six directions as illustrated in Figure 4.3. By using Ohm?s law
and the I-V property of the capacitor, equation (4.19) becomes:
(?Vx+R
x
? ?Vx?R
x
) + (?Vy+R
y
? ?Vy?R
y
)
+(?Vz+R
z
? ?Vz?R
z
) = C@Va@t +I0; (4.20)
where ?Vi+ and ?Vi? (i is x, y, or z) are the voltage difierences in the positive
and negative direction i, and Ri (i is x, y, or z) are the resistance in direction i.
There is a great similarity between equations (4.18) and (4.20), and thus, (4.20)
can be made the same as (4.18) by letting the circuit parameters satisfy:
Rx = ?x
2
ajx ; (4.21)
152
Ry = ?y
2
ajy ; (4.22)
Rz = ?z
2
ajz ; (4.23)
C = bj; (4.24)
I0 = fj: (4.25)
Note that I0 is normally variable due to the dependence of fj on u1;u2;:::;uM.
Since equation (4.18) comes from (4.15), the RC circuit shown in Figure 4.3 with
parameters selected following equation (4.21)-(4.25) describes the same reaction-
difiusion system as equation (4.15) when ?i ! 0 (i = x;y;or z).
Equations (4.21)-(4.25) is not a unique solution for mapping a reaction-difiusion
system into the RC circuit in Figure 4.3. Multiplying by a constant on both sides
of equation (4.15) creates another equation that describes exactly the same system.
Hence the circuit mapping equations provided by (4.21)-(4.25) can be expanded
to be the following:
8>
><
>>:
RiC = (?i)2 ? bjaji;
RiI0 = (?i)2 ? fjaji;
(i = x; y;or z): (4.26)
The 3-dimensional RC network shown in Figure 4.3 can be changed to implement
reaction-difiusion systems of 1-dimension or 2-dimension by removing the parts
of the circuit in the direction(s) that does not exist in the corresponding system
as shown in Figure 4.4. Figure 4.4(a) illustrates a 1-dimensional RC circuit with
3 discretized non-ground nodes. Figure 4.4(b) depicts the RC circuit for a 2-
dimensional case, and there are 4 nodes in the flgure. Figure 4.4(a) and (b) describe
the reaction-difiusion systems formulated with:
ajx@
2uj
@x2 = bj
@uj
@t +fj; (4.27)
153
Rx
Ry
C
I0
IC
C
I0IC
Rx a
a
(a) (b)
Figure 4.4: Using an RC network to implement reaction-difiusion systems in the
case of (a) 1-dimensional, (b) 2-dimensional.
and
ajx@
2uj
@x2 +ajy
@2uj
@y2 = bj
@uj
@t +fj; (4.28)
with the circuit parameters being still restricted by equation (4.26).
The feasibility of using a circuit to represent a reaction-difiusion system re-
lies on the fact that an RC network is also a reaction-difiusion system, and all
reaction-difiusion systems actually describe scenarios of the same nature, which is
the movement of a certain physical substance over a particular media through a
difiusion process. The boundary conditions are similar for most reaction-difiusion
systems. Two common boundary conditions are a) the no- ux condition, a situa-
tion where substances can not move across the boundary, and b) the flx amount
condition, which re ects the situation when the substance amount is flxed by ex-
ternal conditions [107].
The proposed circuit mapping methodology is limited by the requirement that
the space is uniformly gridded for individual directions. The limitation arises from
154
the derivation of equation (4.20), where the resistances Ri (i is x, y, or z) are
constants and do not change by varying positions. Due to the relation between Ri
and ?i (i is x, y, or z) deflned by equation (4.26), ?i needs to be flxed too. The
grid sizes of difierent directions, however, do not need to be equal, i.e. Rx, Ry, and
Rz can be difierent.
4.3.2 Transistor Circuit Simulation of 1-Dimensional Car-
diac Propagation
The normal pattern of propagation in the heart, from the sinoatrial node, through
the atrium, to the atrioventricular node, to the bundle of His to the ventricular
muscle, is an essentially 1-dimensional sequence [83]. Thus the simulation of a
1-dimensional propagation system, i.e. a single cardiac flber, gains the knowledge
of normal cardiac activity and is helpful to understand the electrical system of the
whole heart. A 1-dimensional cardiac model is simpler than a multi-dimensional
one and is usually constructed by fewer discretized nodes and, thus, consumes
less computational resources and time to simulate. Hence it is widely used by
researchers to obtain prototyping results. In this section we show the simulation
of the 1-dimensional cardiac propagation using an RC circuit.
Figure 4.5 illustrates the circuit setup for simulating the electrical signal propa-
gation with a 1-dimensional model following the method presented in the previous
section, and, in this case, there is only one difiusion variable, i.e. M = 1 in
equation (4.15). Each pair of the membrane capacitor Cmi and the ionic cur-
rent Ii (i = 1;2:::N) composes a single ionic membrane model that represents a
segment of the cardiac flber. N identical segments are connected with resisters
(R1 = R2 = ::: = RN?1) in series to construct a 1-dimensional model of the heart
155
Cm1Istimulus
Vm1? Ca2+
Na+
K+ I1 Cm2
Vm2?
?
?
Cm, N-1
Vm ,N-1?
IN-1 Cm,N
Vm ,N?
IN
R1 R2 RN-1
Ca2+
Na+
K+
Ca2+
Na+
K+
Ca2+
Na+
K+I2
a segment of cardiac fiber
Figure 4.5: Circuit of 1-dimensional cardiac propagation.
muscle. On the left of Figure 4.5, the current source Istimulus is an external stimulus
and is responsible to create a current to trigger the activation of the flrst cardiac
segment consisting of Cm1 and I1. The excitation is then passed down through the
coupling resistors to the entire chain until it reaches the last node N. The circuit,
without any current  owing into node N from its right side, is naturally consistent
with the no- ux boundary condition of the cardiac propagation model, formulated
by equation (4.13).
The difierential equation of the 1-dimensional propagation system is modifled
from equation (4.14) by making dx = 0 and dy = 0. For the other coe?cients in
equation (4.14), we use the values used by Bray and Wikswo in [108] in our propa-
gation model. The surface-to-volume ratio fl takes 3000 cm?1, dz = 1:863 mS=cm,
and Cm = 1?F=cm2. Iion is deflned by the Beeler-Reuter model provided in equa-
tion (2.41)- (2.60). Equation (4.14) then is rewritten as:
6:21?10?8@
2V 0
m
@z2 = 10
?7@V
0
m
@t +I
0
ion; (4.29)
where we, again, scale down the magnitude of the membrane potential by 100, the
same way as we did in chapter 3, and use scaled V 0m to distinguish it from the
original variable Vm, I0ion is scaled Iion by a factor of 106.
The circuit parameters of the RC implementation for equation (4.29) are se-
156
lected following equation (4.26). ?x is set to 600 ?m, that is, two adjacent nodes
in Figure 4.5 represent a section of 1-dimensional cardiac tissue with 600?m in
length. Following the equations in (4.26), the circuit parameters - the resisters,
the capacitors and the current sources in Figure 4.5, are set as follows:
8
>>>>
>><
>>>
>>>:
Ri = 5:76?104 ?;
Cmi = 0:1 ?F; (i = 1;2:::N);
Ii = I0ion;
(4.30)
Note that the initial conditions of the reaction difiusion equation can be also shifted
to zeros with the method presented in section 3.2.3, for which equation (4.29) stays
unchanged after the initial value shifting process.
The circuit in Figure 4.5 is realized with VLSI by replacing each node with
the implementation circuit of the Beeler-Reuter model presented in chapter 3 (the
top level circuit diagram is shown in Figure 3.19) and substituting the resisters
with their VLSI realization. The implementation of resistors with VLSI devices
has been introduced in section 3.3.1 and shown in Figure 3.5(b). As the action
potential is in the range of [?85 mV;50 mV], the membrane voltage on each node
V 0mi (i = 1;2:::N) is in the range of [0 V;1:35 V] after being scaled and shifted
for the initial conditions. Hence, the transistor-implemented resisters are required
to have linear V-I properties when inputs V+ and V? are independently changing
within [0 V;1:35 V]. The resistor circuit, being non-ideal two-terminal devices, can
experience current difierences in two terminals, that is the current  owing into V+
is difierent from the current  owing out of V? due to the impact from the exter-
nal circuits. This happens when the currents in the output transistors in Figure
3.5(b) do not mirror properly the currents decided by the input stage. Hence,
good current mirrors are important to ensure the resistor properties of the circuit,
157
z
t [ms]
Vm
?[
V]
Figure 4.6: PSpice simulation results of 1-dimensional cardiac propagation con-
structed by transistor circuit.
and, thus, cascode current mirrors are employed for our heart implementation. In
addition, making the output transistors still working in the saturation region when
applied with input voltages (which are also the voltages at the outputs) is another
key point to necessitate I+ = ?I?. The transistor circuit of implementing the
resistors Ri = 5:76?104 ? (i = 1;2:::N) is provided in the schematics page Sch-32
in appendix A.
In our transistor simulation circuit, totally 20 nodes are established, i.e. N = 20
in Figure 4.5. The stimulation Istimulus is set to give a square pulse (impulse) of
2.2 mA for 0.1 ms starting at 10 ms. The corresponding simulation results of
the 1-dimensional model are shown in Figure 4.6. The z axis indicates the node
number, from 1 to 20. The flrst node is excited at t = 10 ms, and its extremely
high upstroke is because of the external stimulus. The excitation is propagated
down to the other nodes, and the last (20th) node is excited at about t = 70 ms.
158
4.3.3 2-Dimensional Cardiac Propagation: Spiral Reentry
In this section, we will provide the circuit simulation for 2-dimensional cardiac
propagation and present the phenomenon of reentry spiral that can be visualized
only in 2-dimensional or 3-dimensional models of the heart. The spiral waves are
life threatening because they can lead to abnormal rapid heart beat, that is not
controlled by the pace maker of the heart, and then to ventricular flbrillation, an
arrhythmia2 leading to sudden cardiac death.
Though our heart muscles are thick enough that the electrical propagation over
the heart is a 3-dimensional process, we treat the problem as 2-dimensional as do
most studies of spatial patterns of electrical activity in the heart [107]. We will
still adopt equation (4.14) in our simulation modifled for the 2-dimensional case.
The VLSI circuit realization is similar to the 1-dimensional case presented in the
previous section, both using the implementation circuit of the Beeler-Reuter model
discussed in chapter 3 and the same transistor circuit for resisters, and, thus, will
not be repeated here. Due to the extremely expensive time cost in simulating the
VLSI circuits with PSpice for the propagation model (for instance, the simulation
of Figure 4.5 with N = 20 took 20 hours for modeling 0.6 s of the heart electrical
activity), the simulation of the transistor circuits to demonstrate a 2-dimensional
spiral wave is unafiordable. Since we have showed the validity of using VLSI
circuits to model the cardiac propagation by applying it in the 1-dimensional case,
in the following simulation, we no long adopt transistor circuits and take use of
ideal component circuits to avoid the overwhelming simulation time demanded by
the transistor heart model.
2Arrhythmia is a group of conditions in which the muscle contraction of the heart is irregular,
faster, or slower than normal [26].
159
Reentry and Refractory Period
Under certain circumstances, a cardiac impulse may re-excite some portion of the
heart through which it had passed previously. This phenomenon is called reentry
[26], which is believed to be the primary mechanism underlying the tachycardia,
an abnormal fast cardiac rhythm, and flbrillation, a state where cardiac muscles
undergo an irregular type of contraction without proper coordination [109]-[115].
The tachycardia and flbrilation can be in the ventricle or the atria, and especially
the ventricular flbrillation is most severe and is the most common cause of sudden
cardiac death which is responsible for about 1 out of 6 deaths in the Western world
[116]. Therefore, the dynamics of the reentry waves, the mechanisms of remedy
methods and drug therapies are of interest in understanding sudden cardiac death
in patients.
Reentry is tightly related to the refractory period. During much of the action
potential, the membrane is unable to flre a second action potential, no matter
how strong an external stimulus is applied to the cell. The state in which a
membrane is unexcitable is called the absolute refractory period. The refractory
state occurs because a large fraction of the sodium channels is inactivated and
cannot be reopened until the membrane is repolarized [26]. This can be seen in
Figure 3.34(c), which shows the sodium current INa is most of the time around
zero because of the closed channels, except at the overshoot of the action potential.
Figure 4.7 illustrates the membrane potential when stimulated within and out of
the refractory period with a second stimulus. The thick line depicts a normal
action potential, activated by a stimulus occurring at t = 100 ms with a strength
0:3 mA=cm2 for 0.1 ms. A second stimulation, the same strength, is applied at
t = 200 ms, t = 300 ms, t = 350 ms, and t = 380 ms separately, and these do
160
Time [ms]
Vm
[mV
]
(1)
(2)
(3)
(4)
(5)
(6)(7)
Figure 4.7: Action potential applied with a second stimulus at difierent moment:
(1) t = 200 ms, (2) t = 300 ms, (3) t = 350 ms, (4) t = 380 ms, (5) t = 400 ms,
(6) t = 410 ms, and (7) t = 420 ms. Stimulus strength: 300 mA=cm2 for 0.1 ms.
not cause a second action potential, as depicted by the dotted curves (1) to (4),
because the stimulations happens during the refractory period. The curves (5) to
(7) correspond to the stimulations, the same strength as the flrst one, applied at
t = 400 ms, t = 410 ms and t = 420 ms respectively, and the waveforms show
that the activated action potentials are more and more like normal ones as the
stimulation is more and more postponed.
There are two types of reentry: anatomical reentry and functional reentry.
Functional reentry involves waves circulating around refractory cardiac tissue, and
this will be discussed further with our 2-dimensional propagation simulation in the
next section, where we show a spiral wave, the form of reentrant pattern being
in 2-dimension. Here we explain reentry with an example of anatomical reentry,
in which a conduction obstacle is involved. The example [26] is shown in Figure
4.8. In each of the two panels, a single bundle of cardiac flbers splits into two
branches, 1 and 2, and they are connected again by bundle 3. Normally, as shown
161
1 2
3
(a) (b)
1 2
3
Figure 4.8: Example of reentry. (a) Action potential travels normally. (b) Reentry
is caused by a unidirectional block.
in Figure 4.8(a), an action potential moving down travels in each branch 1 and 2,
and then enters bundle 3 from both sides. The wave from the left side can not
proceed further because the tissues on the right side are just activated by the wave
from branch 2 and thus are absolutely refractory. The wave from the right side
can not travel farther either for the same reason. Figure 4.8(b) shows the same
bundle structure except having a unidirectional block, depicted with the dark area
in branch 2. An action potential travels down branch 1 normally but is blocked
in branch 2. The wave transmitted from the left side then is able to propagate
to the right side and splits to travel both upwards and downwards. If when the
wave penetrates the unidirectional block from the bottom, the upper portion of
branch 2 has exited its refractory period and becomes excitable again; the wave
can continue to travel up to where it originally comes from and also to branch
1 again, and this causes a reentry. The presented anatomical reentry example is
simulated with PSpice using the 1-dimensional propagation model. A modifled
transistor circuit of the VLSI implementation of resistors is employed to simulate
the unidirectional block. The circuit setup and the PSpice simulation results are
provided in Appendix B.
162
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .x
y
30        29        28                              3          2 1        0 
0
1
2
29
30
31 impulses stimulate the nodes of x = [0,30] , y = 0 at t = 10 ms.
11 impulses stimulate 
the nodes at x = 30,       
y = [0,10] at t = 335 ms.
Istimulus,1 
Istimulus,2 
Figure 4.9: Circuit for 2-dimensional cardiac propagation.
Circuit Simulation of a Spiral Wave
The circuit setup for 2-dimensional cardiac propagation is illustrated in Figure
4.9. A 31 ? 31 array is shown constructed with each node being composed of a
membrane capacitor and a non-linear ionic current (not depicted here) as shown in
Figure 4.4(b). Resistors are connected in between each two adjacent nodes, which
are also not shown in Figure 4.9. The nodes are identifled with their positions
deflned by the x axis and y axis, with integer points in the intervals x 2 [0;30],
y 2 [0;30]. Equation (4.31) is the mathematical description of the 2-dimensional
propagation system we are trying to model using the circuit in Figure 4.9, and it
is modifled for two space dimensions and rewritten as:
6:21?10?9(@
2V 0
m
@x2 +
@2V 0m
@y2 ) = 10
?7@V
0
m
@t +I
0
ion; (4.31)
where V 0m and I0ion are scaled membrane potential and ionic current. Here we still
use the values used by Bray and Wikswo in [108] for the surface-to-volume ratio
fl = 3000 cm?1, the conductances dx = dy = 0:186 mS=cm, and the membrane ca-
pacitance Cm = 1?F=cm2. Discretizing equation (4.31) by meshing the continuous
spatial domain into grids, this difierential equation can be mapped into the RC
network shown in Figure 4.9 by setting proper values for the resistors, capacitors,
163
and the ionic currents in the flgure; this has been introduced in section 4.3.1. We
set the grids in the x and y direction to be ?x = ?y = 160 ?m. The circuit pa-
rameters are selected under the restriction of equation (4.26), and are summarized
in equation (4.32):
8>
>>>>
><
>>>>
>>:
Rx = Ry = 4:1?104 ?;
C = 0:1 ?F;
I0 = I0ion;
(4.32)
where Rx and Ry are resistances between nodes in the x and y directions, C is
node capacitance, and I0 is the ionic current of each node (refer to Figure 4.4 for
the names of circuit components). Here we still use the Beeler-Reuter model to
formulate the ionic current, provided in equation (2.41)-(2.60).
In order to activate the 2-dimensional RC network and create action potentials,
external stimulus is needed to be applied to the circuit. The stimulating elements
are conflgured for generating a spiral reentry wave on the RC circuit and can be
divided into two groups, as illustrated in Figure 4.9. The flrst group of stimulating
current sources Istimulus;1 are connected to the nodes at y = 0 (x 2 [0;30]), and
they all flre a pulse at time t = 10 ms, which creates a plane action potential wave
that propagates in the increasing y direction. The second group of stimulating
currents Istimulus;2, connected to x = 30 and y 2 [0;10], send an impulse at time
t = 335 ms. This combination triggers a spiral wave as we will see soon. Istimulus;2,
8 mA in magnitude and lasting for 0.1 ms, is set stronger than Istimulus;1, 6.2 mA
of the same duration, in order to successfully re-activate the elements that exit the
refractory period a short while ago.
The simulation results are provided in Figure 4.10. The x axis and the y
axis are the same as in Figure 4.9, and the z axis is the membrane potential.
164
As shown in Figure 4.10(a), at t = 10 ms, the flrst group of stimulus triggers
the components at an edge of the array. After 30 ms, that is at t = 40 ms,
the generated action potential is transmitted to the nodes at y = 16, and this is
depicted in Figure 4.10(b). At t = 300 ms, shown in Figure 4.10(c), the action
potentials have been propagated to every node in the array and the components
close to the origins (at y = 0) are largely repolarized, i.e. back to resting potential.
At t = 335 ms, the second group of stimulus trigger the circuit as illustrated in
Figure 4.10(d), and this situation is known to happen in the heart when a group
of cells flres abnormally [37]. At this moment, many components have not been
fully repolarized from the last action potential. Figure 4.10(e) shows that 10 ms
later, the action potential wave triggered by the second stimulation is propagated
to the neighborhood. Note that a part of the region where the second stimulus
flres is still in the refractory state due to the recent passage of the plane action
potential wave. Therefore, at the beginning the second wave can only propagate in
the decreasing x direction to where the components become excitable. Soon more
components exit the refractory period of the flrst action potential and recover the
ability to be activated, the second wave can enter this region, propagating in the
y direction flrst, and then going in the positive x direction. This is illustrated in
Figure 4.10(f)(g)(h). The wave propagates in the negative y direction as shown
in Figure 4.10(i)(j), when the components where the second wave originates exit
the refractory period and flre a third action potential, giving a reentrant action
potential. This process can continue, forming a rotating spiral wave pattern.
165
t = 10 ms
x y
(a)
Vm
?[V
]
t = 40 ms
x y
(b)
Vm
?[V
]
t = 300 ms
x y
(c)
Vm
?[V
]
t = 335 ms
x y
(d)
Vm
?[V
]
t = 345 ms
x y
(e)
Vm
?[V
]
t = 375 ms
x y
(f)
Vm
?[V
]
166
t = 395 ms
x y
(g)
Vm
?[V
]
t = 415 ms
x y
(h)
Vm
?[V
]
t = 435 ms
x y
(i)
Vm
?[V
]
t = 455 ms
x y
(j)
Vm
?[V
]
Figure 4.10: A spiral wave in 2-dimension cardiac model.
167
Chapter 5
Conclusions and Future Work
5.1 Conclusions
In this dissertation, we have presented analog VLSI circuits which simulate the
electrical activities of the heart in chapter 3 and 4. The VLSI design of the heart
model is based on the Beeler-Reuter model, describing the ionic membrane mech-
anism of ventricular cells, along with the continuous cardiac propagation model
which evolved from cable theory. In order to implement the mathematical de-
scriptions of the heart model in a VLSI circuit, we reformulate the Beeler-Reuter
equations, and discretize the reaction-difiusion equation of the propagation model.
The PSpice simulation of the transistor circuit of a cardiac cell and a segment of
1-dimensional cardiac flber has been carried out, and the circuit simulated action
potential has been discussed.
In chapter 2, we presented a review for the kinetics of the membrane ionic cur-
rents and provided the equations of the Beeler-Reuter model. The Beeler-Reuter
model describes a degree-eight flrst order ordinary difierential system, and for-
mulates the membrane potential, four ionic currents, six gating variables, and a
168
varying calcium concentration. The time-based difierential system is realized using
the integration properties of capacitors and nonlinearities of MOS and NPN tran-
sistors, presented in chapter 3. The circuit simulated action potential of a cardiac
cell showed satisfactory accuracy compared with the action potential calculated
from the original equations in the Beeler-Reuter model.
In chapter 4, we discussed the cardiac propagation model, a partial difierential
equation that describes the variation of the membrane potential in space. The
space-dependent partial difierential equation is realized with a non-linear RC net-
work. The propagation of the cardiac electrical activity has been simulated in
1-dimensional model with a transistor circuit. The important reentry phenomena
and the spiral wave have been simulated successfully with a 2-dimensional model.
It is worth to emphasize that we have proposed a design  ow and a set of
steps for constructing analog VLSI circuits to calculate mathematical equations (in
chapter 3), and this methodology can be applicable to systems other than the heart.
We have introduced equation reformulation schemes to change a mathematical
descriptionforeasiercircuitrealization; thisusestheMatlabCurveFittingToolbox
to transform equations to expressions that can be simulated by simpler circuits, a
parameter scaling procedure to convert numerical ranges of variables to the feasible
working region of circuit currents and voltages, and an initial value shifting scheme
to change the initial conditions of difierential equations to zero. We have designed
a set of transistor circuit blocks of mathematical function units, implemented a
time-based ordinary difierential system (for a cardiac cell), and proposed the circuit
structuretomapreaction-difiusionsystems. Alltheintroducedmethodscontribute
to the fundamentals for a novel technique of obtaining numerical solutions by
analog VLSI circuits. Analog circuits realized on hardware have the potential
169
for fast computation due to their convenient and parallel calculation method for
difierentiation compared to the iterative method used by digital computers, and
this makes it possible to produce application-specialized fast analog processing
elements.
5.2 Future Work
This work can be seen as a flrst step towards the development of a full electrical
conduction system of the entire heart using analog VLSI circuits. To construct a
complete cardiac electrical system with circuits, more work is necessary. Listed
below are several additional aspects that need to be studied.
? Power Characterization
Power e?ciency is an important issue in terms of the satisfactory battery span in
portable electronics. In the case of VLSI implementation of the heart model, the
potential applications, such as portable cardiac processing devices, implantable
medical electronics, etc., all demand low power dissipation. To simulate the whole
heart, millions of spatial nodes may be needed, and, thus, when a full heart model
is realized in silicon and membrane ionic mechanism details are implemented for
each node, power consumption can become critical to necessitate the proper func-
tioning of the whole system. Therefore, the power consumption of the VLSI heart
model needs to be characterized and minimized, and low power techniques (such
as decreasing the biasing currents, simplifying the circuit topologies to reduce the
number of transistors, changing the working region of transistors, lowering the sup-
ply voltages, changing fabrication technology etc.) are desirable to be incorporated
into the circuit design.
170
? Integration of Capacitors
We have shown in chapter 3 that the maximum capacitor (1 ?F) used in the VLSI
design of heart model can be realized by a capacitor 104 times smaller, that is 100
pF. Taking the selected AMI 1.5 ?m IC process as an example, in which the avail-
able capacitor option PiP (poly2 over poly) is about 600 aF=?m2, a capacitance of
100 pF needs a silicon area of about 408?408 ?m2. Thus, the required capacitor
area is still very large, and further study is needed to flnd smaller VLSI elements
that calculate integrals. One way is to improve the presented capacitance ampli-
flcation schemes specifled for the VLSI cardiac circuit (refer to chapter 3 section
3.4.7). As mentioned in section 3.4.7, minimizing the leakage currents of bipolar
devices can reduce the limitation on the capacitor ampliflcation factor. Shrinking
the time scale is another way to avoid large capacitance, and this needs the investi-
gation of high-frequency performance of the VLSI heart circuits of this dissertation.
Also, switched-capacitor circuits [117](pp. 417) could be investigated for a possible
method of magnifying capacitances, in which clock-controlled switches are used to
mediate charging currents. Developing new semiconductor technology for higher
capacitivity, and inventing novel integration devices to replace area-consuming ca-
pacitors could be other directions for solving the problem.
? Enhancement with Heterogeneous Membrane Properties
The presented VLSI heart model only simulates the electrical activities of the ven-
tricular cells. The variation of the action potential generated at difierent portions
of the heart needs to be investigated, and the development of new circuits are
needed to generate these variable action potentials. The current VLSI design can
be extended to include the controlabilities for adjusting the action potential for
171
difierent types of cardiac cells. Due to the great similarities of the action poten-
tials at various kinds of cardiac cells (such as ventricle, atrial, sinus node, Purkinje
flbers, etc.) and the same ionic basis of all the models, it is possible to implement
the action potentials of difierent cells with similar circuits for changing gating
opening/closing rates and channel ionic properties. A control mechanism added
in the circuit implementation may be able to let the same circuit simulate various
types of action potentials. In addition, the added control circuit can also serve as
a handler to tune the system for simulating cardiac response in various conditions
such as drug efiects, emotions regulation, etc.
? ErrorAnalysisofImplementingtheContinuousExcitationPropagationModel
The same as for all the numerical method, the circuit implementation of discretized
reaction-difiusion system sufiers with discretization errors. The spatial discretiza-
tion may introduce errors in the action potentials and in uence the propagation
of the action potential in space, and the wave conduction velocity is a typical and
easily recognizable indicator for coarse discretization [38]. The efiects of spatial
segmentation in the continuous model of the cardiac excitation propagation have
been discussed in previous research works [38][100]. Similar analysis of the seg-
mentation errors introduced by flnite space steps on the circuit implementation
can be carried out in our heart model by taking account the continuous nature
of the time domain in the VLSI heart simulation. Furthermore, as mentioned in
chapter 4 section 4.2.1, the continuous cable model takes the intracellular space as
a medium with continuous difiusive properties and does not consider the discrete
properties of the gap junctions. Therefore, a study on how close the circuit imple-
mented excitation propagation model is to reality is necessary to provide a helpful
guidance for creating more accurate VLSI heart models.
172
? Real 3-dimensional - Adaption to Heart Natural Topology
Though the heart can be simplifled by 1-dimensional (for the specialized exci-
tatory and conductive system from the sinus node to the Purkinje flbers) and
2-dimensional (for the atrial and ventricular walls) models, the heart has a 3-
dimensional topology, with curved and rotated muscle flbers and non-zero thickness
of the atrial and ventricular walls, and, thus, 3-dimensional models are unavoidable
to describe more accurately the complicated structure of the heart. However, in the
VLSI implementation of the heart model, 3-dimension means huge entangled wire
connections and densely occupied circuits, because the 3-dimensional network has
to be projected to 2-dimension surface due to the 2-dimension nature of the cur-
rent VLSI technology; this provokes the need for new semiconductor techniques for
implementing heart circuits. The concepts for 3-D VLSI [118] and  exible (bend-
able) transistors [119] are possible directions to implement real 3-dimensional heart
models with circuits. The study of the detailed anatomic structure of the heart
and the orientation of the cardiac flbers is necessary to incorporate into the VLSI
implementation. Irregular meshing for simulating the complex geometry of the
cardiac muscle is also desirable to be investigated.
5.2.1 More - Application Perspectives
Up to now, we have been concentrating on the simulations of the cellular ionic
mechanisms and the resulting electrophysiological activities in the normal heart.
The eventual purpose of modeling the cardiac electrical behavior is to understand
thecardiac physicalphenomena and themechanismsofheartfailure, andapply this
knowledge to biomedical research to contribute to the development of therapeutic
strategies, drugs and medical devices. Therefore, the capability of simulating ab-
173
normal activities in the heart and the inclusion of other interactions with the heart
electrical behaviors are imperative to make the VLSI heart model more thorough
and, hence, beneflcial for real life cardiology. In the following, we give an overview
on what else can be included in the circuit cardiac simulation.
? Models for Diseased Cardiac Tissues
Acute ischemia in cardiac muscles can critically impair the mechanical and electri-
cal behavior of the heart and facilitate the development of reentry waves, cardiac
arrhythmias, and ventricular flbrillation [26][120]. The enhancement of the model
to be able to handle ischemia cardiac tissues necessitates the simulation of the
abnormal electrophysiological properties in the heart caused by ischemia. Sophis-
ticated models that include the formulation of K+ concentrations, ATP-sensitive
K+ channels, and efiects of pH are needed to simulate this kind of diseased cardiac
cells [121]. The relative previous work can be found in [120][122] for the simulation
of ischemia in the heart. It is even more exciting to model the mutated cardiac cells
and inherited heart disorders. Such cell models provide us insight into the genetic
bases of arrhythmias and heart diseases. [123] is an example of using simulation
to study the action potential changes in mutated cardiac tissue.
? Arrhythmias and Their Underlying Mechanisms
In this work we only simulated reentry waves, which is the most common mech-
anism of arrhythmia and plays a crucial role in the initiation of ventricular flb-
rillation, the most severe arrhythmia. The simulation of the progression and sus-
tainment of ventricular flbrillation is not performed in this work but is important
to gain knowledge of the dynamics of ventricular flbrillation and invent e?cient
174
treatment. Besides, to enrich our understanding of the arrhythmia caused by reen-
try (i.e. abnormal impulse conduction), more precise cardiac conduction models
may be favorable to extract more information from the simulation. One aspect
of improving the model is to include the discrete properties of gap junctions and
their dynamic resistances. Gap junctions are clusters of intercellular channels,
which exhibit several conductance states [124], in which the probability of a chan-
nel being in the open state is determined by the transjunction voltage. Hence, the
conductivity of gap junctions changes during the transmission process of electrical
signals. [99] is a previous work that studied the in uence of varying gap junction
resistance on the electrical propagation in cardiac tissues.
Furthermore, there are two other basic mechanisms causing arrhythmias. Un-
like reentry-caused arrhythmias, which result from abnormal impulse conduction,
these two are both related to abnormal impulse generation, which arises from
abnormal automaticity or triggered activity [121]. Automaticity refers to the abil-
ity of the heart to initiate its own beat, and it allows the heart to beat even
when it is removed from the body [26]. Suppressed automaticity can lead to si-
nus node dysfunction, and enhanced automaticity can result in arrhythmias [125].
Triggered activity is the spontaneous multiple depolarization triggered by early
afterdepolarizaiton (EAD) and delayed afterdepolarizaiton (DAD). Previous works
that simulate abnormal automaticity and triggered activity have been reported
[126]-[128].
? Medication and Treatment
It is important to include medication and other methods of treatment in the sim-
ulation model in order to study the side efiects of medicines and validate new
medical devices. Therefore, to include the drug efiects on the ionic channels and
175
to model implantable electronic devices are necessary for research of investigating
cures for heart diseases. The anti-arrhythmic drugs are divided into four classes
according to the Vaughan-Williams classiflcation [129], and most of the drugs work
by regulating the ionic channels and adjusting the depolarization or repolarization
process. There have been literatures that formulate the equations for drug efiects
[130]-[132] and they can be easily added to the circuit simulation by the VLSI
implementation schemes presented in chapter 3. Implantable devices for heart
diseases include artiflcial pacemakers and implantable cardioverter/deflbrillators
(ICD). Artiflcial pacemakers sense the heartbeat and send electrical charge when
the heartbeat is slower than a threshold for triggering more heartbeat. ICD works
in a similar way, but it \actions" when detecting heart rate being too fast and
delivers an electrical shock to return the rhythm to normal. Therefore, the im-
plantable devices can also modify the heart electrical behavior and need to be
included in appropriate simulations.
? Electrocardiogram
An electrocardiogram (ECG) is a graphic produced by recording the electrical
voltageonthesurfaceofthebodygeneratedbythesmallcardiaccurrentscausedby
an impulse that spreads into the tissues surrounding the heart [2]. Since abnormal
electrical activity of the heart alters the shapes of the waves in an ECG, an ECG is
the prime clinical tool for the diagnosis of cardiovascular diseases. The capability
of constructing the ECG associated with the simulated arrhythmias would enable
us to link between the electrical activities inside the heart and the measurement on
the surface of the body and would assist in better understanding of the mechanisms
underlying the heart diseases observed clinically. Previous works such as [133] and
[134] have proposed methods of modeling ECG in simulations.
176
? Coupling with Nervous Regulation and Mechanical Systems
Sofarwehavebeendiscussing theelectrical activityofthe heart alone. It isobvious
that the heart, being embedded in the body, can be also in uenced by other parts
of the body. The automaticity of the heart can be regulated by extrinsic factors
such as temperature, oxygen tension in the blood, nervous control, hormones, etc.
There have been researches carried out to formulate mathematical models for the
control of heart rate by autonomic nervous system [135]-[137]. In addition, the
heart rate can be also increased by the mechanical stretch of the right atrial wall
[2]. The mechanical behaviors are tightly coupled with the electrical system in the
heart. On one hand, the electrical activity causes the mechanical contraction of the
heart; on the other hand, the mechanics in uences the electrical behavior by the
mechano-electric feedback. The study of the electrical and mechanical components
together allows us to expand the exploration from the cardiac electrophysiology
to the full functioning of the heart, and the cardiovascular physiology and the
circulatory system. The mathematical models that describe the cardiac mechanics
have been studied by researchers and the simulation of the electro-mechanical
coupling have been performed. The related literatures are [5][138][139].
177
Appendix A
Transistor Circuits of
Beeler-Reuter Model
Implementation
In Appendix A, we provide the complete transistor circuits for implementing the
Beeler-Reuter model. The circuit schematics are divided into six categories, as
listed in Table A.1: Vm, the top level circuit; IK1, Ix1, INa, and Is, the circuits for
creating four ionic currents; and Miscellaneous circuits. The block diagrams of
the circuits and their simulation results have been presented in chapter 3, section
3.4.
The schematics of the cardiac cell VLSI design are shown hierarchically in a top-
down manner, i.e. an upper level of circuit is presented flrst, and its sub-circuits
are provided later. The legends related to the hierarchy are illustrated in Figure
A.1. Figure A.1(a) is a hierarchic module, whose module name and the schematics
page number of its transistor implementation are provided below the icon. In a
schematics of a module, the input and output terminals, bearing the same names
178
Vin     Vout
V-Bufferpp. Sch-27 Vin Vout VbiasVbias
(a) (b) (c)
Figure A.1: Legends for hierarchical schematics.
as shown in its the module block, are represented by the symbols shown in Figure
A.1(b). To help organize the schematics circuits, on-page connectors, depicted in
Figure A.1(c), are used to make wire connections without drawing wires among
the nodes that are distributed far away to each other in the same schematic page.
The schematics? contents are summarized in Table A.1. In each schematic, the
schematic?s page number is shown in the middle below the circuit, and a brief
description is provided in the right bottom corner of each page.
179
Table A.1: List of schematics
Category Description Schematics Page
Vm Top level of circuit for Vm Sch-1
Top level diagram for IK1 Sch-2
IK1 Sigmoid circuit 1 for IK1 Sch-3
Sigmoid circuit 2 for IK1 Sch-4
Top level diagram for Ix1 Sch-5
Ix1 fix1, and flx1 Sch-6, 7
Ix1 Sch-8
Multiplier for IK1 Sch-9
Top level diagram for INa Sch-10
INa fim, flm, fih, flh, fij, and flj Sch-11 to Sch-16
Multi-input multiplier for INa Sch-17
2-input multiplier and transconductor for INa Sch-18
Top level diagram for Is Sch-19
fid, fld, fif, and flf Sch-20 to Sch-23
Is [Ca]i and Es Sch-24
Transconductors for Is Sch-25
3-input multiplier for Is Sch-26
NPN, Vbias, V-Bufier, Transconductor G1 Sch-27
Transconductors G2, G3, G4, and G5 Sch-28
Miscellaneous Transconductor G6, and pool circuit 1 & 2 Sch-29
Tanno multiplier 1 Sch-30
Tanno multiplier 2 and current duplicator Sch-31
Resistor Sch-32
180
Pa
ge
 S
ch
-1
Ca
teg
ory
De
scr
ipti
on
Vm To
p le
vel
 cir
cui
t fo
r V
m
Eq
uat
ion
(3.
99)
181
Pa
ge
 S
ch
-2
Ca
teg
ory
De
scr
ipti
on
IK1 To
p le
vel
 cir
cui
t fo
r IK
1
Eq
uat
ion
(3.
100
)
182
Pa
ge
 S
ch
-3
Ca
teg
ory
De
scr
ipti
on
IK1 Sig
mo
id c
ircu
it 1
 for
 IK1
Eq
uat
ion
(3.1
00)
 2n
d t
erm
183
Pa
ge
 S
ch
-4
Ca
teg
ory
De
scr
ipti
on
IK1 Sig
mo
id c
ircu
it 2
 for
 IK1
Eq
uat
ion
(3.
100
) 3r
d t
erm
184
Pa
ge
 S
ch
-5
Ca
teg
ory
De
scr
ipti
on
Ix1 To
p le
vel
 cir
cui
t  fo
r Ix
1
Eq
uat
ion
(3.
103
), (
3.1
17)
185
Pa
ge
 S
ch
-6
Cat
ego
ry
De
scr
iptio
n
Ix1 ?x1
Equ
atio
n
(3.1
13)
186
Pa
ge
 S
ch
-7
Cat
ego
ry
De
scr
iptio
n
Ix1 ?x1
Equ
atio
n
(3.1
14)
187
Pa
ge
 S
ch
-8
Ca
teg
ory
De
scr
ipti
on
Ix1 Ix1
Eq
uat
ion
(3.1
16)
188
Pa
ge
 S
ch
-9
Cat
ego
ry
De
scr
iptio
n
Ix1 Mu
ltip
lier
 for
 x1
 an
d Ix
1
Equ
atio
n
(3.1
17)
189
Pa
ge
 S
ch
-1
0
Ca
teg
ory
De
scr
ipti
on
INa Top
 lev
el c
ircu
it fo
r IN
a
Eq
uat
ion
(3.1
03)
, (3
.11
8)
190
Pa
ge
 S
ch
-1
1
Ca
teg
ory
De
scr
ipti
on
INa ?m
Eq
uat
ion
?m
 in
 Ta
ble
 3.5
191
Pa
ge
 S
ch
-1
2
Ca
teg
ory
De
scr
ipti
on
INa ?m
Eq
ua
tion
?m
 in
 Ta
ble
 3.
5
192
Pa
ge
 S
ch
-1
3
Ca
teg
ory
De
scr
ipti
on
INa ?h
Eq
ua
tion
?h
 in
 Ta
ble
 3.5
193
Pa
ge
 S
ch
-1
4
Ca
teg
ory
De
scr
ipti
on
INa ?h
Eq
uat
ion
?h 
in T
abl
e 3
.5
194
Pa
ge
 S
ch
-1
5
Ca
teg
ory
De
scr
ipti
on
INa ?j
Eq
uat
ion
?j i
n T
abl
e 3
.5
195
Pa
ge
 S
ch
-1
6
Ca
teg
ory
De
scr
ipti
on
INa ?j
Eq
uat
ion
?j i
n T
abl
e 3
.5
196
Pa
ge
 S
ch
-1
7
Ca
teg
ory
De
scr
ipti
on
INa Mu
ltip
lier
 for
 INa
Eq
uat
ion
(3.1
18)
, Fi
gur
e 3
.30
197
Pa
ge
 S
ch
-1
8
Ca
teg
ory
De
scr
ipti
on
INa (a)
 2-
Inp
ut m
ulti
plie
r fo
r In
a
(b)
 V-
I co
nve
rte
r fo
r IN
a
Eq
uat
ion
(3.
118
)
198
Pa
ge
 S
ch
-1
9
Ca
teg
ory
De
scr
ipti
on
Is Top
 lev
el c
ircu
it fo
r Is
Eq
uat
ion
(3.1
28)
199
Pa
ge
 S
ch
-2
0
Cat
ego
ry
De
scr
iptio
n
Is ?d
Equ
atio
n
(3.1
23)
200
Pa
ge
 S
ch
-2
1
Cat
ego
ry
De
scr
iptio
n
Is ?d
Equ
atio
n
(3.1
24)
201
Pa
ge
 S
ch
-2
2
Ca
teg
ory
De
scr
ipti
on
Is ?f
Equ
atio
n
(3.1
26)
202
Pa
ge
 S
ch
-2
3
Cat
ego
ry
Des
crip
tion
Is ?f
Equ
atio
n
(3.1
27)
203
Pa
ge
 S
ch
-2
4
Ca
teg
ory
De
scr
ipti
on
Is [Ca
]ia
nd 
Es
Eq
uat
ion
(3.
129
), (
3.1
30)
204
Pa
ge
 S
ch
-2
5
Ca
teg
ory
De
scr
ipti
on
Is Tra
nsc
ond
uct
ors
 
Eq
uat
ion
(3.
128
), (
3.1
30)
205
Pa
ge
 S
ch
-2
6
Ca
teg
ory
De
scr
ipt
ion
Is Mu
ltip
lier
 fo
r Is
Eq
ua
tion
(3.
12
8)
206
Pa
ge
 S
ch
-2
7
Cat
ego
ry
De
scr
iptio
n
Mis
cel
lan
eou
s
Vbi
as1
, Vb
ias2
, V-
buf
fer,
 NP
N, G
1
207
Pa
ge
 S
ch
-2
8
Cat
ego
ry
De
scr
iptio
n
Mis
cel
lan
eou
s
G2
, G
3, G
4, G
5.
208
Pa
ge
 S
ch
-2
9
Ca
teg
ory
De
scr
ipti
on
Mis
cel
lan
eou
s
G6
, po
ol c
ircu
it 1
 &  
2.
209
Pa
ge
 S
ch
-3
0
Ca
teg
ory
De
scr
ipt
ion
Mi
sc
ell
an
eo
us
Mu
ltip
lier
 1 
210
Pa
ge
 S
ch
-3
1
Cat
ego
ry
De
scr
iptio
n
Mis
cel
lan
eou
s
Mu
ltip
lier
 2, 
cur
ren
t du
plic
ato
r. 
211
Pa
ge
 S
ch
-3
2
Ca
teg
ory
De
scr
ipti
on
Mis
cel
lan
eou
s
Re
sist
or
212
Appendix B
PSpice Simulation of Anatomical
Reentry Using Unidirectional
Block
In Appendix B, we present a 1-dimensional circuit that simulates the anatomical
reentry phenomena presented in section 4.3.3. The triangular structure of the
reentry cardiac tissues, shown in Figure 4.8, is to be modeled by a circuit ring,
which is similar to that in Figure 4.5, except the flrst node and the last node is
connected with a special resistor that models a unidirectional block. The simplifled
diagram of the circuit is illustrated in Figure B.1, where the ring represents a 1-
dimensional RC circuit (refer to Figure 4.5) composed of 120 nodes, starting from
node #1 at x = 0;y = 1, with the node ID number arranged in the increasing
order in the direction of anti-clockwise, and #1, #31, #61, and #91 are the four
nodes that are on the x or y axis. There is a unidirectional block between node
#1 and #120, which blocks the electrical propagation from node #1 to #120 and
only allows the passage of signals in the opposite direction. We have presented
213
x
y
Stimulate
node #1
#1 [0,1]
#31 [1,0]#61 [0,-1]
#91 [-1,0]
Action potential 
propagation direction
Uni-directional 
block
Figure B.1: A circuit ring composed of 120 node for modeling anatomical reentry.
in detail in section 4.3.2 the construction of an RC circuit to model the reaction-
difiusion system of the cardiac propagation, and all circuit parameters in the ring
circuit follow section 4.3.2, and, thus, we will not repeat the circuit conflguration
here. In the following discussion, we focus on the introduction of implementing
the unidirectional block and show the simulation results of the anatomical reentry.
Constant resistors can be realized with VLSI devices using the difierential pair
circuit shown in Figure 3.5(b), as introduced in section 3.3.1. On the top of
Figure B.2, (a1) depicts a simplifled version of the transistor-implemented resistor
presented in section 3.3.1, where R+ and R? are the two terminals of the equivalent
resistor. The block is the difierential pair circuit, whose voltage inputs V+ and V?
and current outputs I+ and I? are connected like in Figure B.2 (and Figure 3.5(b))
for a normal resistor. The relation of the resistor currents vs. the difierence of the
input voltages of the resistor (we call it a symmetric resistor) is shown in Figure
B.2(a2).
There are two candidate circuits for implementing the unidirectional block,
and both are based on a difierential circuit. The flrst candidate is shown in Figure
B.2(b1), in which the terminal I? does not connect to R? and hence there is no
214
Symmetric
Resistor
IR+ 
VR+ -VR-I+ I-
V-V+
R+ R-
Non-symmetric
Resistor IR+
VR+ -VR-
I+ I-
V-V+
R+ R-
Half-symmetric
Resistor
VR+ -VR-
I+ I-
V-V+
R+ R-
IR+
IR-
VR+ -VR-
VR+ -VR-
IR-
IR-
VR+ -VR-
NMOS
I-mirror
PMOS
I-mirror
(a1) (a2)
(b1) (b2)
(c1) (c2)
Figure B.2: Implementation of bidirectional and unidirectional resistors using a
modifled difierential pair circuit.
current allowed at the R? terminal, while the R+ terminal remains the same as for
(a1). Since the I-V properties of the resistor implementation in (b1) is difierent in
two terminals, we call it a non-symmetric resistor. This non-symmetric resistor
does not in uence the external circuit connected with R?, but afiects the circuit
connected to R+. Therefore, a signal can pass from R? to R+, but not from R+ to
R?. Figure B.2(c1) shows the second circuit that realizes the unidirectional block,
where the currents sourced from I+ and I? are \cut" to leave a half of the I-V
branch by adding a NMOS and a PMOS current mirror; and this is shown in its
I-V curves in Figure B.2(c2). Due to the \halved" I-V properties of this resistor,
we call the resister a half-symmetric resistor. A half-symmetric resistor operates
like a normal resistor when V+ > V?, and when V+ < V? it shows no conductivity.
If a non-symmetric resistor is employed in our ring circuit, it leads to large
re-polarizing current (caused by big voltage difierence across the resistor) that
rapidly discharges the membrane and distorts the action potential when the cardiac
node at one side of the resistor is activated and the node at the other side can
not be due to the isolation efiect of the non-symmetric resistor. Therefore, in
215
the simulation of the anatomic reentry, we adopt the half-symmetric resistor to
implement the unidirectional block. Since an excited membrane potential is higher
than the resting potential, the properties of the half-symmetric resistor necessitate
the propagation of the active potential in one direction but not the other. An
external stimulus is applied to excite the ring of the 1-dimensional cardiac tissue.
It sends a square pulse of current to node #1 at time t = 10 ms, with 2 mA in
the amplitude, and 0.1 ms of duration. Ideal components are used to construct
all the nodes for saving the longitude simulation time. The simulation results are
provided in Figure B.3. At t = 11 ms, shown in Figure B.3(a), the flrst node is just
activated by the stimulus. In (b), at t = 180 ms, the action potential is propagated
over about a half of the circle. Note that the propagation is anti-clockwise due to
the existence of the unidirectional block between #1 and #120. At t = 380 ms,
all nodes on the ring have been activated at lease once, and the action potential
is transmitted back to node #1 and creates a reentry wave around there; this is
shown in Figure B.3(c). Figure B.3(d) and (e) show that the reentry wave is passed
further down to other nodes. The big notch on the voltage potential contour at
the positive x and y quarter is due to the shorter duration of the newly reentered
action potential (at nodes number ? 1) compared to the flrst action potential of
the same node. This results from the not-fully-recovered refractory state of the
nodes when the second action potential is activated. This can be seen clearly in
Figure B.4, which we will discuss shortly. In Figure B.3(f), at t = 780 ms, a
third circulation of the action potential starts. The reentry wave will turn around
and around and never stop, and it is responsible for many clinical arrhythmias
(disturbances of cardiac rhythm) [26].
Figure B.4 illustrates the membrane potential of node #1 changed by time.
216
t = 11 ms t = 180 ms
t = 380 ms t = 460 ms
t = 560 ms t = 780 ms
(a) (b)
(c) (d)
(e) (f)
Figure B.3: Reentry wave caused by unidirectional block in 1-dimensional model.
217
Figure B.4: The time-course of the membrane potential of node #1.
A full action potential is generated at t = 10 ms when the stimulus sends a
current impulse to the node. The second action potential, occurring at about
t = 360 ms, is much smaller than the flrst one due to its activation being too
close to the repolarization stage of the flrst action potential, and consequently its
magnitude is degraded by the refractory period of the flrst action potential. The
third action potential is a full action potential again, and it happens at t = 730 ms.
The excitation of node #1 will continue on and on, and the shape of each action
potential depends on how well the node is recovered from the previous action
potential.
218
BIBLIOGRAPHY
219
BIBLIOGRAPHY
[1] \Causes of Deaths", http : ==www:deardeath:com=causes of death:htm.
[2] A. C. Guyton, Textbook of Medical Physiology, W. B. Saunders Company,
Harcourt Brace Jovanovich, Inc., 8th edition, 1991.
[3] \Heart", http : ==www:answers:com=topic=heart?4.
[4] U.S. Department of Health, Education, and Welfare, A Handbook of Heart
Terms, U.S. Government Printing O?ce, May 1964, pp. 42.
[5] F. B. Sachse, Computational Cardiology: Modeling of Anatomy, Electrophys-
iology, and Mechanics, Springer Berlin/Heidelberg, 2004.
[6] J. W. Cain and D. G. Schaefier, \Two-Term Asymptotic Approximation of a
Cardiac Restitution Curve", SIAM Review, Vol. 48, No. 3, March 2006, pp.
537-546.
[7] A. L. Hodgkin and A. F. Huxley, \A Quantitative Description of Membrane
Current and its Application to Conduction and Excitation in Nerve", Journal
of Physiology, Vol. 117, No. 4, August 1952, pp. 500-544.
[8] R. A. FitzHugh, \Impulses And Physiological States In Theoretical Models Of
Nerve Membrane", Biophysical Journal, Vol. 1, No. 4, July 1961, pp. 445-466.
[9] J. Nagumo, S. Animoto, and S. Yoshizawa, \An Active Pulse Transmission
Line Simulating Nerve Axon", Proceedings of the Institute of Radio Engineers,
Vol. 50, No. 10, October 1962, pp. 2061-2070.
[10] B. van der Pol, \On Relaxation Oscillations", Philosophical Magazine, Vol.
2, No. 11, July 1926, pp. 978-992.
[11] H. Glaze, \Mathematical Models for Heart Rate Responses to Vagel Nerve
Simulation", Ph.D. Dissertation, Electrical Engineering Department, Stan-
ford University, 1971.
220
[12] R. H. Clayton, \Computational Models of Normal and Abnormal Action
Potential Propagation in Cardiac Tissue: Linking Experimental and Clinical
Cardiology", Physiological Measurement, Vol. 22, No. 3, August 2001, pp.
15-34.
[13] D. Noble, \A Modiflcation of the Hodgkin-Huxley Equations Applicable to
Purkinje Fibre Action and Pace-maker Potentials", Journal of Physiology,
Vol. 160, No. 2, February 1962, pp. 317-352.
[14] R. E. McAllister, D. Noble, and R. W. Tsien, \Reconstruction of the Electrical
Activity of Cardiac Purkinje Fibres", Journal of Physiology, Vol. 251, No. 1,
September 1975, pp. 1-59.
[15] G. W. Beeler and H. Reuter, \Reconstruction of the Action Potential of
Ventricular Myocardial Fibres", Journal of Physiology, Vol. 268, No. 1, June
1977, pp. 177-210.
[16] L. Ebihara and E. A. Johnson, \Fast Sodium Current in Cardiac Muscle:
A Quantitative Description", Biophysical Journal, Vol. 32, No. 2, November
1980, pp. 779-790.
[17] C. Luo and Y. Rudy, \A Model of the Ventricular Cardiac Action Potential -
Depolarisation, Repolarisation and Their Interaction", Circulation Research,
Vol. 68, No. 6, June 1991, pp. 1501-1526.
[18] C. Luo and Y. Rudy, \A Dynamic Model of the Cardiac Ventricular Ac-
tion Potential - Simulations of Ionic Currents and Concentration Changes",
Circulation Research, Vol. 74, No. 6, June 1994, pp. 1071-1097.
[19] D. Noble and Y. Rudy, \Models of Cardiac Ventricular Action Potentials:
Iterative Interaction Between Experiment and Simulation", The Royal Society
of London Philosophical Transactions. Series A,Vol.359, No.1783, June2001,
pp. 1127-1142.
[20] O. V. Aslanidi, V. N. Biktashev, M. Chen, R. H. Clayton, A. V. Holden, J.
V. Tucker, and H. Zhang, \Computational Biology of Cardiac Arrhythmias:
from Basic Science to Clinical Application", Proceedings of UK e-Science All
Hands Meeting, Nottingham UK, September 2003, pp. 889-896.
[21] L. F. Pavarino and P. C. Franzone, \Parallel Solution of Cardiac Reaction-
Difiusion Models", Proceedings of the 15th International Conference on Do-
main Decomposition Methods, Berlin, Germany, July 2003, pp. 669-776.
[22] J. B. Pormann, C. S. Henriquez, J. A. Board (Jr), D. J. Rose, D. M. Harrild,
and A. P Henriquez, \Computer Simulations of Cardiac Electrophysiology",
221
Proceedings of the 2000 ACM/IEEE Conference on Supercomputing, Dallas
Texas USA, November 2000, pp. 24.
[23] R. Sarpeshkar, \Analog Versus Digital: Extrapolating from Electronics to
Neurobiology", Neural Computation, Vol. 10, No. 7, October 1998, pp. 1601-
1638.
[24] M. R. Boyett, A. Clough, J. Dekanski, and A. V. Holden, \Modelling Cardiac
Excitation and Excitability", Computational Biology of the Heart, John Wiley
& Sons, chapter 1, March 1997, pp. 1-46.
[25] S. P. Ellner and J. Guckenheimer, \Membrane Channels and Action Poten-
tials", Dynamic Models in Biology, Princeton University Press, chapter 3,
March 2006, pp. 71-106.
[26] R. M. Berne, M. N. Levy, B. M. Koeppen, and B. A. Stanton, Physiology,
Mosby, Inc., 4th edition, 1998.
[27] \Ventricular Action Potential", http : ==www:answer:com=topic=ventricular
?action?potential.
[28] F. Bezanilla, \The Voltage Sensor in Voltage-Dependent Ion Channels", Phys-
iological Reviews, Vol. 80, No. 2, April 2000, pp. 555-592.
[29] E. C. Cooper and L. Y. Jan, \Ion Channel Genes and Human Neurological
Disease: Recent Progress, Prospects, and Challenges", Proceedings of the
National Academy of Sciences, Vol. 96, No. 9, April 1999, pp. 4759-4766.
[30] \Ion Channel", http : ==en:wikipedia:org=wiki=Ion channel.
[31] J. Malmivuo and R. Plonsey, \The Heart", Bioelectromagnetism, Oxford
University Press, New York, 1995, pp. 119-130.
[32] D. DiFrancesco and D. Noble, \A Model of Cardiac Electrical Activity In-
corporating Ionic Pumps and Concentration Changes", The Royal Society
of London Philosophical Transactions. Series B, Vol. 307, No. 1133, January
1985, pp. 353-398.
[33] R. J. Ramirez, S. Nattel, and M. Courtemanche, \Mathematical Analysis
of Canine Atrial Action Potentials: Rate, Regional Factors, and Electrical
Remodeling", American Journal of Physiology - Heart and Circulatory Phys-
iology, Vol. 279, No. 4, October 2000, pp. 1767-1785.
[34] H. Zhang, A. V. Holden, I. Kodama, H. Honjo, M. Lei, T. Varghese, and
M. R. Boyett, \Mathematical Models of Action Potentials in the Periphery
and Center of the Rabbit Sinoatrial Node", American Journal of Physiology
- Heart and Circulatory Physiology, Vol. 279, No. 1, July 2000, pp. 397-421.
222
[35] K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, and A. V. Panfllov, \A Model
for Human Ventricular Tissue", American Journal of Physiology - Heart and
Circulatory Physiology, Vol. 286, No. 4, 2004, pp. 1573-1589.
[36] C. S. Henriquez and A. A. Papazoglou, \Using Computer Models to Under-
stand the Roles of Tissue Structure and Membrane Dynamics in Arrhythmo-
genesis", Proceedings of the IEEE, Vol. 84, No. 3, March 1996, pp. 334-354.
[37] N. F. Otani, \Computer Modeling in Cardiac Electrophysiology", Journal of
Computational Physics, Vol. 161, No. 1, June 2000, pp. 21-34.
[38] O. Blanc, \A Computer Model of Human Atrial Arrhythmial", http :
==ltswww:epfl:ch=pub files=blanc=thesis 2537:pdf, 2002.
[39] H. Lodish, A. Berk, S. L. Zipursky, P. Matsudaira, D. Baltimore, and J.
Darnell, \Transport Across Cell Membranes", Molecular Cell Biology W.H.
Freeman & Company; 4th edition, October 1999.
[40] A. Hernndez-Cruz, M. Daz-Muoz, M. Gmez-Chavarn, R. Canedo-Merino, D.
A. Protti, A. L. Escobar, J. Sierralta, and B. A. Surez-Isla, \Properties of the
Ryanodine-sensitive Release Channels that Underlie Cafieine-induced Ca2+
Mobilization from Intracellular Stores in Mammalian Sympathetic Neurons",
European Journal of Neuroscience , Vol. 7, No. 8, August 1995, pp. 1684-1699.
[41] M. R. Guevara, \Dynamics of excitable cells", Nonlinear Dynamics in Physi-
ology and Medicine, Springer-Verlag, New York, chapter 4, January 2003, pp.
87-121.
[42] A. L. Hodgkin and A. F. Huxley, \Currents Carried by Sodium and Potas-
sium Ions through the Membrane of the Giant Axon of Loligo", Journal of
Physiology, Vol. 116, No. 4, April 1952, pp. 449-472.
[43] M. J. Ackermann and D. E. Clapham, \Ion Channels - Basic Science and
Clinical Disease", New England Journal of Medicine, Vol. 336, No. 22, May
1997, pp. 1562-1567.
[44] S. S. Demir, J. W. Clark, C. R. Murphey, and W. R. Giles, \A Mathematical
Model of a Rabbit Sinoatrial Node Cell", American Journal of Physiology,
Vol. 266, No. 3, March 1994, pp. 832-852.
[45] L. Priebe and D. J. Beuckelmann, \Simulation Study of Cellular Electric
Properties in Heart Failure", Circulation Research, Vol. 82, No. 11, June
1998, pp. 1206-1223.
[46] A. F. Strassberg and L. J. DeFelice, \Limitations of the Hodgkin-Huxley
Formalism: Efiects of Single Channel Kinetics on Transmembrane Voltage
Dynamics", Neural Computation, Vol. 5, No. 6, November 1993, pp. 843-855.
223
[47] H. Reuter, \Divalent Cations as Charge Carriers in Excitable Membranes",
Progress in biophysics & Molecular Biology, Vol.26, 1973, pp. 1-43.
[48] H. Reuter and H. Scholz, \A Study of the Ion Selectivity and the Kinetic
Properties of the Calcium-Dependent Slow Inward Current in Mammalian
Cardiac Muscle", Journal of Physiology, Vol.264, No. 1, January 1977, pp.
17-47.
[49] J. B. Bassingthwaighte and H. Reuter, \Calcium Movements and Excitation
Contraction Coupling in Cardiac Cells", Electrical Phenomena in the Heart,
Academic Press New York NY, 1972, pp. 353-395.
[50] \Numerical Analysis", http : ==en:wikipedia:org=wiki=Numerical analysis.
[51] K. Ravindran, K. Ramarao, E. Vidal, and M. Ismail, \Compact Low Voltage
Four Quadrant CMOS Current Multiplier", Electronics Letter, Vol. 37, No.
24, November 2001, pp. 1428-1429.
[52] O. Oliaei and P. Loumeau, \Four-quadrant class AB CMOS current multi-
plier", Electronics Letter, Vol. 32, No. 25, December 1996, pp. 2327-2329.
[53] R. J. Wiegerink, \A CMOS Four-Quadrant Analog Current Multiplier", 1991
IEEE International Symposium on Circuits and Systems, Vol. 4, Singapore,
June 1991, pp. 2244-2247.
[54] K. Tanno, O. Ishizuka, and Z. Tang, \Four-quadrant CMOS Current-mode
Multiplier Independent of Device Parameters", Circuits and Systems II: Ana-
log and Digital Signal Processing, IEEE Transactions on, Vol. 47, No. 5, May
2000, pp. 473-477.
[55] B. Maundy and M. Maini, \A Comparison of Three Multipliers Based on the
V 2gs Technique for Low-voltage Applications", IEEE Transactions on Circuits
and Systems I: Fundamental Theory and Applications, Vol. 50, No. 7, July
2003, pp. 937-940.
[56] G. Han and E. Sanchez-Sinencio, \CMOS Transconductance Multipliers: a
Tutorial", Circuits and Systems II: Analog and Digital Signal Processing,
IEEE Transactions on, Vol. 45, No. 12, December 1998, pp. 1550-1563.
[57] S. Liu and C. Chang, \CMOS Analog Divider and Four-Quadrant Multiplier
Using Pool Circuits", Solid-State Circuits, IEEE Journal of, Vol. 30, No. 9,
September 1995, pp. 1025-1029.
[58] W. Gai, H. Chen, and E. Seevinck, \Quadratic-Translinear CMOS Multiplier-
Divider Circuit", Electronics Letters, Vol. 33, No. 10, May 1997, pp. 860-861.
224
[59] A. Cichocki and R. Unbehauen, \New Analogue Four-Quadrant Voltage Di-
vider", Electronics Letters, Vol. 26, No. 19, September 1990, pp. 1580-1582.
[60] A. Motamed, C. Hwang, and M. Ismail, \CMOS Exponential Current-to-
Voltage Converter", Electronics Letters, Vol. 33, No. 12, June 1997, pp. 998-
1000.
[61] Q. Duong, T. K. Nguyen, and S. Lee, \CMOS Exponential Current-to-Voltage
Circuit Based on Newly Proposed Approximation Method", Proceedings of
the 2004 International Symposium on Circuits and Systems, Vol. 2, May 2004,
pp. 865-868.
[62] \Curve Fitting Toolbox", http : ==www:mathworks:com=access=helpdesk=
help=toolbox=curvefit=index:html?=access=helpdesk=help=toolbox=curvefit=
cftool:html.
[63] L. Wang, Y. Jiang, and R. W. Newcomb, \A Current Based VLSI Degree-Two
Chaos Generator", Chaos in Circuits and Systems, World Scientiflc Series on
Nonlinear Science, Series B, Vol. 11, 2002.
[64] C. Wang and J. Wang, \Design of Linear Transimpedance Ampliflers", Pro-
ceedings. 4th International Conference on ASIC 2001, Shanghai, China, Oc-
tober 2001, pp. 232-235.
[65] I. Mucha, \Low-output-impedance CMOS Voltage Bufier", Electronics Let-
ters, Vol. 28, No. 22, October, 1992.
[66] M. Neag and O. McCarthy, \Low-output-impedance CMOS Voltage Bufier",
Semiconductor Conference, 1998. CAS ?98 Proceedings. 1998 International,
Vol. 1, Sinaia, Romania, October, 1998, pp. 175-180.
[67] S. W. Tsay and R. W. Newcomb, \A Neuro-type Pool Arithmetic Unit",
Proceedings of IEEE ISCAS 1991, 1991, Singapore, pp. 2518-2521.
[68] R. Lee, \Measures of Correlation for Waveform Data in Clinical Research",
http : ==biomech:brighton:ac:uk=help=cmc=.
[69] R. E. McAllister and D. Noble, \The Time and Voltage Dependence of the
Slow Outward Current in Cardiac Purkinje Fibres", Journal of Physiology,
Vol. 186, No. 3, October 1966, pp. 632-662.
[70] G. W. Beeler and H. Reuter, \Membrane Calcium Current in Ventricular
Myocardium Fibres", Journal of Physiology, Vol.207, No. 1, March 1970, pp.
191-209.
225
[71] H. Reuter, \Localization of Beta Adrenergic Receptors and Efiects of No-
radrenaline and Cyclic Nucleotides on Action Potentials, Ionic Currents and
Tension in Mammalian Cardiac Muscle.", Journal of Physiology, Vol. 242,
No. 2, October 1974, pp. 429-451.
[72] L. S. Gettes and H. Reuter, \Slow Recovery from Inactivation of Inward
Currents in Mammalian Myocardial Fibres", Journal of Physiology, Vol. 240,
No. 3, Augest 1974, pp. 703-724.
[73] W. Trautwein, T. F. Mcdonald, and O. Tripathi, \Calcium Conductance and
Tension in Mammalian Ventricular Muscle.", P gers Archiv European Journal
of Physiology, Vol. 354, No. 1, March 1975, pp. 55-74.
[74] G. A. Rincon-Mora, \Active Capacitor Multiplier in Miller-Compensated
Circuits", IEEE Transactions on Solid-State Circuits, Vol. 35, No. 1, January
2000, pp. 26-32.
[75] Y. Tang, M. Ismil, and S. Bibyk, \Active Miller Capacitor Multiplier for
Compact On-chip PLL Filter", IEEE Electronics Letters, Vol. 39, No. 1,
January 2003, pp. 43-45.
[76] C. J. Chang and K. H. Chen, \Bidirectional Current-Mode Capacitor Mul-
tiplier in DC-DC Converter Compensation", Proceedings of the Fifth Inter-
national Workshop on System-on-Chip for Real-Time Applications, Alberta
Canada, July 2005, pp. 111-116.
[77] Y. Tsividis, Operation and Modeling of the MOS Transistor, 2nd edition,
McGraw-Hill, 1999.
[78] Q. Zhang, J. J. Liou, J. R. McMacken, J. Thomson, and P. Layman, \SPICE
Modeling and Quick Estimation of MOSFET Mismatch Based on BSIM3
Model and Parametric Tests", IEEE Journal of Solid-State Circuits, Vol.
36, No. 10, October 2001, pp. 1592-1595.
[79] O. S. Unsal, J. W. Tschanz, K. Bowman, V. De, X. Vera, A. Gonzalez, and O.
Ergin, \Impact of Parameter Variations on Circuits and Microarchitecture",
IEEE Micro, Vol. 26, No. 6, November/December 2006, pp. 30-39.
[80] H. Onodera, \Variability: Modeling and Its Impact on Design ", IEIE Trans-
actions on Electronics, Vol. E89-C, No. 3, March 2006, pp. 342-348.
[81] \Wafer Electrical Test Data and SPICE Model Parameters AMIS ABN
(1.50 micron)", http : ==www:mosis:com=Technical=Testdata=ami?abn?
prm:html.
226
[82] L. Clerc, \Directional Difierence of Impulse Spread in Trabecular Muscle from
Mammalian Heart", Journal of Physiology, vol. 255, No. 2, February 1976,
pp. 335-346.
[83] A. V. Holden and A. V. Panfllov, \Modelling Propagation in Excitable Me-
dia", Computational Biology of the Heart, John Wiley & Sons, chapter 3,
March 1997, pp. 65-100.
[84] A. N. Kolmogorov, I. G. Petrovskii, and N. S. Piskunov, \A Study of the
Difiusion Equation with Increase in the Amount of Substance, and Its Appli-
cation to a Biological Problem", Selected Works of A. N. Kolmogorov, Kluwer
Academic Publishers, 1991, pp. 242-270.
[85] C. V. Pao, \Nonlinear Parabolic and Elliptic Equations", Plenum Press, New
York, 1992, pp. 214-219.
[86] A. M. Turing, \The Chemical Basis of Morphogenesis", The Royal Society of
London Philosophical Transactions. Series B, Vol. 237, No. 641, August 1952,
pp. 37-72.
[87] Y. B. Pesin and A. A. Yurchenko, \Some Physical Models of the Reaction-
Difiusion Equation, and Coupled Map Lattices", Russian Mathematical Sur-
veys, Vol. 59, No. 3, 2004, pp. 481-513.
[88] I. R. Epstein and J. A. Pojman, \An Introduction to Nonlinear Chemical
Dynamics", Oxford University Press, 1998, pp. 347-348.
[89] R. C. Barr, R. Plonsey, and C. R. Johnson, \Membrane Current from Trans-
membrane Potentials in Complex Core-conductor Models", IEEE Transac-
tions on Biomedical Engineering, Vol. 50, No. 4, April 2003, pp. 405-411.
[90] M. S. Spach, \The Discontinuous Nature of Electrical Propagation in Cardiac
Muscle", Journal Annals of Biomedical Engineering, Vol. 11, No. 3, May 1984,
pp. 208-261.
[91] W. Thomson (Lord Kelvin), \On the Theory of the electric Telegraph", The
Proceeding of the Royal Society of London, Vol. 7, 1854/1855, pp. 382-399.
[92] J. Clark and R. Plonsey, \A Mathematical Evaluation of the Core Conductor
Model", Biophysical Journal, Vol. 6, No. 1, January 1966, pp. 95-112.
[93] B. Chen, \Mathematical Models of Motion Detection in the Fly?s Vi-
sual Cortex", http : ==etd:lib:ttu:edu=theses=available=etd ? 11152005 ?
114857=unrestricted=Chen Baili diss:pdf, MathematicalDepartment, Texas
Tech University, December 2005.
227
[94] I. R. Eflmov, \Cardiac Bioelectricity: Imaging and Electrotherapy", http :
==efimov:wustl:edu=class=BME140=handout:pdf, Department of Biomedi-
cal Engineering, Washington University, September 2005.
[95] A. G. Kleber and Y. Rudy, \Basic Mechanisms of Cardiac Impulse Propa-
gation and Associated Arrhythmias", Physiological Reviews, Vol. 84, No. 2,
April 2004, pp. 431-488.
[96] N. F. Hooke, \E?cient Simulation of Action Potential Propagation in a Bido-
main", Department of Computer Science, Duke University, September 1992.
[97] N. Sperelakis and L. Ramasamy, \Modeling Electric Field Transfer of Excita-
tion at Cell Junctions", IEEE Engineering in Medicine and Biology, Vol. 21,
No. 6, November-December 2002, pp. 130-143.
[98] E. Entcheva and F. J. Claydon, \A Discrete Model of the Dynamic Behavior
of the Cardiac Muscle", Computers in Cardiology (1996 IEEE Conference on
Computers in Cardiology), Indianapolis USA, September 1996, pp. 609-612s.
[99] A. P. Henriquez, R. Vogel, B. J. Muller-Borer, C. S. Henriquez, R. Weingart,
and W. E. Cascio, \In uence of Dynamic Gap Junction Resistance on Impulse
Propagation in Ventricular Myocardium: A Computer Simulation Study",
Biophysical Journal, Vol. 81, No. 4, October 2001, pp. 2112-2121.
[100] J. Wu and D. P. Zipes, \Efiects Of Spatial Segmentation In The Continuous
Model Of Excitation Propagation In Cardiac Muscle", Journal of Cardiovas-
cular Electrophysiology, vol. 10, No. 7, 1999, pp. 965-972.
[101] A. M. Pertsov, \Three-Dimensional Organization of Reentry in Fibrillating
Ventricular Wall", Computer Simulation and Experimental Assessment of
Cardiac Electrophysiology, Blackwell Publishing, 1st edition, New York, July
2001, pp. 63-68.
[102] M. C. Trudel, R. M. Gulrajani, and L. J. Leon, \Simulation of Propagation
in a Realistic-Geometry Computer Heart Model with Parallel Processing",
IEEE 2001 Proceedings of the 23rd Annual EMBS International Conference,
October 25-28, Istanbul, Turkey, October 2001, pp. 359-362.
[103] P. C. Franzone, P. Deu hard, B. Erdmann, J. Lang, and L. F. Pavarino,
\Adaptivity in Space and Time for Reaction-Difiusion Systems in Electro-
cardiology", SIAM Journal on Scientiflc Computing , Vol. 28, No. 3, March
2006, pp. 942-962.
[104] O. H. Schmitt, \An electrical theory of nerve impulse propagation", The
American Journal of Physiology, Vol. 119, No. 2, June 1937, pp. 399-400.
228
[105] J. Nagumo, S. Yoshizawa, and S. Arimoto, \Bistable Transmission Lines",
IEEE Transactions on Circuit Theory, Vol. CT-12, No. 3, September 1965,
pp. 400-412.
[106] M. R. Roussel, \Reaction Difiusion Equations", http : ==people:uleth:ca= ?
roussel=nld=Turing:pdf.
[107] S. P. Ellner and J. Guckenheimer, \Spatial Patterns in Biology", Dynamic
Models in Biology, Princeton University Press, chapter 7, March 2006, pp.
217-242.
[108] M. A. Bray and J. P. Wikswo, \Examination of Optical Depth Efiects on
Fluorescence Imaging of Cardiac Propagation", Biophysical Journal, Vol. 85,
No. 6, December 2003, pp. 4134-4145.
[109] F. H. Fenton, A. Karma, H. M. Hastings, and S. J. Evans, \Transition
from Ventricular Tachycardia to Ventricular Fibrillation as Function of Tissue
Characteristics in a Computer Model", Proceedings of the 22nd Annual EMBS
International Conference, July 23-28, 2000, Chicago U.S.A., pp. 391-394.
[110] R. Groningen, \Molecular Adaptations in Human Atrial Fibrillation: Mech-
anisms of Protein Remodeling", http : ==disertations:ub:rug:nl=FILES=
faculties=medicine=2000=b:j:j:m:brundel =titlecon:pdf, 2000.
[111] S. Sinha and D. J. Christini \Termination of Reentry in an Inhomogeneous
Ring of Model Cardiac Cells", Physical Review E, Vol. 66, No. 6, December
2002, pp. 061903.1-061903.7.
[112] M. Courtemanche, \Complex Spiral Wave Dynamics in a Spatially Dis-
tributed Ionic Model of Cardiac Electrical Activity", Chaos, Vol. 6, No. 4,
July 1996, pp. 579-600.
[113] A. Xu and M. R. Guevara, \Two Forms of Spiral-Wave Reentry in an Ionic
Model of Ischemic Ventricular Myocardium", Chaos, Vol. 8, No. 1, November
1997, pp. 157-174.
[114] C. R. Johnson and R. C. Barr, \Termination of Reentrant Propagation by
a Single Extracellular Stimulus in an Atrial Ring Model", Proceedings of the
23rd Annual EMBS International Conference, October 25-28, 2001, Istanbul
Turkey, pp. 398-401.
[115] F. Xie, Z. Qu, A. Garflnkel, and J. N. Weiss, \Electrical Refractory Period
Restitution and Spiral Wave Reentry in Simulated Cardiac Tissue", American
Journal of Physiology - Heart and Circulatory Physiology, Vol. 283, No. 5,
November 2002, pp. 448-460.
229
[116] M. Philip, \Mathematical Modelling of Cardiac Arrhythmias", http :
==www:f:kth:se= ? f97?mph=rapport=rapport:html.
[117] P. R. Gray and R. G. Meyer, \Analysis and Design of Analog Integrated
Circuits", John Wiley & Sons, Inc., 3rd edition, 1993.
[118] P. Zavracky, M. Zavracky, D-P Vu, and B. Dingle, \Three-Dimensional
Processor Using Transferred Thin Film Circuits", US patent 5,656,548, Patent
and Trademark O?ce, Washington, D.C., August 1997.
[119] D. Bartholomew, \Winning Technologies: Flexible Transistors", http :
==www:industryweek:com=ReadArticle:aspx?ArticleID = 771&Section
ID = 38, December 2000.
[120] J. Ferrero (Jr), V. Torres, F. Montilla, and E. Cplpmar, \Simulation of
Reentry during Acute Myocardial Ischemia: Role of ATP-Sensitive Potassium
Current and Acidosis", Computers in Cardiology (2000 IEEE Conference on
Computers in Cardiology), Cambridge MA USA, September 2000, pp. 239-
242.
[121] N. V. Thakor, J. Ferrero (Jr), J. Salz, B. I. Gramatikov, and J. Ferrero (Sr),
\Electrophysiologic Models of Heart Cells and Cell Networks", Engineering
in Medicine and Biology Magazine, IEEE, Vol. 17, No. 5, September/October
1998, pp. 73-83.
[122] R. M. Shaw and Y. Rudy, \Electrophysiologic Efiects of Acute Myocardial
Ischemia : a Theoretical Study of Altered Cell Excitability and Action Poten-
tial Duration", Cardiovascular Research, Vol. 35, No. 2, 1997, pp. 256-272.
[123] S. Seven, S. Vecchietti, I. Rivolta, C. Napolitano, S. G. Priori, and S. Cav-
alcanti, \Action Potential Changes Due to Y1795H Mutation in Brugada
Syndrome Patients: A Simulation Study", Computers in Cardiology (2003
IEEE Conference on Computers in Cardiology), Thessaloniki Greece, Septem-
ber 2003, pp. 437-440.
[124] R. Vogel and R. Weingart, \Mathematical Model of Vertebrate Gap Junc-
tions Derived from Electrical Measurements on Homotypic and Heterotypic
Channels", Journal of Physiology, Vol. 510, No. 1, 1998, pp. 177-189.
[125] F. J. Jaeger, \Cardiac Arrhythmias", http : ==www:clevelandclinicmeded
:com=DISEASEMANAGEMENT=cardiology=arrhythmias=arrhythmias
:htm.
[126] J. Saiz, J. M. Ferrero, L. M. Roa, and N. V. Thakor, \Simulation of triggered
Activity and Abnormal Automaticity Inventricular Myocytes", Computers
230
in Cardiology (1995 IEEE Conference on Computers in Cardiology), Vienna
Austria, September 1995, pp. 345-348.
[127] C. H. Luo and Y. Rudy, \A Dynamic Model of the Cardiac Ventricular Ac-
tion Potential. II. Afterdepolarizations, Triggered Activity, and Potentiation",
Circulation Research, Vol. 74, No. 6, June 1994, pp. 1097-1113.
[128] J. Zeng and Y. Rudy, \Early Afterdepolarizations in Cardiac Myocytes:
Mechanism and Rate Dependence", Biophysical journal, Vol. 68, No. 3, March
1995, pp. 949-964.
[129] \Anti-arrhythmic Drugs", http : ==www:anaesthetist:com=icu=organs=
heart=rhythm=Findex:htm#arhythm:htm.
[130] C. F. Starmer, \How Antiarrhythmic Drugs Increase the Rate of Sudden
Cardiac Death", International Journal of Bifurcation and Chaos, vol. 12, No.
9, 2002, pp. 1953-1968.
[131] C. F. Starmer, A. O. Grant, and T. J. Colatsky, \What Happens When
Cardiac Na Channel Function is Compromised? 2: Numerical Studies of the
Vulnerable Period in Tissue Altered by Drugs", Cardiovascular Research, vol.
57, No. 4, March 2003, pp. 1062-1071.
[132] O. Bernus, B. V. Eyck, H. Verschelde, and A. V. Panfllov, \Transition From
Ventricular Fibrillation to Ventricular Tachycardia: a Simulation Study on
the Role of Ca2+-channel Blockers in Human Ventricular Tissue", Physics in
Medicine and Biology, Vol. 47, No. 23, December 2002, pp. 4167-4173.
[133] A. Cimponeriu, C. F. Starmer, and A. Bezerianos, \A Theoretical Analysis
of Acute Ischemia and Infarction Using ECG Reconstruction on a 2-D Model
of Myocardium", IEEE Transactions on Biomedical Engineering, Vol. 48, No.
1, 2001, Janurary 2001, pp. 41-54.
[134] R. H. Clayton and A. V. Holden, \Computational Framework for Simulating
the Mechanisms and ECG of Re-entrant Ventricular Fibrillation", Physiolog-
ical Measurement, Vol. 23, No. 4, November 2002, pp. 707-726.
[135] H. W. Chiu and T. Kao, \A Mathematical Model for Autonomic Control of
Heart Rate Variation", IEEE Engineering in Medicine and Biology Magazine,
Vol. 20, No. 2, March/April 2001, pp. 69-76.
[136] H. R. Warner and A. Cox, \A Mathematical Model of Heart Rate Control by
Sympathetic and Vagus Efierent Information", Journal of Applied Physiology,
Vol. 17, No. 2, March 1962, pp. 349-355.
231
[137] A. C. Sanderson, \Input-output Analysis of an IPFM Neural Model: Efiects
of Spike Regularity and Record Length", IEEE Transactions on Biomedical
Engineering, Vol. 27, No. 3, March 1980, pp. 120-131.
[138] P. J. Hunter, A. D. McCulloch, and H. E. D. J. ter Keurs, \Modelling
the Mechanical Properties of Cardiac Muscle", Progress in Biophysics and
Molecular Biology, Vol. 69, No. 2, March 1998, pp. 289-331.
[139] J. J. Rice, R. L. Winslow, and W. C. Hunter, \Comparison of Putative
Cooperative Mechanisms in Cardiac Muscle: Length Dependence and Dy-
namic Responses", American Journal of Physiology - Heart and Circulatory
Physiology, Vol. 276, No. 5, May 1999, pp. H1734-H1754.
232