Abstract
Background
Brain, heart and skeletal muscle share similar properties of excitable tissue, featuring both discrete behavior (allornothing response to electrical activation) and continuous behavior (recovery to rest follows a temporal path, determined by multiple competing ion flows). Classical mathematical models of excitable cells involve complex systems of nonlinear differential equations. Such models not only impair formal analysis but also impose high computational demands on simulations, especially in largescale 2D and 3D cell networks. In this paper, we show that by choosing Hybrid Automata as the modeling formalism, it is possible to construct a more abstract model of excitable cells that preserves the properties of interest while reducing the computational effort, thereby admitting the possibility of formal analysis and efficient simulation.
Results
We have developed CellExcite, a sophisticated simulation environment for excitablecell networks. CellExcite allows the user to sketch a tissue of excitable cells, plan the stimuli to be applied during simulation, and customize the diffusion model. CellExcite adopts Hybrid Automata (HA) as the computational model in order to efficiently capture both discrete and continuous excitablecell behavior.
Conclusions
The CellExcite simulation framework for multicellular HA arrays exhibits significantly improved computational efficiency in largescale simulations, thus opening the possibility for formal analysis based on HA theory. A demo of CellExcite is available at http://www.cs.sunysb.edu/~eha/ webcite.
Background
An excitable cell has the ability to propagate an electrical signal—known at the cellular level as the Action Potential (AP)—to surrounding cells. An AP corresponds to a change of potential across the cell membrane, and is caused by the flow of ions between the inside and outside of the cell. The major ion species involved in this process are sodium, potassium and calcium; they flow through multiple voltagegated ion channels as poreforming proteins in the cell membrane. Excitation disturbances can occur in the behavior of these ion channels at the cell level, or in the propagation of the electrical waves at the cellnetwork level. Examples of excitable cells are neurons, cardiac myocytes and skeletal muscle cells. Excitable cell networks are important in the normal functioning and in the pathophysiology of many biological processes. In particular, the rhythmic, pumplike function of the heart is driven by muscle contractions, which are in turn triggered by cellgenerated electrical signals (excitations). Of special interest are cardiac arrhythmias: disruptions of the normal excitation process due to faulty processes at the cellular level, single ionchannel level, or at the level of celltocell communication. The clinical manifestation is a rhythm with altered frequency (tachycardia or bradycardia) or the appearance of multiple frequencies (polymorphic Ventricular Tachycardia) with subsequent deterioration to a chaotic signal (Ventricular Fibrillation). VF [1] is a typically fatal condition in which there is uncoordinated contraction of the cardiac muscle of the ventricles in the heart. As a result, the heart fails to adequately pump the blood, and hypoxia may occur.
Excitable tissue is modeled in terms of reactiondiffusion systems. Thus, a typical continuous representation would involve partial differential equations (PDEs) for the diffusing species (typically the transmembrane potential), and a system of nonlinear ordinary differential equations describing all other state variable that are normally considered nondiffusing. These may include ionchannel gating variables and ion concentrations. The first mathematical model of ionic processes that underly cell excitation was empirically developed in 1952 by Hodgkin and Huxley (HH) for a squid giant axon [2]. This provided the basis for subsequent models of increasing complexity, using multiple continuous state variables (voltage, ionchannel gates, ion concentrations) to describe APs in different cell types [35]. Current models of cardiac cells include more than 20 such state variables and a large number of fitted parameters. Detailed models of cardiac excitation are perceived as overdetermined systems and, as such, make both qualitative—i.e. checking general properties—and quantitative analysis—i.e. by simulation—at the organ or even tissue level impractical. At the opposite end of the spectrum, completely discrete models based on cellular automata (CA) have emerged [6,7].
The first generation of CA models used nearestneighbor diffusion modeling (Neumann and Moore neighborhoods) and a small number of discrete states, resulting in unrealistic AP shape and wave propagation. Secondgeneration CA models [7] focused on correct representation of wavefront curvature effects by employing more complex neighborhood functions, such as Gaussian, circular templates or randomized lattices. Furthermore, the transitions rules for the relaxation states were updated to reflect a higher threshold for excitation and to effectively represent the relative and absolute refractory period. The latest generation is exemplified by Barkley's model [8], in which a standard finitedifference method is used to calculate the diffusive term, but CAlike rules govern the kinetics of the two model variables, with adjustable thresholds. Recently, modified CA models have been used to study cardiac excitability and for comparison with experimental data [6,9]. A body of literature provides clear links between the classical continuous PDE representation and the more ad hoc CAbased approach as an alternative description of reactiondiffusion systems. The purely discrete nature of CA presents some difficulties in capturing subtle nonstepwise features of excitation.
One way to reduce the complexity of models based on differential equations while preserving the fundamental features of these systems is to construct a more abstract model that preserves the properties of interest. One promising approach is based on the use of Hybrid Automata (HA) [10,11] as a modeling formalism for complex biological processes. More formally, HA are an extension of finite automata that allows one to associate a continuous behavior with each state. The approach of [12] demonstrated the feasibility of using HA as a modeling formalism for excitable cells. The biological behavior of such cells is intrinsically hybrid in nature, featuring both discrete (allornothing response to electrical behavior) and continuous behavior (recovery to rest follows a temporal path, determined by multiple competing ion flows). Starting from a biological interpretation of their APs, 4state HA models have been derived in [12] for several classical excitablecell types. In this paper, we present CellExcite, a sophisticated simulation environment for excitablecell networks based on these 4state models.
Results and discussion
CellExcite allows the user to sketch a tissue of excitable cells, plan the stimuli to be applied during simulation, and customize the diffusion model. As Figure 1 illustrates, the architecture of CellExcite consists of two main components:
Figure 1. General Architecture of CellExcite.
CellExcite Graphical User Interface (GUI) This component provides several panels that the user can access in order to customize the features of an excitablecell network. Using the Tissue Panel, the user can specify the tissue size, simulation time, and list of stimuli to be applied during simulation. The Parameters Panel allows the user to select the desired 4state hybrid automaton representing the behavior of a particular type of excitable cell: HH for neuron, NNR for cardiac myocytes. Furthermore a user can specify the distance between two neighboring cells, the membrane capacitance, and the time step of the simulation. The Diffusion Panel enables the user to select from among different lattices that could better approximate the disposition of cells in a 2D tissue. Through this panel, the user can also choose the radius of the voltage influence a cell has on its neighbors.
EventDriven HAbased Simulator This component has been implemented by extending the eventdriven approach described in [13] with the following new features:
• The event priority queue is optimized: duplicate events, i.e. multiple events of the same type for the same cell and time step generated by neighboring cells, are eliminated.
• Several neighborhood functions are added: cell networks can be represented both as a triangular and square lattice. The diffusion and electrical propagation within a cell network are modeled with an exponential neighborhood function.
• Colored snapshots and video representing both the discrete and continuous behavior of the system can be generated.
An example insilico simulation of excitable cells
In this section, we provide an example simulation of an excitablecell network using CellExcite. The goal is to first simulate ventricular fibrillation and then defibrillation on a tissue of NNR (neonatal rat) cardiac myocytes arranged in a 400x400 array. We wish to simulate this network for the first 500 ms with a time step of 1^{−3} ms.
Sketching a tissue and planning the stimuli to be applied using the Tissue Panel
To carry out the simulation, we must first perform the following operations:
• Resize the tissue to 400x400 cells.
• Set the simulation time to 500 ms.
• Apply at the beginning of the simulation, for 1 ms, a first stimulus of 800 μA/cm^{2} to a rectangular area of the tissue, from row 310 on the top of the tissue to row 395 on the bottom, and from column 5 on the left to column 15 on the right.
• Apply 145 ms after the beginning of the simulation, for 1 ms, a second stimulus of 1000 μA/cm^{2} to a rectangular area of the tissue, from row 235 on the top of the tissue to row 245 on the bottom, and from column 0 on the left to column 150 on the right.
• Apply 400 ms after the beginning of the simulation, for 1 ms, a third stimulus of 800 μA/cm^{2} to a rectangular area of the tissue, from row 4 on the top of the tissue to row 394 on the bottom, and from column 4 on the left to column 394 on the right.
As Figure 2 shows, the Tissue Panel of the CellExcite GUI allows the user to insert all these data using a convenient visual framework. Figure 2 contains a snapshot of the GUI while the user is planning the second stimulus on the tissue.
Figure 2. Planning a stimulus using Tissue Panel.
Setting the cell parameters
The next step is to set the singlecell parameters, such as the choice of 4state HA model, the distance between two neighboring cells, the cell's membrane capacitance, and the time step of the simulation. For our example, we proceed as follows:
• Choose NNR as the 4state HA model representing the behavior of a single cell
• Set the distance between two neighboring cells to 0.01 cm
• Set the membrane capacitance to 1 μF/cm
• Set the simulation time step to 0.001 msec
Figure 3 a) is a snapshot of the GUI while the user is setting the singlecell parameters.
Figure 3. Setting parameters and diffusion panel.
Selecting a diffusion model
To complete the experiment, we need to specify the lattice on which we would like to dispose the cell network. The DiffusionPanel of CellExcite provides two different lattices: triangular and square. Furthermore, a user can specify a cell's radius of influence with respect to neighbor cells. A gradient of colors, from yellow to red, indicate respectively the minor or major weight, based on a normalized exponential function, of neighbor cells with respect to the cell disposed in the middle of the panel. The higher the weight, the larger the number of neighbor cells taken into account. In this case, the simulation of the diffusion process is more precise but also more computationally intensive, as the number of computations is increased. As Figure 3 b) shows, for our experiment, we choose a triangular lattice with a radius of 4.
Results of simulation
Figure 4 shows the results obtained with CellExcite, by simulating our insilico experiment. The top of the figure depicts the continuous behavior of the tissue during the simulation. This is a quantitative analysis showing the voltage of each cell of the tissue at a particular time step. The bottom of the figure depicts the discrete behavior of the tissue. This is a qualitative analysis showing, for each cell of the tissue and for each time step, the AP state of the hybrid automaton: resting, stimulated, upstroke and plateau. As the sequence of pictures shows, the first stimulus has a normal propagation, while the second generates two counterclockwise spirals because there is not enough space for normal propagation. After 150 ms, the heads of the two spirals collide and generate other spirals. This degeneration of the electrical propagation is the cause of fibrillation. A third stimulus at is applied at 400 ms to all cells in the tissue in order to recover them the resting state. This phenomena is comparable to defibrillation.
Figure 4. Results of the simulation.
Performance comparison
To illustrate the performance gains achieved in an HAbased simulation framework for excitable cells versus a traditional ODEbased one, consider Table 1. There we give the CPU times required for the two frameworks on a 2second simulation of NNR tissue, with the tissue size ranging from a 2x2 cell array to a 400x400 cell array. On the more substantial cases (i.e., array sizes of 16x16 on up), it can be seen that the HAbased simulation framework enjoys an almost 8fold speedup over the ODEbased one. All results, which were originally reported in [12], were produced on a PC equipped with a Pentium Intel 4 CPU 3.00GHz with 1GB of memory.
Table 1. Performance comparison for 2second simulation
Availability and requirements
• Project Name: Excitable Hybrid Automata (EHA)
• Project HomePage:[14]
• Operating System: Windows, Linux
• Programming Language: C, Java
• Licence: The CellExcite software package is available under the GNU Less General Public License (LGPL) license. Please contact the first author for details.
Conclusions and future work
We developed CellExcite, a hybridautomatabased visualization framework for excitablecell networks. CellExcite provides a user interface that allows the user to sketch a tissue of excitable cells, plan the stimuli to be applied during simulation, and customize the diffusion model. The use of multicellular HA to model networks of excitable cells is a reasonable compromise to reduce the computational demands of large scale 2D cellnetwork simulations, without losing the ability to capture important features of the action potential such as restitution and refractoriness.
Besides quantitative analysis obtained by simulation, this formal model opens up the possibility of qualitative analysis. To this end, as future work, we aim to extend CellExcite with tools that allow the runtime verification of specific spatiotemporal patterns such as the creation of dangerous spirals. This could be helpful to find automatically an energyefficient strategy to counteract degeneration of the normal electrical propagation. This modelbased control could be exploited for the construction of the nextgeneration cardio defibrillator. Furthermore, we would like to exploit agentbased technology to better distribute the computation on a gridcomputing environment.
Methods
Action potential
The electrical signal at the cellular level is known as an action potential (AP). Action potentials for ventricular cells (the major cells in the heart) are externally triggered events: a cell fires an action potential as an allornothing response to a suprathreshold electrical signal, and each AP follows more or less the same sequence of events and has the same magnitude regardless of applied stimulus. An AP lasts for a couple of hundred milliseconds in most mammals. During the AP, no reexcitation can occur, which is a safety mechanism to ensure the reliable working of the heart. The early portion of an AP is known as the “absolute refractory period” due to its nonresponsiveness to further simulation. The later portion of an AP is known as the “relative refractory period”, during which an altered secondary excitation event is possible if the stimulation threshold is raised.
Despite differences in AP duration, morphology and underlying ion currents between different species and different regions in the heart, the following major AP phases can be identified: resting phase, rapid upstroke, early repolarization phase, plateau or later repolarization phase, and final repolarization (identical to the resting state due to the cyclic nature of an AP). The resting state features a constant transmembrane potential (difference between the inside and outside potential of the cell) of about −80 mV for most species; i.e. the membrane is polarized at rest. During the AP upstroke, the transmembrane potential rapidly changes (over the course of a couple of milliseconds) from negative to positive; i.e. the membrane depolarizes. This is followed by an early repolarization phase. A slower, plateau phase is present in most mammalian action potentials, during which calcium influx facilitates the muscle contraction. After this phase, a faster initial repolarization brings the potential back to the resting state. Because of the universal nature of these AP features between species and regions, we use them as a guide in the construction of HA models [12].
The HodgkinHuxley Model
The first quantitative description of cellular excitation was empirically developed by Hodgkin and Huxley (HH) for a squid giant axon [2]. The HH model includes three ionic currents: fast inward sodium, outward potassium, and a timeindependent linear (leak) current. The generalized form of HH model is as follows:
where V is the transmembrane voltage [mV], whose variation forms the AP; are the maximum channel conductance [mS/μ F] for the sodium (Na), potassium (K) and the leakage channel (L), respectively; E_{Na}, E_{K}, E_{L} are the reversal potentials [mV] for the sodium, potassium and the leakage channel, respectively; m, h, n are voltage and timedependent ion channel gates, following the same general differential equation in y, where y_{inf} and τ_{inf} represent the steadystate and the time constant of a gate; C is the cell capacitance [μF], and I_{st} is the stimulation current [μA/μF].
LuoRudy Guinea Pig Ventricular Cell Model
In a series of papers, Y Rudy et al. have developed some of the most detailed cardiac cell models to date, targeting guinea pig [5,15]. The ionchannel description in these models follows the same framework as the HH model, but a much larger number of ion currents is included. The complexity of this class of models is further increased by the addition of active ion pumps, intracellular compartments for calcium transport and calcium buffers. The detailed description of the LRd model is omitted here.
Neonatal Rat Ventricular Cell Model
Among the mammalian species, the mouse and the rat have a substantially different AP morphology—much more triangular with almost absent plateau phase—compared to the AP simulated by the LRd model. Neonatal rats are often used as an experimental model in cardiac electrophysiology, and a computational model is a derable tool. A neonatal rat model (NNR) is being developed by Entcheva et al., derived from the LRd model. In [12], Pye et al. use a hybridautomaton formulation (following the same structure as for LRd) with adjusted parameters to replicate the behavior of this detailed ionic model.
Modeling Action Potential using Hybrid Automata
To permit formal analysis and to increase simulation efficiency, it could be useful to perform abstraction on a set of nonlinear differential equations describing the behavior an excitable cell to obtain a hybrid automaton. Formally, a hybrid automaton is defined as follows.
A Hybrid automaton H consists of the following [11]:
 As finite set X = x_{1}, ⋯, x_{n} of realnumbered variables. The number n is called the dimension of H. We write of dotted variables (which represent first derivatives of variables in X), and X′ = x′_{1},⋯, x′_{n} of primed variables (which represent values of variables X at the conclusion of discrete steps).
 A finite, discrete control graph (V,E). The vertices in V are called control modes. The edges in E are called control switches.
 Three vertex labeling functions init, inv and flow that assign to each control mode e υ ∈ V three predicates. Each initial condition init(v) and invariant condition inv(v) is a predicate whose free variables are from X. Each flow condition flow(v) is a predicate whose free variables are from .
 An edge labeling function jump that assigns to each control switch e ∈ E a predicate. Each jump condition jump(e) is a predicate whose free variables are from X ∪ X′
 A finite set Σ of events, and an edge labeling function event that assigns to each control switch an event.
The HA models chosen have four control modes: resting and final repolarization (FR), stimulated, upstroke and, plateau and early repolarization (ER). Initially, the cell is in the resting and final repolarization mode. When (externally) stimulated with the event V_{S}, the cell enters the stimulated mode and updates its voltage according to the stimulus current. Upon termination of the stimulation, via event , with a subthreshold voltage, the cell returns back to resting mode without firing AP. If the stimulus is suprathreshold, i.e., V_{s} >V_{T} holds, the excited cell will generate an action potential by progressing to the upstroke mode. The recovery course of the cell follows the transitions to mode plateau and early repolarization and then to resting and final repolarization. The guards on the control switches monitor the transmembrane potential, rather than imposing a rigid timing scheme. This approach allows for AP adaptation (response to various pacing frequencies).
HA for the HH model
The HA for the HH model is shown in Figure 5. Variables υ_{x} and υ_{y} define a secondorder system of linear differential equations in each control mode. They are an abstraction of the ionic currents and the gates. I_{st} is the excitation current and V_{S} is the stimulation event. The membrane voltage is V = υ_{x} − υ_{y}. The initial control mode is q_{0}. The mode invariants are given below the differential equations describing the transmembrane voltage. Like the switch guards, they depend on three modelspecific (see Table 2): threshold voltage V_{T}, overshoot voltage V_{O}, and repolarization voltage V_{R}. The transition guards are bracketed.
Figure 5. Hybrid automata for HH model.
HA for the LRd model
The modeling framework for the LRd model is similar to that for HH. However, to properly represent the longer maintained plateau phase of the cardiac AP and to capture its frequency adaptation, we extend the hybrid model with additional variables. A third continuous variable, υ_{z}, is added. The need for such a variable in the LRd and NNR models can be explained by the major difference in the ion fluxes between neurons and cardiac cells; namely, calcium flux plays a profound role in the maintenance of the AP plateau for proper cardiac muscle contraction to take place.
Additionally, a new restitutionrelated continuous variable υ′ is added to the LRd model, which is used to modify the overall voltage by reflecting changes in the diastolic interval (DI). DI is the time period between AP's recovery and a new stimulation, i.e., the rest period between successive stimulations. It is known that the immediate memory of an excitable cell is directly linked to the DI: a shorter DI results in a shorter following AP, while a longer DI produces a longer AP [16]. This simple memory model helps capture the proper response of AP to pacing frequency, which is an essential feature of the cardiac excitation. The resulting system of differential equations in the corresponding modes is however no longer linear. The hybrid automaton used to simulate LRd model is defined in Figure 6 with .
Figure 6. Hybrid automata for LRd model.
HA for the NNR model
The HA for NNR is defined in Figure 7. Here, f(θ) = 1 + 2θ. For improved modeling of cellcell interactions in cardiac excitation, the threshold voltages do not remain constant during simulation; instead, they also become a function of θ = V′/V_{R}, defined as follows: and
Figure 7. Hybrid automata for NNR model.
Parameters definitions
The values of the coefficients and constants occurring in the HH, LRd and NNR HA models are summarized in Table 2. They were obtained either from the cited literature or empirically through experimentation.
Table 2. Parameters for HH LRd and NNR 4state HA models
Diffusion Modeling in Multicellular HA
The spatially extended model of multiple connected excitable cells is defined as below, where the lefthand side represents the diffusion term (propagation of the transmembrane voltage), while the righthand side describes the activity at each cell, according to the selected model (detailed, PDEbased (1) or hybrid (2)):
(1)
(2)
where ∇ is the Laplacian Operator in 2D; σ is the 2x2 conductivity tensor; β [cm^{−1}] is the surface area to volume ratio for the cells; I_{ion}[mA/cm^{2}] is the sum of all ion currents and I_{st} is the stimulus current. In a preliminary study we coupled the HA description of each cell to classical Laplacianbased diffusion modeling equations, solved by a finite element numerical scheme [13]. Due to the matrix operations involved, this solution is not the most computationally efficient. In this study we extended the spatial modeling of coupled hybrid automata cells by employing alternative diffusion modeling operators. For correct modeling of the diffusion of the electrical potential, simulation data have to be consistent with the above equation for the medium of interest. The Laplacianbased solution of the diffusion term has been substituted with optimized neighborhood functions to match the experimentally obtained dispersion curve, since the latter requires extended neighborhoods, thus slowing computations. The exponential neighborhood function that we used is defined as follows:
where r_{max} is maximum radius of influence of a cell respect to its neighbor cells and d <r_{max} is the distance of a cell from a particular neighbor cell. This function provides the weight that assume the voltage of each cell respect to the others. The closer a cell is to its neighbors, the greater the influence of its voltage on its neighbor cells' voltage in the next time step. This function is dependent on the chosen lattice topology. As Figures 8 and 9 illustrate, in CellExcite both triangular and square lattices are supported. Table 3 and 4 report the weights of the cells for each lattice as a function of the distance.
List of abbreviations used
• AP: Action Potential
• CA: Cellular Automata
• DI: Diastolic Interval
• GUI: Graphical User Interface
• HA: Hybrid Automata
• HH: Hodgkin Huxley
• LGPL: Less General Public License
• NNR: NeoNatal Rat
• PDE: Partial Differential Equation
• VF: Ventricular Fibrillation
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
EB wrote the paper and participated in the tool design and development. FC, EE, RG and SAS conceived the study, coordinated the tool design, and helped to draft the manuscript. All authors read and approved the final manuscript.
Acknowledgements
Research supported in part by the Italian FIRBMIUR LITBIO: Laboratory for Interdisciplinary Technologies in Bioinformatics, by UNICAM ASSICOS and by NSF Faculty Early Career Award CCR0133583 and NSF Award CCF0523863.
This article has been published as part of BMC Bioinformatics Volume 9 Supplement 2, 2008: Italian Society of Bioinformatics (BITS): Annual Meeting 2007. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/9?issue=S2
References

Karma A: New paradigm for drug therapies of cardiac fibrillation.
Proc Natl Acad Sci U S A 2000, 97(11):56875689. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hodgkin AL, Huxley AF: A quantitative description of membrane currents and its application to conduction and excitation in nerve.
J Physiol 1952, 117:500544. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Barkley D, Knees M, Tuckerman L: Spiralwave dynamics in a simplemodel of excitable media  the transition from simple to compound rotation.

Di Francesco D, Noble D: A model of cardiac electrical activity incorporating ionic pumps and concentration changes.
Philos Trans R Soc Lond B Biol Sci 1985, 307:353398. PubMed Abstract

Luo CH, Rudy Y: A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes.
Circ Res 1994, 74(6):10711096. PubMed Abstract

Bub G, Glass L, Publicover NG, Shrier A: Bursting calcium rotors in cultured cardiac myocyte monolayers.
Proc Natl Acad Sci USA 1998, 95(17):1028310287. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gerhardt M, Schuster H, Tyson JJ: A cellular automation model of excitable media including curvature and dispersion.
Science 1990, 247(4950):15631566. PubMed Abstract  Publisher Full Text

Barkley D: A model for fast computersimulation of waves in excitable media.

Bub G, Shrier A: Propagation through heterogeneous substrates in simmple excitable media models.
Chaos 2002, 12(3):747753. PubMed Abstract  Publisher Full Text

Alur R, Courcoubetis C, Henzinger T, Ho P: Hybrid automata: an algorithmic approach to the specification and verification of hybrid systems.
In In Proceedings of Hybrid Systems I, LNCS 736 1993, 209229.

Henzinger TA: The theory of Hybrid Automata.
In Proceedings of the 11th IEEE Symposium on Logic in Computer Science 1996, 278293.

Ye P, Entcheva E, Smolka SA, Grosu R: Modelling excitable cells using cyclelinear hybrid automata.
IET Systems Biology 2008, 2(1):2432.
In Press
PubMed Abstract  Publisher Full Text 
True MR, Entcheva E, Smolka SA, Ye P, Grosu R: Efficient EventDriven Simulation of Excitable Hybrid Automata.
Engineering in Medicine and Biology Society, 2006. EMBS '06. 28th Annual International Conference on the IEEE 2006, 31503153.

Site of EHA [http://www.cs.sunysb.edu/~eha webcite].

Luo CH, Rudy Y: A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction.
Circ Res 1991, 68(6):15011526. PubMed Abstract

Berger R: Electrical restitution hysteresis  good memory or delayed response?
Circulation Research 2004, 94:567569. PubMed Abstract  Publisher Full Text