Volume 9 Supplement 2
Italian Society of Bioinformatics (BITS): Annual Meeting 2007
CellExcite: an efficient simulation environment for excitable cells
 Ezio Bartocci^{1, 3}Email author,
 Flavio Corradini^{1},
 Emilia Entcheva^{2},
 Radu Grosu^{3} and
 Scott A Smolka^{3}
DOI: 10.1186/147121059S2S3
© Bartocci et al.; licensee BioMed Central Ltd. 2008
Published: 26 March 2008
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/.
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 [3–5]. 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 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.
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
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
Performance comparison
Performance comparison for 2second simulation
cell array size  original  hybrid 

2 × 2 cell array  5 s  3 s 
4 × 4 cell array  9 s  3 s 
8 × 8 cell array  26 s  6 s 
16 × 16 cell array  93 s  14 s 
32 × 32 cell array  365 s  51 s 
64 × 64 cell array  1460 s  198 s 
400 × 400 cell array  17 h 10 m 33 s  2 h 13 m 38 s 
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
where V is the transmembrane voltage [mV], whose variation forms the AP; ${\overline{g}}_{Na},\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}{\overline{g}}_{K}\phantom{\rule{0.5em}{0ex}},\phantom{\rule{0.5em}{0ex}}{\overline{g}}_{L}$ 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.

As finite set X = x_{1}, ⋯, x_{ n } of realnumbered variables. The number n is called the dimension of H. We write $\dot{X}\phantom{\rule{0.5em}{0ex}}=\phantom{\rule{0.5em}{0ex}}{\dot{x}}_{1},\phantom{\rule{0.5em}{0ex}}\cdots ,\phantom{\rule{0.5em}{0ex}}{\dot{x}}_{n}$ 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 $X\cup \dot{X}$.

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 ${\overline{V}}_{S}$, 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
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.
HA for the NNR model
Parameters definitions
Parameters for HH LRd and NNR 4state HA models
V _{ R }  V _{ T }  V _{ O }  α ^{ 0 } _{ x }  α ^{ 0 } _{ y }  α ^{ 0 } _{ z }  α ^{ 1 } _{ x }  α ^{ 1 } _{ y }  α ^{ 1 } _{ z }  α ^{ 2 } _{ x }  α ^{ 2 } _{ y }  α ^{ 2 } _{ z }  α ^{ 3 } _{ x }  α ^{ 3 } _{ y }  α ^{ 0 } _{ z }  

HH  10  10  83  −0.98  −0.16  N/A  N/A  −0.16  N/A  1.4  15  N/A  −0.98  −0.16  N/A 
LRd  20  20  138  −0.1  −0.1  −0.1  N/A  −0.1  −0.1  200  0  100  −0.001  0.036  0.008 
NNR  20  30  120  −0.025  −0.07  −0.2  N/A  −0.07  −0.2  250  200  125  −0.025  −0.07  −0.2 
Diffusion Modeling in Multicellular HA
Weights calculated using neighborhood function for triangular lattice
#  0  1  2  3  4  5  6  7  8 

d  0  r  $\sqrt{3}r$  2r  $\sqrt{7}r$  3r  $2\sqrt{3}r$  $\sqrt{13}r$  4r 
$\stackrel{\Delta}{{\nu}_{r}}\left(d\right)$  −6  1  N/A  N/A  N/A  N/A  N/A  N/A  N/A 
${\stackrel{\Delta}{\nu}}_{2r}\left(d\right)$  −42.96  4.48  1.68  1  N/A  N/A  N/A  N/A  N/A 
${\stackrel{\Delta}{\nu}}_{3r}\phantom{\rule{0.5em}{0ex}}\left(d\right)$  −191.70  14.39  7.38  5.29  1.94  1  N/A  N/A  N/A 
${\stackrel{\Delta}{\nu}}_{2r}\left(d\right)$  −726.24  42.52  25.79  20.08  9.48  5.75  2.72  2.12  1 
Weights calculated using neighborhood function for square lattice
#  0  1  2  3  4  5  6  7  8  9 

d  0  r  $\sqrt{2}r$  2r  $\sqrt{5}r$  $2\sqrt{2}r$  3r  $\sqrt{10}r$  $\sqrt{13}r$  4r 
${\stackrel{\square}{\nu}}_{r}\left(d\right)$  −4  1  N/A  N/A  N/A  N/A  N/A  N/A  N/A  N/A 
${\stackrel{\square}{\nu}}_{2r}\left(d\right)$  −32.8  4:48  2.72  1  N/A  N/A  N/A  N/A  N/A  N/A 
${\stackrel{\square}{\nu}}_{3r}\left(d\right)$  −159.88  14:39  10.31  5.29  3.79  1.40  1  N/A  N/A  N/A 
${\stackrel{\square}{\nu}}_{4r}\left(d\right)$  −636.52  42:52  33.12  20.08  15.64  12.18  5.75  4.48  2.12  1 
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
Declarations
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
Authors’ Affiliations
References
 Karma A: New paradigm for drug therapies of cardiac fibrillation. Proc Natl Acad Sci U S A 2000,97(11):5687–5689.PubMed CentralView ArticlePubMedGoogle Scholar
 Hodgkin AL, Huxley AF: A quantitative description of membrane currents and its application to conduction and excitation in nerve. J Physiol 1952, 117: 500–544.PubMed CentralView ArticlePubMedGoogle Scholar
 Barkley D, Knees M, Tuckerman L: Spiralwave dynamics in a simplemodel of excitable media  the transition from simple to compound rotation. Physica Review A 1990,42(4):2489–2492.View ArticleGoogle Scholar
 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: 353–398.View ArticleGoogle Scholar
 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):1071–1096.View ArticlePubMedGoogle Scholar
 Bub G, Glass L, Publicover NG, Shrier A: Bursting calcium rotors in cultured cardiac myocyte monolayers. Proc Natl Acad Sci USA 1998,95(17):10283–10287.PubMed CentralView ArticlePubMedGoogle Scholar
 Gerhardt M, Schuster H, Tyson JJ: A cellular automation model of excitable media including curvature and dispersion. Science 1990,247(4950):1563–1566.View ArticlePubMedGoogle Scholar
 Barkley D: A model for fast computersimulation of waves in excitable media. Physica D 1991,49(1–2):61–70.View ArticleGoogle Scholar
 Bub G, Shrier A: Propagation through heterogeneous substrates in simmple excitable media models. Chaos 2002,12(3):747–753.View ArticlePubMedGoogle Scholar
 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, 209–229.View ArticleGoogle Scholar
 Henzinger TA: The theory of Hybrid Automata. In Proceedings of the 11th IEEE Symposium on Logic in Computer Science 1996, 278–293.View ArticleGoogle Scholar
 Ye P, Entcheva E, Smolka SA, Grosu R: Modelling excitable cells using cyclelinear hybrid automata. IET Systems Biology 2008,2(1):24–32. In PressView ArticlePubMedGoogle Scholar
 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, 3150–3153.View ArticleGoogle Scholar
 Site of EHA . [http://www.cs.sunysb.edu/~eha]
 Luo CH, Rudy Y: A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction. Circ Res 1991,68(6):1501–1526.View ArticlePubMedGoogle Scholar
 Berger R: Electrical restitution hysteresis  good memory or delayed response? Circulation Research 2004, 94: 567–569.View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.