Model Checking Implantable Cardioverter Defibrillators
Houssam Abbas, Kuk Jin Jang, Zhihao Jiang, Rahul Mangharam
Department of Electrical and Systems Engineering
University of Pennsylvania, Philadelphia, PA, USA

arXiv:1512.08083v1 [cs.SY] 26 Dec 2015

{habbas, jangkj, zhihaoj, rahulm}@seas.upenn.edu

ABSTRACT
Ventricular Fibrillation is a disorganized electrical excitation of the heart that results in inadequate blood flow to
the body. It usually ends in death within seconds. The
most common way to treat the symptoms of fibrillation is
to implant a medical device, known as an Implantable Cardioverter Defibrillator (ICD), in the patient’s body. Modelbased verification can supply rigorous proofs of safety and
efficacy. In this paper, we build a hybrid system model
of the human heart+ICD closed loop, and show it to be a
STORMED system, a class of o-minimal hybrid systems that
admit finite bisimulations. In general, it may not be possible to compute the bisimulation. We show that approximate
reachability can yield a finite simulation for STORMED systems, which improves on the existing verification procedure.
In the process, we show that certain compositions respect
the STORMED property. Thus it is possible to model check
important formal properties of ICDs in a closed loop with
the heart, such as delayed therapy, missed therapy, or inappropriately administered therapy. The results of this paper
are theoretical and motivate the creation of concrete model
checking procedures for STORMED systems.

1.

INTRODUCTION

Implantable Cardioverter Defibrillators (ICDs) are life-saving
medical devices. An ICD is implanted under the shoulder,
and connects directly to the heart muscle though two electrodes and continuously measures the heart’s rhythm (Fig.
1). If it detects a potentially fatal accelerated rhythm known
as Ventricular Tachycardia (VT), the ICD delivers a highenergy electric shock or sequence of pulses through the electrodes to reset the heart’s electrical activity. Without this
therapy, the VT can be fatal within seconds of onset. In the
US alone, 10,000 people receive an ICD every month. Studies have presented evidence that patients implanted with
ICDs have a mortality rate reduced by up to 31% [19].
Unfortunately, ICDs suffer from a high rate of inappropriate therapy due to poor detection of the current rhythm
on the part of the ICD. In particular, a class of rhythms

Sense

Therapy

Right Atrium
Electrode

ICD

Can (Shock)
Electrode
Left Atrium
Left Ventricle
Right Atrium
Right Ventricle

Shock
Coils

Right Ventricular Electrode
Atrial
Sensed Event (AS)
Ventricular
Sensed Event (VS)

AS

AS
VS

VS

Atrial
Signal
Ventricular
Signal
Shock
Signal

Figure 1: ICD connected to a human heart via two
electrodes. The ICD monitors three electrical signals
(known as electrograms) traversing the heart muscle.

known as SupraVentricular Tachycardias (SVTs) can fool
the detection algorithms. Inappropriate shocks increase patient stress, reduce their quality of life, and are linked to
increased morbidity [22]. Depending on the particular ICD
and its settings, the rates of inappropriate therapy can range
from 46% to 62% of all delivered therapy episodes [9]. Current practice for ICD verification relies heavily on testing
and software cycle reviews. With the advent of computer
models of the human heart, Model-Based Design (MBD) can
supply rigorous evidence of safety and efficacy. This paper
presents hybrid system models of the human heart and of
the common modules of ICDs currently on the market, and
shows that the closed loop formed by these models is formally verifiable. The objective is to develop model checkers
for ICDs to further their MBD process.
No work exists on ICD verification. Earlier work on verification of medical devices (formal or otherwise) focuses on
pacemakers. In [14] the authors developed timed automata
models of the whole heart+pacemaker loop which allows verification of LTL properties. In [6] the authors perform probabilistic testing of Hybrid I/O automata models of heart and
pacemaker. However, they can not be symbolically verified.
Later work on pacemakers [18] develops a formalized cellular automata (CA) model of the heart and uses Event-B for
expressing its properties, and in [12] invariants of pacemaker
and cardiac cells are verified. The ICD algorithms are more
complex than a pacemaker’s: an ICD measures the timing of
events, but also measures and processes the morphology of

Atrial Electrode

Sinoatrial
(SA) Node

Three Consecutive
Fast Intervals
(TCFI)

SA

Atria

Peak Tracking

Atrioventricular
(AV) Node

Vector Timing
Correlation
(VTC)

AV

Stability

Ventricles
Right Ventricular
Apex (RVA)

Dynamic Sensitivity
Adjustment
(Exponential Decay)

RVA
Right

Left
Ventricular Electrode

Human Heart

Blanking
Period

Input:
Electrogram
Waveforms

Cellular Automata Model

Output:
Event
Waveforms

ICD Sensing

Other
Discriminators

ICD
Detection

Figure 2: The whole heart is modeled as a 2D mesh of cells (Section 3). The ICD electrodes are shown in the right
atrium and ventricle. The electrogram signals measured through the electrodes are processed by the sensing module
(ICD Sensing, see Section 4). The detection algorithm (Section 5) determines the current rhythm using the processed
signal (ICD Detection).

the electrical signal in the heart to distinguish many types of
arrhythmias. Thus, we need three models for ICD verification: a timing and voltage model of the heart, a model of the
ICD’s algorithms, and a model for voltage measurement by
the ICD electrodes. This takes the model out of the realm
of timed automata and into hybrid automata proper. More
generally, approaches to approximate verification of similar
hybrid systems include falsification of general Metric Temporal Logic properties [5] and δ-reachability [16].
The first contribution of this paper is to develop a hybrid
system model of the heart, the ICD measurement process,
and of the algorithmic components of ICDs from most major manufacturers on the market (Fig. 2). We show that
the composition of these three models admits a finite bisimulation [1]. The ICD models presented here are the first
formalization of ICD operation to the best of our knowledge.
To establish this result we use the theory of STORMED hybrid systems [27], a class of hybrid systems that have finite
bisimulations. Our second contribution is two general results
for STORMED systems. First we prove that parallel compositions of STORMED systems yield STORMED systems.
Secondly, we show that any definable over-approximate reach
tubes can replace the exact trajectories of a STORMED system, yielding a system that still admits a finite simulation
(but no longer a bisimulation). Finally, we show that the
reach sets computed by the reachability tool SpaceEx [8] (a
widely used and scalable reachability tool) are definable and
so can be used to build the simulation. Thus SpaceEx can
be used as part of a model checker for STORMED systems.
Our interest in not simply in a particular manufacturer’s
arrhythmia detection algorithm: rather, we are interested
in those components that are common to most of them,
thus making our results relevant to them. The components we model or some variation on them are included in
the ICDs of Boston Scientific, Medtronic, Saint-Jude Medical and Biotronik. This is the first example of a practical
STORMED system that the authors are aware of.
Organization. Section 2 covers some preliminaries on hybrid systems. Sections 3 presents the heart model, and
Sections 4-5 model the ICD. Sections 6 and 7 prove general results on STORMED systems: namely that a defin-

Algorithm 1 Computing a bismimulation respecting ∼
Require: Transition system T = (Q, Σ, −
→, Q0 ), equivalence relation ∼.
Set S = Q/ ∼
while ∃P, P 0 ∈ S and σ ∈ Σ s.t. ∅ 6= P 0 ∩ P ostσ (P ) 6= P 0
do
Set S = S \ {P 0 } ∪ {P 0 ∩ P ostσ (P ), P 0 \ P ostσ (P )}
end while
Return S

able over-approximation of the flows such as that computed
by SpaceEx preserves finiteness of the simulation, and that
compositions of STORMED systems are STORMED.

2.

HYBRID SYSTEMS AND SIMULATIONS

This section presents fairly standard definitions on hybrid
systems and their simulations [1]. It also defines STORMED
hybrid systems, which admit finite bisimulations [27].

2.1

Transition and hybrid systems

Definition 2.1. A transition system T = (Q, Σ, −
→, Q0 )
consists of a set of states Q, a set of events Σ , a transition relation −
→⊂ Q × Σ × Q, a set of initial states Q0 . We
σ
write q −
→ q 0 to denote a transition element (q, σ, q 0 ) ∈−
→.
σ
→ q0 }
Given P ⊂ Q, we define P ostσ (P ) := {q 0 | ∃q ∈ P.q −
Given an equivalence relation ∼ on Q, the quotient system
∗
T / ∼ is T / ∼= (Q/ ∼, {∗}, −
→∼ , Q0 / ∼) where [q] −
→∼ [q 0 ]
σ
iff q −
→ q 0 for some σ ∈ Σ. Here [q] is the equivalence class
of q and Q/ ∼ is the set of equivalence classes of ∼.
Definition 2.2. Given two transition systems T1 and T2
with the same state space Q, a simulation relation from T1
to T2 is a relation S ⊂ Q × Q such that for all (q1 , q2 ) ∈
σ
σ
S, if q1 −
→1 q10 , there exists a q20 ∈ Q s.t. q2 −
→2 q20 and
0
0
(q1 , q2 ) ∈ S. A bisimulation relation between T1 and T2 is
both a simulation relation from T1 to T2 and from T2 to T1 .
The bisimulation B is said to respect ∼ if (q, q 0 ) ∈ B =⇒
q ∼ q 0 . The following algorithm, if it terminates, yields a
finite bisimulation for T that respects the given equivalence
relation [1]. Moreover, it is the coarsest bisimulation (with
respect to inclusion) that respects ∼. Given a set of atomic
propositions AP , if ∼ is s.t. q ∼ q 0 iff both states satisfy exactly the same set of atomic propositions, then model

checking temporal logic properties can be done on the finite
bisimulation instead of the possibly infinite T .
Definition 2.3. A hybrid automaton is a tuple
H = (X, L, H0 , {f` }, Inv, E, {Rij }(i,j)∈E , {Gij }(i,j)∈E )
where X ⊂ Rn is the continuous state space equipped with
the Euclidian norm k · k, L ⊂ N is a finite set of modes,
H0 ⊂ X × L is an initial set, {f` }`∈L determine the continuous evolutions with unique solutions, Inv : L → 2X defines
the invariants for every mode, E ⊂ L2 is a set of discrete
transitions, Gij ⊂ X is guard set for the transitions (so H
transitions i → j when x ∈ Gij ), Rij : X → X is an edgespecific reset function.
Set H = L × X. Given (`, x0 ) ∈ H, the flow θ` (; x0 ) : R+ →
Rn is the solution to the IVP ẋ(t) = f` (x(t)), x(0) = x0 .
The associated transition system is TH = (H, E ∪ {τ }, −
→
S
e
τ
e
, H0 ) with −
→= ( e∈E −
→)∪ −
→ where (i, x) −
→ (j, y) iff e =
τ
(i, j), x ∈ Gij , y = Rij (x) and (i, x) −
→ (j, y) iff i = j and
there exists a flow θi (·; x) of H and t ≥ 0 s.t. θi (t; x) = y
and ∀t0 ≤ t, θi (t0 ; x) ∈ Inv(i). For a set P ⊂ H,P|X denotes
its projection onto X, and P|L its projection onto L.
Definition 2.4. [Reachability] Let H be a hybrid system
with hybrid state space H, I = [0, b) ⊂ [0, +∞) be a (possibly
unbounded) interval, t ∈ I, and  > 0. The -approximate
continuous reachability operator, Rt : 2H → 2H is given by
Rt (P ) = {(i, x) ∈ X|∃x0 ∈ P|X , t ≥ 0.||θi (t; x0 ) − x|| ≤ }
where P = {i} × W , W ⊂ Inv(i). Define also RI (P ) =
∪t∈I Rt (P ). The (exact) discrete reachability operator is:
Rd (P )

= ∪j:(i,j)∈E Rij (P ∩ Gij )

For a hybrid system, P ostσ computes the forward reach sets,
and is implemented by R0[0,∞) and Rd . Algorithm 1, applied
to TH , implements the following iteration, in which Ft (P)
τ
→1 respecting
is the coarsest bisimulation with respect to −
e
the partition P, and Fd (P) := {(h1 , h2 ) | (h1 −
→ h01 ) =⇒
e0

(∃e0 ∈ E, h02 .h2 −→ h02 ∧ h01 ≡P h02 )} ∩ P [27]:
W0 = Ft (Q/ ∼),

∀i ≥ 0, Wi+1 = Ft (Fd (Wi ))

(1)

This iteration (equivalently, Alg. 1) does not necessarily terminate for hybrid systems because the reach set might intersect a given block of Q/ ∼ an infinite number of times
(see [17] for an example). The class of systems introduced
in the next section has the property that Algorithm 1 does
terminate for it and returns a finite S.

2.2

O-minimality and STORMED systems

We give a very brief introduction to o-minimal structures.
A more detailed introduction can be found in [17] and references therein. We are interested in sets and functions
in Rn that enjoy certain finiteness properties, called orderminimal sets (o-minimal). These are defined inside structures A = (R, <, +, −, ·, exp, . . .). The subsets Y ⊂ Rn
we are interested in are those that are definable using firstorder formulas ϕ: Y = {(a1 , . . . , an ) ∈ Rn | ϕ(a1 , . . . , an )}.
(First-order formulas use the boolean connectives and the
1

I.e., Ft only considers the continuous transition relation.
τ
c
Namely, it is a bisimulation of TH
:= (Q/ ∼, {∗}, −
→, Q0 / ∼).

quantifiers ∃, ∀). The atomic propositions from which the
formulas are recursively built allow only the operations of
the structure A on the real variables and constants, and the
relations of A and equality. For example 2x − 3.6y < 3z
and x = y are valid atomic propositions of the structure
LR = (R, <, +, −, ·), while cosh(x) < 3z is not because cosh
is not in the structure. These structures are already sufficient to describe a set of dynamics rich enough for our
purposes and for various classes of linear systems.
Definition 2.5. A theory of (R, . . .) is o-minimal if the
only definable subsets of R are finite unions of points and
(possibly unbounded) intervals. A function f : x 7→ f (x) is
o-minimal if its graph {(x, y) | y = f (x)} is a definable set.
We use the terms o-minimal and definable interchangeably,
and they refer to Lexp = (R, <, +, −, ·, exp) which is known
to be o-minimal. The dot product between x, y ∈ Rn is
denoted x · y, and d(Y, S) = inf{ky − sk | (y, s) ∈ Y × S}.
Definition 2.6. [27]. A STORMED hybrid system (SHS)
Σ is a tuple (H, A, φ, b− , b+ , dmin , , ζ) where H is a hybrid
automaton, A is an o-minimal structure, dmin , , ζ are positive reals, b− , b+ ∈ R and φ ∈ X such that:
(S) The system is dmin -separable, meaning that for any
e = (`, `0 ) ∈ E and `00 6= `0 ,d(Re (G(`,`0 ) ), G(`0 ,`00 ) ) > dmin 2
(T) The flows (i.e., the solutions of the ODEs) are TimeIndependent with the Semi-Group property (TISG), meaning
that for any ` ∈ L, x ∈ X, the flow θ` starting at (`, x) satisfies: 1) θ` (0; x) = x, 2) for every t, t0 ≥ 0, θ` (t + t0 ; x) =
θ` (t0 ; θ` (t; x))
(O) All the sets and functions of H are definable in the ominimal structure A
(RM) The resets and flows are monotonic with respect to
the same vector φ, meaning that
1) (Flow monotonicity) for all ` ∈ L, x ∈ X and t, τ ≥ 0,
φ · (θ` (t + τ ; x) − θ` (t; x)) ≥ ||θ` (t + τ ; x) − θ` (t; x)||, and
2) (Reset monotonicity) for any edge (`, `0 ) ∈ E and any
x− , x+ ∈ X s.t. x+ = R`,`0 (x− ),
1. if ` = `0 , then either x− = x+ or φ · (x+ − x− ) ≥ ζ
2. if ` 6= `0 , then φ · (x+ − x− ) ≥ ||x+ − x− ||
(ED) Ends are Delimited: for all e ∈ E we have φ · x ∈
(b− , b+ ) for all x ∈ Ge
Intuitively, the above conditions imply the trajectories of the
system always move a minimum distance along φ whether
flowing or jumping, which guarantees that no area of the
state space will be visited infinitely often. This is at the
root of the finiteness properties of STORMED systems. The
following result justifies the interest in STORMED systems:
they admit finite bisimulations.
Theorem 2.1. [27] Let H be a STORMED hybrid system,
2
The original definition of separability [27] required the
guards themselves to be separated, which is insufficient to
guarantee that if H flows, it flows a uniform minimum distance along φ. Indeed assume the guards are separated. If
x ∈ G(`,`0 ) and y = R(`,`0 ) (x), it can be that y ∈ G(`0 ,`00 )
and thus a jump happens, even though G(`,`0 ) and G(`0 ,`00 )
are separated. Therefore we need d(y, G`0 ,`00 ) > dmin for
all y ∈ Re (Ge ), which is the condition we use in Def. 2.6.
The properties of SHS, in particular the existence of finite
bisimulation, are therefore preserved by this change.

End

Phase 4

Phase 0

V̇ (i, j) = a(i, j)| V

V T /SV T ?

ṫ = 0

ṫ = 1
t˙p = 0
V (i, j)  Vmin ?

V (i, j) > Vth ?

Phase 1

V̇ (i, j) = d (d > 0)
ṫ = 1
t˙p = 0

tp

t

g, (g > 0)

ṫ = 1
t˙p = 0
t

Phase 2

Phase 3 - ERP
V (i, j) =

b, (b > 0)

ṫ = 1

Figure 3: When the ICD makes a VT/SVT decision,
all systems transition to mode End.

V (i, j) > Vmax ?

V̇ (i, j) =

t

tp > PD?

tp = 0

Proposition 2.1. If the state space X of a hybrid automaton H is bounded, then its guards have delimited ends.
Proof. For all guard sets G and all x ∈ G, ||φ · x|| ≤ ||φ|| ·
||x|| ≤ ||φ||. max{||x||, x ∈ X} < ∞.

3.

HEART MODEL

For the verification of ICDs, we adopt the cellular automata
(CA)-based heart model developed in [24],[7]. This model
lies in-between high spatial fidelity but slow to compute
PDE-based whole heart models [26], and low spatial fidelity
but very fast-to-compute automata-based models [20]. PDEbased models are not currently amenable to formal verification, both theoretically and practically. Models based on
ionic currents [13] might be more accurate but are likely to
be more computationally expensive. Timed automata models can not simulate the electrograms needed for ICD verification. CA-based models are appealing due to their intuitive
correspondence with the heart’s anatomy and function and
their relative computational simplicity. CA-based models
were used in [18],[2] and [6]. This paper’s model also has the
important advantage of forming the basis of software used
to train electrophysiologists, and allows interactive simulation of surgical procedures like ablation [23]. In particular,
it can simulate fibrillation and other tachycardias.
This paper’s automata:All hybrid automata in this paper
have the whole state space as invariants and transitions are
urgent (taken immediately when the guard is enabled). We
also observe that, as will be seen in Section 5, i) the ICD
will always reach a decision of VT or SVT in finite time, ii)
at which point it resets its controlled (software) variables so
new values are computed for the next arrhythmia episode.
So while the heart can beat indefinitely, for the purposes
of ICD verification, there’s a uniform upper bound on the
length of time of any execution. Let D ≥ 0 be this duration
(D is on the order of 30sec depending on device settings).
Also, the electrogram (EGM) voltage signal s has upper and
lower bounds s and s. Therefore, every mode of every automaton in what follows has a transition to mode End shown
in Fig. 3. We don’t show these transitions in the automata
figures to avoid congestion.

3.1

Cellular automata model

The heart has two upper chambers called the atria and two
lower chambers called the ventricles (Fig. 1) The synchronized contractions of the heart are driven by electrical activity. Under normal conditions, the SinoAtrial (SA) node (a
tissue in the right atrium) spontaneously depolarizes, pro-

tp

Vmax,2 ?
t

V̇ (i, j) = d2 , (d2 < d)
V (i, j)

Vth,2 ?

tp

t

ṫ = 1
t˙p = 0

Upstroke 2

Phase 3 - RRP

We need the following result in what follows.

t

ṫ = 1
t˙p = 0
V (i, j)

V̇ (i, j) = a(i, j)| V
ṫ = 1
t˙p = 0

tp

V̇ (i, j) = 0

V (i, j)  Vth

and let P be an o-minimal partition of its hybrid state space.
Then H admits a finite bisimulation that respects P.

tp > DP h1 ?

Figure 4: Hybrid model Hc of one cell of the heart
model. AP figure from [11]. Vth,2 > Vth , Vmax,2 < Vmax
ducing an electrical wave that propagates to the atria and
then down to the ventricles (Fig.2) In this model, the myocardium (heart’s muscle) is treated as a 2D surface (so it
has no depth), and discretized into cells, which are simply
regions of the myocardium (Fig. 2). Thus we end up with
N 2 cells in a square N -by-N grid. A cell’s voltage changes
in reaction to current flow from neighboring cells, and in response to its own ion movements across the cell membrane.
This results in an Action Potential (AP).
Fig. 4 shows how the AP is generated by a given cell [15]:
in its quiescent mode (Phase 4), a cell (i, j) in the grid has
a cross-membrane voltage V (i, j, t) equal to Vmin < 0. As it
gathers charge, V (i, j, t) increases until it exceeds a threshold voltage Vth . In Phase 0, the voltage then experiences a
very fast increase (Phase 0), called the upstroke, to a level
Vmax > 0, after which it decreases (Phase 1) to a plateau
(Phase 2). It stays at the plateau level for a certain amount
of time PD then decreases linearly to below Vth (Phase 3 ERP). Once below Vth it is said to be in the Relative Refractory Period (Phase 3 - RRP) . In Phase 3 - RRP, the cell
can be depolarized a second time, albeit at a higher threshold Vth,2 , slower and to a lower plateau level Vmax,2 < Vmax
(Upstroke 2). Otherwise, when the voltage reaches Vmin
again, the cell enters the quiescent stage again. This model
is suitable for both pacemaker and non-pacemaker cells, the
main differences being in the duration of the plateau (virtually non-existent for pacemaker cells), and the duration of
phases 0 and 4 (both are shorter for pacemaker cells).
In Fig. 4, V (i, j) ∈ R denotes the voltage in cell (i, j) of
2
the grid, and V = (V (1, 1), . . . , V (N 2 , N 2 ))T in RN groups
the cross-membrane voltages of all cells in the heart. The
whole heart model HCA is the parallel composition of these
N 2 single-cell models. The (i, j)th cell’s voltage at time t
in Phase 4 depends on that of its neighbors and its own as
follows [24]
V̇ (i, j, t)

=

1
[V (i − 1, j, t) + V (i + 1, j, t) − 2V (i, j, t)]
Rh
1
+
[V (i, j − 1, t) + V (i, j + 1, t) − 2V (i, j, t)]
Rv

= a(i, j)T V (t), a(i, j) ∈ RN

2

(2)

where Rh , Rv are conduction constants that can vary across
the myocardium. Thus V evolves according to a linear ODE

ẏ == 0 ^ ÿ < 0?

V̇ = AV where A is the matrix whose rows are the a(i, j).
The two states t and tp are clocks. Clock tp keeps track of
the value of the last discrete jump. We will use this arrangement in all our models: it avoids resetting the clocks which
preserves Reset Monotonicity.

Peak Tracking

y

f

1

ṫ = 1
t˙p = 0

(y

where k·k is the Euclidian norm, p0 and p1 are the electrodes’
positions and pi,j is the position of the (i, j)th cell on the
2D myocardium (p0 , p1 , pi,j ∈ R2 ). Positions p0 , p1 should
be chosen different from pi,j to avoid infinities.

T h, T h0

T h)^
(t

eF

tp

M inDecP )?

tp

t

yM

0

f

0

(

t, tp

0

Blanking
ṫ = 1
t˙p = 0
T˙h = 0

M inT P ^ f == 1?
minT h
eF
( 1/3)⇤ln{
}
Th
T h, T ho
(3/4) ⇤ yM
t

T˙h = 0

ICDs observe the electrical activity through three channels
(Fig. 1). Each signal is called an electrogram (EGM) signal.
The signal read on a channel is given by [7]:


1
1
1 X
−
V̇ (i, j, t) (3)
s(t) =
K i,j ||pi,j − p0 ||
||pi,j − p1 ||

tp

tp

t

t

tp

blankingP eriod ?
tp

t

ṫ = 1
t˙p = 0

1
minT h
) ⇤ ln{
} T h = max{minT h, decT h}
3
Th
minT h
Exponential Decay
eF
⇤ t}
TC
˙
y := |s|; eF = 0; T˙h = 0

decT h := T ho ⇤ exp{

Figure 5: HSense . States not shown in a mode have a 0
˙ = 0 in all modes.
derivative, e.g., eF
Ventricular Signal

500

Signal
Event
Threshold
Threshold Minimum
Blanking Period

450
400
350

Amplitude (mV)

Extensions. The Action Potential Duration (APD) restitution mechanism of heart cells as modeled in [24] can be
included in this model without changing its formal properties. More detailed APD restitution models exist [10]. Also,
note that cell topology (the way cells are connected to each
other) is not a factor in determining the STORMED property, so other topologies than a rectangular mesh may be
used.

yM

300
250

Blanking
Period

200

Threshold

150

We now state and prove the main result of this section.

Minimum Threshold

100
50

Theorem 3.1. Let HCA be the whole heart cellular automaton model obtained by parallel composition of N 2 models Hc
2
with state vector x = [V, t, tp , s] ∈ RN × R3 . Assume that
all executions of the system have a duration of D ≥ 0. Then
HCA is STORMED.
Proof. We verify each property of STORMED. In this and
all the proofs that follow, the approach is the same: (ED)
holds by Prop. 2.1 because our state spaces are bounded. After establishing properties (S), (T ) and (O), we draw up the
constraints on φ and ε imposed by reset and flow monotonicity (property (RM)). Then we argue that these constraints
can be solved for φ and ε. Often there is more than one
solution and we just point to one.
(S) Separability holds because Vmin < Vth < Vth,2 < Vmax,2 <
Vmax and P D > 0, DP h1 > 0. For example, on transition
Phase 4 → Phase 0, V (i, j) = Vth , which is separated from
the next guard {V (i, j) > Vmax } by |Vmax − Vth |.
(T) All flows are linear or exponential and thus are TISG.
(O) The flows, resets and guard sets are all definable in Lexp .
In particular the flow of V̇ = AV is exponential with real
exponent, and s is a sum of exponentials and linear terms.
2
(RM) We seek a vector φ = (φV , φt , φp , φs )T ∈ RN +3 such
that resets and flows are monotonic along φ. Only transitions p → q 6= p are to be found in HCA , during which only
−
tp is reset. Always, t+
p = t ≥ tp , thus the reset is indeed
monotonic as can be seen by choosing any ε > 0 and φp > ε.
Monotonic flows: φ must also be such that in all modes:
φ · (θ` (t + τ ; x) − θ` (t; x)) ≥ ε||θ` (t + τ ; x) − θ` (t; x)||
Decomposing, we want
φV · (V (t + τ ) − V (t)) + φt τ + φp · 0
(4)
+φs · (s(x, t + τ ) − s(x, t)) ≥ ε||θ` (x, t + τ ) − θ` (x, t)||

0

9100

9200

9300

9400

9500

9600

9700

9800

9900

10000

Seconds (ms)

Figure 6: Example of dynamic threshold adjustment in
ICD sensing algorithm. The shown signal is rectified.

Now note that all flows have bounded derivatives in every
bounded duration of flow and are thus Lipschitz. Let LV be
the Lipshitz constant of V (t) and Ls that of s(t). Then on
the LHS of the above inequality we have φV · (V (t + τ ) −
V (t)) + φs · (s(t + τ ) − s(t)) ≥ −φV LV τ − φs Ls τ . On the
RHS we have ε(LV τ + Ls τ + τ ) ≥ ε(||V (t + τ ) − V (t)|| +
||s(t + τ ) − s(t)|| + τ ) ≥ ε(||θ` (x, t + τ ) − θ` (x, t)||) Thus (4)
is satisfied if the stronger inequality
−φV LV τ − φs Ls τ + φt τ ≥ ε(LV τ + Ls τ + τ )
is satisfied. But this can be achieved by, for example, choosing φV = φs = 0 and φt ≥ ε(LV + Ls + 1).
(ED) Our system has bounded state spaces: V and s are
voltages typically in the range [−80, 60] mV and tp ≤ t ≤ D.
So (ED) holds by Lemma 2.1.

4.

ICD SENSING

Sensing is the process by which cardiac signals s measured
through the leads of the ICD are converted to cardiac timing events. The ICD sensing algorithm is a threshold-based
algorithm which declares events when the signal exceeds a
dynamically-adjusted threshold T h.
Fig. 5 shows the model HSense of the sensing algorithm,
and Fig. 6 illustrates its operation. The sensing takes place
on the rectified EGM signal y = |s|. After an event is
declared at the current threshold value (y(t) ≥ T h(t) in
Fig. 5), the algorithm tracks the signal in order to measure
the next peak’s amplitude (Peak Tracking). For a duration
M inT P (min tracking period) the latest peak is saved in
yM . A variable f indicates that a peak was found. Af-

ter a peak is found (f == 1) and after the end of the
tracking period, the algorithm enters a fixed Blanking Period (Blanking), during which additional events are ignored.
On the transition to Blanking, T h and T h0 are set to 3/4
the current value of yM and the exponential factor of deh
). At the end of
cay is updated (eF = (−1/3) ∗ ln minT
TH
the blanking period, the algorithm then transitions to the
Exponential Decay mode in which T h decays exponentially
from T h0 to a minimum level (Exponential Decay): T h(t) =
max(minT h, T h0 · exp(−(eF/T C)t)). The algorithm stays
in the Exponential Decay mode for at least a sampling period of M inDecP . Correspondingly, there is a de facto Maximum Decay Period M axDecP after which the system transitions again to PeakTracking since the signal y is bound
to exceed the minimum threshold minT h. Different manufacturers may use a step-wise decay instead of exponential,
but the principle is the same. Local peak detection is modeled via the ẏ = 0 ∧ ÿ < 0 transition. While y = |s| is
non-differentiable at 0, the peak will occur away from 0, as
shown in Fig. 6. The other states in Fig. 5 are t, tp (clocks).
minT h and T C are constant parameters.
Theorem 4.1. HSense is STORMED.
Proof. (S) By definition, we only need to consider transitions between different modes to establish separability. For
all such transitions, there is a minimum dwell time in the
mode before taking the transition, namely M inT P in PeakTracking, BlankingP eriod in Blanking, and M inDecP in
mode ExponentialDecay. So the system is separable since
there is a uniform minimum flow before jumping.
(T) Flows are either constant, (piece-wise) linear, or piecewise linear and exponential (in the case of y and its derivatives) and therefore are TISG.
(O) All the flows, resets and guard sets are definable in
Lexp . (The absolute value and max functions can be broken
down into boolean disjunctions of definable functions, and
t 7→ ln(t) is o-minimal by o-minimality of exp).
(RM) The state is x = (t, tp , y, yM , f, T h, T h0 , eF ) ∈ R8 ,
and let φ = (φt , φp , φy , φm , φf , φT h , φ0 , φeF ) be the corresponding φ vector. Recall that the EGM voltage s, and so
y = |s|, is upper-bounded by VM .
ExponentialDecay → PeakTracking. Only tp , yM and
f are modified, so monotonicity produces the constraint
W ant

φp (t−tp )+φm (0−yM )+φf (0−1) ≥ ε(|t−tp |+|yM |+1).
We require the stronger constraint to hold:
φt M inDecP − φm VM − φf

W ant

≥ ε(M axDecP + VM + 1)

PeakTracking → PeakTracking. Only yM and f are
W ant

reset. Algebraic manipulation yields −2VM φm + φf ≥ ζ
PeakTracking → Blanking. tp , eF, T h and T h0 are reset,
so we get
φp (t − tp ) + φeF (−(1/3) ln(minT h/T h) − eF )
+φT h (3yM /4 − T h) + φ0 (3yM /4 − T h0 )
minT h
1
) − eF |
≥ ε(|t − tp | + | − ln(
3
Th
3yM
3yM
+|
− T h| + |
− T h0 |)
4
4
T h is lower-bounded by minT h at all times, and it is naturally upper-bounded by VM as the threshold should never

exceed the largest possible attainable voltage. By the same
token, 0 ≤ eF ≤ (1/3) ln(VM /minT h). Then we want the
stronger inequality
φp M inT P

+ φeF (0 − (1/3) ln(VM /minT h)
+ φT h (−VM ) + φ0 (−VM )
1
VM
≥ ε(M axT P + | ln(
)| + |VM | + |VM |)
3
Th

Blanking → ExponentialDecay. Only tp is reset and
therefore we want, φp (t−tp ) ≥ ε(|t−tp |), thus the transition
yields φp ≥ ε.
The above equations can be simultaneously satisfied. The
simplest thing would be to set all φ terms that appear above
to 0 except for φt , φp which are calculated accordingly.
The flows can be shown to be monotonic along the same φ
and with the same ε. For example, in mode ExponentialDecay, only t, y and T h flow. Making use of the VM bound on y,
we get the constraint φt τ −2VM φy +φT h (T h(t+τ )−T h(t)) ≥
ε(τ +2VM +|T h(t+τ )−T h(t)|), which yields φt ≥ ε, φy ≤ −ε
and φT h ≥ ε. Similarly for the rest.

5.

ARRHYTHMIA DETECTION

Ventricular Tachycardia (VT) is an example of a tachycardia
originating in the ventricles, in which the ventricles spontaneously beat at a very high rate. If the VT is sustained,
or degenerates into Ventricular Fibrillation (VF), it can be
fatal. A tachycardia that originates above the ventricles is
referred to as a SupraVentricular Tachycardia (SVT) and is
a diseased but non-fatal condition. In what follows, we will
refer to sustained VT and VF together as VT. The ICD’s
main task is to discriminate VT from SVT and deliver therapy to the former only.
Most VT/SVT detection algorithms found in ICDs today
are composed of individual discriminators. A discriminator is a software function whose task is to decide whether
the current arrhythmia is SVT or VT. No one discriminator can fully distinguish between SVT and VT. Thus a
detection algorithm is often a decision tree built using a
number of discriminators running in parallel. The detection algorithm of Boston Scientific is shown in Fig. 7 [3].
We have modeled each discriminator in this detection algorithm as a STORMED hybrid system. The algorithm
itself is then a hybrid system. The ICD system is thus
HICD = HSense ||HDetection−Algo where HDetection−Algo
is the parallel composition of the discriminator systems. In what follows, we present three of these discriminators we modeled, which are found in most ICDs and model
them as hybrid systems, and prove they are STORMED.

5.1

Three Consecutive Fast Intervals

Our first module simply detects whether three consecutive
fast intervals have occurred, where ‘fast’ means the interval
length, measured between 2 consecutive peaks on the EGM
signal, is shorter than some pre-set amount. See Fig. 8.
States t and tp are clocks as before. The vector L3 is threedimensional, and stores the values of the last three intervals.
The event VEvent? is shorthand for the transition y(t) ≥ T h
being taken by the HSense automaton. In other words, it
indicates a ventricular event. Then L3 gets reset to L+
3 =

No

8/10 intervals faster
than VT threshold
!

Atrial Origin
Electrogram

Initial Detection

Last 10
Ventricular
intervals

NSR Template

No

8/10 intervals faster
than VF threshold
!

Yes

Ventricular Origin
Electrogram

Yes

VT Duration

VT#

VF Duration

RhythmID
No

V rate>A rate
by at least 10bpm

No

Yes

VT#

No

NSR Template

Yes

VTC
Correlated

SVT#

A-fib Rate=TRUE
& V rate unstable

Figure 9: EGMs of different origin have different mor-

Yes

phologies. The correlation of an EGM with respect to a
stored EGM template is used to determine the origin.

Figure 7: Boston Scientific’s detection algorithm

~⌫

µ

µ + s(t)

Circulate(~⌫ , 1)

↵

↵ + s(t)sm (t)

0, w

tp

w
1
R2i

R1

t

ṫ = 0
µ̇ = ↵˙ = ˙ = 0

ṫ = 1
µ̇ = ↵˙ = ˙ = 0

Circulate(L3 , t)

Figure 8: Three Consecutive Fast Intervals HT CF I

⇢˙ = 0
ẇ =

(z1 , z2 , z3 )+ := Circulate(L3 , t − tp ) where

 



0 1 0
0
z2
 z3  = 0 0 1 L3 +  0 
L+
3 =
0 0 0
t − tp
t − tp

+ s(t)2

1

V Event?
L3

R3

, ( ⌧ 1)

L3  th?
DurationEnds?

Calculate VTC

⇢˙ = 0
ẇ = 0

IDLE

W indowEnds & ⇢new < th?
~⌫

(5)

Lemma 5.1. HT CF I is STORMED.
Proof. We show that the reset are monotonic - the other
properties are easily checked. For reset monotonicity, we
invoke the fact that there is a minimum beat-to-beat separation: heartbeats can’t follow one another with vanishingly
small delays. In other words, there exists m > 0 such that
t − t−
p > m. Similarly, there’s a maximum delay between
two heartbeats, call it B. Now, we seek a vector φ ∈ R5 s.t.




t−t
z2 − z1
W ant
φ·  t − tp  = φp (t−tp )+φL3 ·  z3 − z2  ≥ ζ > 0
+
t − t p − z3
L3 − L3
{z
}
|
δ

(6)
p
Now |δ| is upper bounded by 3 · (2B)2 since each element
is the difference of intervals shorter than B. Also, t − t−
p >
m > 0. So choose φL3 = (φz,1 , φz,2 , φz,3 ) > 0 elementwise. (6) is satisfied if the following stronger inequality is
satisfied, which can√be achieved by an appropriate choice of
P
φz,i : φp m ≥ ζ + 12B 2 31 φz,i

5.2

t == i · Ts ?

th?

t, µ, ↵,

ṫ = 1
t˙p = 0
L˙3 = 0

W indowEnds & ⇢new

Vector Timing Correlation

It has been clinically observed that a depolarization wave
originating in the ventricles (as produced during VT for example) will in general produce a different EGM morphology
than a wave originating in the atria (as produced during
SVT) [3]. See Fig. 9. A morphology discriminator measures the correlation between the morphology of the current
EGM and that of a stored template EGM acquired during
normal sinus rhythm. If the correlation is above a pre-set
threshold for a minimum number of beats, then this is an
indication that the current arrhythmia is supraventricular
in origin. Otherwise, it might be of ventricular origin.

t, µ, ↵,

Circulate(~⌫ , 1)
0, w

1

Figure 10: VTC calculation. iTs is the sampling time
for the ith fiducial point, i = 1, . . . , 8. R21 , . . . , R28 are
the corresponding resets. For clarity of the figure, 8
transitions are represented on the same edge.
Boston Scientific’s implementation of a morphology discriminator is called Vector and Timing Correlation (VTC). VTC
first samples 8 fiducial points si , i = 1, . . . , 8 on the current
EGM s at pre-defined time instants. Let sm,i be the corresponding points on the template EGM. The correlation is
then calculated as [3]
P
P
P
(8 i si sm,i − ( i si )( i sm,i ))2
P
P
P
P
ρnew =
(8 i s2i − ( i si )2 )(8 i s2m,i − ( i sm,i )2 )
Note that sm is a constant for the purposes of this calculation: it does not change during an execution of VTC. If
3 out of the last 10 calculated correlation values exceed the
threshold, then SVT is decided and therapy is withheld.
The system of Fig. 10 implements the VTC discriminator.
As before, t is a local clock. µ accumulates the values of the
current EGM, α accumulates the product si sm,i , β accumulates s2i . State w is an auxiliary state we need to establish
the STORMED property. ~ν is a 10D binary vector: νi = −1
if the ith correlation value fell below the threshold, and is
+1 otherwise. L3 is the state of HT CF I : the guard condition L3 ≤ th indicates that all its entries have values less
than the tachycardia threshold, which is when HV T C starts
computing. W indowEnds indicates the ‘end’ of an EGM,
measured as a window around the peak sensed by HSense .
Lemma 5.2. HV T C is STORMED.
Proof. Separability obtains by observing that a uniform
minimum time passes between beats and between samples.
TISG is immediate. O-minimality is established by observing that all sets and functions are definable in Lexp . ED

Idle

holds because the state space is bounded. We now show
monotonicity. The state of the system is x = (t, µ, α, β, ~ν , w)T ∈
R4+10+1 . Let φ = (φc , φµ , φα , φβ , φ1 , . . . , φ10 , φw )T ∈ R15 be
the corresponding vector. For flows in mode CalculateVTC,
we seek a φ and ε > 0 such
p that φ · (t + τ − t, 0, −γ(t + τ ) +
γt) = φc τ +φw (−γτ ) ≥ ε τ 2 + γ 2 τ 2 , which is equivalent to
p
φc − φw γ ≥ ε 1 + γ 2 . Reset monotonicity for resets R1,
R2, R3 provides three more constraints on φ and ε:

L˙2 = L˙1 = ̇ =
2

Finalize
ṫ = 1
L˙2 = L˙1 = ̇ =
2

=

−φc t − φµ µ − φα α − φβ β +
W ant

+φw (1 − w) ≥ ζ
(R2) φ · (t − t, s, ssm , s2 , 0, 1 − w)
W ant

φµ s + φα ssm + φβ s2 + φw (1 − w) ≥ ζ

(R3) −φc t − φµ µ − φα α − φβ β +

10
X

φi (νi+1 − νi )

i=1
W ant

+φw (1 − w) ≥ ζ
where ν11 := −1 in R1 and ν11 := 1 in R3. Combine R1
and R3 by choosing φ1 = . . . = φ10 = φµ = φα = φβ = 0:
(R1, 3) − φc t + φw (1 − w) ≥ ζ
(R2) φw (1 − w) ≥ ζ
Now note that when a reset occurs, 0 < w ≤ 1 − γTs := wm
where Ts is the smallest sampling period, and that t ≤ 10B,
B = the maximum peak-to-peak interval, so (R2), (R1, 3)
can be jointly satisfied if −φc 10B + φw (1 − wm ) ≥ ζ . The
2 boxed equations can be jointly satisfied.

5.3

ṫ = 1
L˙2 = L˙1 = ̇ =
2

=0

V Event?
t

0

L2

L 2 + t2

L1

L1 + t
+1

DurationEnds?
1
[L2 L21 /]
2


=0

Figure 11: Stability discriminator.

φi (νi+1 − νi )

i=1

=

=0

DurationBegins?



(R1) φ · (−t, −µ, −α, −β, ν2 − ν1 , ν3 − ν2 , . . . , −1 − ν10 , 1 − w)
10
X

Accumulate

ṫ = 1

Stability discrimination

Stability refers to the variability of the peak-to-peak cycle
length. A rhythm with large variability (above a pre-defined
threshold) is said to be unstable, and is called stable otherwise. The Stability discriminator is used to distinguish between atrial fibrillation, which is usually unstable, and VT,
which is usually stable.
The Stability discriminator shown in Fig. 11 simply calculates the variance of the cycle length over a fixed period called a Duration (measured in seconds). Let DL ≥ 0
be the Duration length. The events DurationBegins? and
DurationEnds? indicate the transitions of a simple system
that measures the lapse of one Duration (not shown here).
State t is a clock, L1 accumulates the sum of interval lengths
(and will be used to compute the average length), L2 accumulates the squares of interval lengths, and κ is a counter
that counts the number of accumulated beats. σ2 is assigned
the value of the variance given by κ1 [L2 − L21 /κ]
Lemma 5.3. HStab is STORMED.
The proof is in the Appendix.
Now that each system was shown to be STORMED, it remains to establish that their parallel composition is STORMED.
This result does not hold in general - Thm. 6.1 gives conditions under which parallel composition respects the STORMED

property. Intuitively, we require that whenever a sub-collection
of the systems jumps, the remaining systems that did not
jump are separated from all of their respective guards by a
uniform distance. This is a requirement that can be shown
to hold for our systems by modeling various minimal delays
in the systems’ operation. We may now state:
Theorem 5.1. Consider the collection of systems HCA ,
HICD = HSense ||HDetection−Algo where the latter is the parallel composition of the discriminator systems. This collection satisfies the hypotheses of Thm. 6.1 (Section 6) and
therefore the parallel system HCA ||HICD is STORMED and
has a finite bisimulation.

6.

COMPOSING STORMED SYSTEMS

The results in this section and the next apply to STORMED
systems in general, including those with time-unbounded
operation. We write [m] = {1, . . . , m}. Given hybrid systems H1 , . . . , Hm in this section, xi , Gi , θi , . . . etc refer to
a state, guard, flow . . . of system Hi , i ≤ m. We show
that the parallel composition of SHS is still a SHS. Recall that θ` (t; x) is the flow starting at (`, x). Given hybrid systems H1 , . . . , Hm , their parallel composition H =
H1 || . . . ||Hm is defined in the usual way: H.X = Πi X i ,
H.L = Πi Li , H.H0 = Πi H0i , Inv(`) = Πi Inv i (`i ), θ` (x, t) =
[θ`11 (x1 , t)(t), . . . , θ`mm (xm , t)(t)]T . The system jumps if any
of its subsystems jumps, so its guard sets are of the form
A1 × . . . × Am where for at least one i, Ai is a guard of Hi ,
and for the rest Aj = X j . When a guard of a subsystem
is satisfied, the state of that subsystem is reset according to
its reset map. The guards are made disjoint to avoid nondeterminism. A system H is deterministic if to every initial
state (`, x), H produces a unique trajectory starting there.
In general H is not separable: indeed for any candidate value
of dmin , one could find a transition (i, j) of H due to, say, a
jump of H1 , s.t. at that moment x2 is closer than dmin to
one of its own guards, say G2(j 2 ,k2 ) . This causes H to further
jump j → k without having traveled the requisite minimum
distance, thus violating the separability of Rij (Gij ) and Gjk .
Therefore we need to impose an extra condition on minimum
separability across sub-systems.
Theorem 6.1. Let Σi = (Hi , A, φi , bi,− , bi,+ , dimin , εi , ζ i ),
i = 1, . . . , m be deterministic SHS defined using the same
underlying o-minimal structure, and where each state space
X i is bounded by BX i .
Define parallel composition Σ = (H, A, φ, b− , b+ , dmin , ε, ζ)
where H = H1 || . . . ||Hm , φ = (φ1 , . . . , φm )T ∈ Rmn , bi,− =
i
inf x∈X φ·x, bi,+ = supx∈X φ·x, ε = min(mini εi , mini Bζ i ),
X

ζ = mini ζ i and

Let ε := min(mini εi , mini

dmin = min (min dimin ,
I⊂[m] i∈I

min
i∈I,j∈[m]\I

dij
min )

X

φ · (x+ − x− ) ≥

ζi
BX i

). Then

ε||xi,+ − xi,− || ≥ ε||x+ − x− ||

i∈K

Assume that the following Collection Separability condition holds: for all i, j ≤ m, 6= j there exists dij
min > 0 s.t.
if x ∈ X is in the reachable set of H and xi ∈ Gie ∧ xj ∈
/
0
j
Gje0 ∀e0 ∈ E j then d(xj , Gje0 )) > dij
min for all e ∈ E where
E j is the edge set of Σj and Gje0 is a guard of Σj on edge
e0 ∈ E j . Then Σ is STORMED.
Proof. (S) In H, let y = (y 1 , . . . , y m ) = Re ((x1 , . . . , xm ))
and assume that it was H1 that caused the jump. Thus
y j = xj , j > 1. Write e = (`, `0 ). By Collection Separability,
j
j
d(y j , Gjej ) > d1j
min for all j > 1, e ∈ E , and by separability
1
1
1
1
1
of H1 d(y , Ge1 ) > dmin for all e ∈ E . So by d(y, G`0 ,`00 ) >
min(d1min , minj>1 d1j
min ) > dmin for any guard leading out of
`0 , and we have separability. The argument can be repeated
for any subset I ⊂ [m] of systems jumping simultaneously.
(T): The H flow θ` (t; x) is TISG because the component
flows θ`ii (t; xi ) are TISG.
(O) The cartesian product of definable sets is definable, so
the system H is o-minimal.
(RM) First we show that resets of H are monotonic, then
that the flows of H are monotonic. Let p, q ∈ L be two
modes of H, p 6= q.
Case 1: H jumps p → p. So any subsystem Hi either jumped
pi → pi or didn’t jump at all. If x+ = x− , then (RM)
is satisfied. Else, define φ := (φ1 , . . . , φm ) ∈ Rn·m , where
φi is the φ vector of system Hi . Then φ · (x+ − x− ) =
P
i
i,+
− xi,− ), where K ⊂ [m] is the set of indices
i∈K φ · (x
of sub-systems that jumped with xi,− 6= xi,+ . Note that K
depends
on x− , x+ . For all x− , x+ pairs (and so for all K)
P
i
i
i∈K ζ ≥ mini∈[m] ζ := ζ > 0. So by (RM) for each Hi ,
X

φ · (x+ − x− ) =

φi · (xi,+ − xi,− ) ≥

i∈K

X

ζi ≥ ζ > 0

i∈K

Thus (RM) is satisfied.
Case 2: H jumps p → q. At least one syb-system Hi jumped
P
pi → q i 6= pi . Then φ·(x+ −x− ) = i∈[m] φi ·(xi,+ −xi,− ) =
P
i
i,+
− xi,− ), where K = K= ∪ K6= ⊂ [m] and
i∈K φ · (x
K= is the index set of subsystems that jumped pi → pi
with xi,+ 6= xi,− , and K6= is the index set of subsystems
that jumped pi → q i 6= pi with xi,+ 6= xi,− . Subsystems that didn’t jump or jumped without changing their
continuous state don’t contribute to the sum. Note that
K= , K6= depend on x− , x+
have φ · (x+ − x− ) ≥
P
P. So we
i
i,+
i,−
i
− x || + i∈K= ζ .
i∈K6= ε ||x
For all X i , ||xi,+ − xi,− || ≤ BX i for all xi,− , xi,+ ∈ X i .
i,+
i,−
Therefore ζ i ||x B−xi || ≤ ζ i for all i ∈ K. So
X

+

−

φ · (x − x ) ≥
X
X
( min εi )||xi,+ − xi,+ || +
i∈K6=

X
i∈K6=

i∈[m]

( min εi )||xi,+ − xi,− || +
i∈[m]

i∈K=

X
i∈K=

ζi
||xi,+ − xi,− || ≥
BX i

So H has monotonic resets.
The flows of H are also monotonicP
along φ. Indeed for any
i
i
i
q ∈ L, φ · (θq (t + τ ; x) − θq (t; x)) = m
i=1 φ · (θq i (t + τ ; x ) −
P i
i
i
i
i
i
i
θqi (t; x )) ≥ i ε ||(θqi (t + τ ; x ) − θqi (t; x ))|| ≥ ε||(θq (t +
τ ; x) − θq (t; x))||
(ED) By Prop. 2.1.

7.

In general it is not possible to compute the reach sets required in Alg. 1 exactly unless the underlying o-minimal
theory is decidable. The HICD ||HCA closed loop is definable in Lexp , and the latter is not known to be decidable.
The authors in [21] proposed approximating the flows and
resets by polynomial flows and resets in the decidable theory LR . However, the approximation process is typically
iterative and requires manual intervention, or is restricted
to subclasses of STORMED systems [21].
Here we show that if an approximate reachability tool with
definable over-approximations is available for the continuous
dynamics, it can be used in Algo 1 (instead of exact reachability) to yield a finite simulation (rather than a bisimulation). Intuitively, the additional intersections of approximate reach sets with blocks of Q/ ∼ do not destroy finiteness of the procedure. Since we only have a simulation,
counter-examples on the abstraction should be validated in
a CEGAR-like fashion.
Lemma 7.1. Let Σ = (H, . . .) be a SHS and ∼ and equivalence relation on X. For any mode ` of H, its dynamical
sub-system D with state space X = H.X and flow θ` admits
a finite simulation S` that respects ∼, returned by Alg. 1.
The proof is in the Appendix. Let Ft (P) := ∩` S`∈L where
P = X/ ∼. Ftε refines all the S` ’s, and it is a finite simulaτ
tion of H by itself w.r.t. the continuous transition −
→. It is




clear that Ft (·) is idempotent: Ft (Ft (P)) = Ft (P)
Theorem 7.1. Let H be a STORMED hybrid system, and
P be a finite definable partition of its state space. Define
W0 = Ft (P),

i∈[m]

∀i ≥ 0, Wi+1 = Ft (Fd (Wi ))

Then there exists U ∈ N s.t. WU +1 = WU and
simulation of H by itself.

Ft (WU )

(7)
is a

Proof. By Lemma 10 of [27] there exists a uniform bound
U on the number of discrete transitions of any execution of
the STORMED system H, so Fd (Wk ) = Wk for all k ≥ U .
Moreover WU +1 = Ft (Fd (WU )) = Ft (WU ) and WU +2 =
Ft (Fd (WU +1 )) = Ft (WU +1 ) = Ft (Ft (WU )) = Ft (WU ) =
WU +1 , so the iterations reach a fixed point. The fact that
Ft (WU ) is a simulation then yields the desired result.

7.1

( min

FINITE SIMULATION FOR STORMED
SYSTEMS

Example: SpaceEx reachable sets

ζi
)||xi,+ − xi,− || Lemma 7.1 required that the over-approximation sets Rt ({x})
BX i
be definable for every x and t (see proof). In practice,

we need to show that the over-approximation actually computed by the reachability tool (which may not be the full
ball Rt (x)) is definable. In this section we show that the
over-approximations computed by SpaceEx [8] are definable. Given the set X ⊂ Rn and finite V ⊂ Rn , parameter
λ ∈ [0, 1] a time step δ > 0, and (i, j) ∈ E, SpaceEx overapproximates Rij (X) by K(V, X) := Rij (T HV (X) ∩ Gij ) ∩
Inv(j) and Rλδ (X) by [8]:
Ωλ (X, δ)

=

(1 − λ)X ⊕ eδA X

+
−
⊕ (λEΩ
(X, δ) ∩ (1 − λ)EΩ
(X, δ))

(8)

where T HV (X) := {x ∈ Rn | ∧~a∈V ~a · x ≤ ρ(~a, X)} is
+
the template hull of X and ρ its support function, EΩ
=
−
2
2 δA
(Φ2 (A X), EΩ = (Φ2 (A e X)), ⊕ is the Minkowski
sum, S = [−|x1 |, |x1 |] × . . . × [−|xn |, |xn |] is the box hull
with |xi | := max{|xi | s.t. x = (x1 , . . . , xn ) ∈ S}.
Theorem 7.2. For all definable polytopes X ⊂ Rn , the sets
K(V, X) and Ωλ (X, δ) is definable are Lexp .
Proof. Let S, Y ⊂ Rn be two definable sets in some ominimal structure A. Let λ ∈ R and let A be a real matrix.
Then the following sets are also o-minimal: λS, AS, S ∩ Y ,
S ⊕ Y , S ∩ Y , T HV (S) and S. Now the result follows
by noting that K(V, X) and Ωλ (X, δ) are constructed by
composing the above definability-preserving operations.

8.

CONCLUSION

In this paper, we presented the first formalization of a hybrid system model of the human heart and ICD closed loop
and showed that it admits a finite bisimulation, and that
definable approximate reachability yields a finite simulation
for STORMED systems.

9.

REFERENCES

[1] R. Alur, T. A. Henzinger, G. Lafferriere, and G. J.
Pappas. Discrete abstractions of hybrid systems.
Proceedings of the IEEE, 88(2), 2000.
[2] E. Bartocci, F. Corradini, M. D. Berardini,
E. Entcheva, S. Smolka, and R. Grosu. Modeling and
simulation of cardiac tissue using hybrid I/O
automata. Th. Com. Sci., 410(33), 2009.
[3] Boston Scientific Corporation. The Compass Technical Guide to Boston Scientific Cardiac Rhythm
Management Products. Device Documentation, 2007.
[4] T. Brihaye and C. Michaux. On the expressiveness
and decidability of o-minimal hybrid systems. Journal
of Complexity, 21(4):447 – 478, 2005.
[5] F. Cameron, G. Fainekos, D. Maahs, and
S. Sankaranarayanan. Towards a verified artificial
pancreas: Challenges and solutions for runtime
verification. In E. Bartocci and R. Majumdar, editors,
Runtime Verification, volume 9333 of Lecture Notes in
Computer Science, pages 3–17. Springer International
Publishing, 2015.
[6] T. Chen, M. Diciolla, M. Kwiatkowska, and
A. Mereacre. Quantitative verification of implantable
cardiac pacemakers over hybrid heart models.
Information and Computation, 236:87 – 101, 2014.
[7] D. D. Correa de Sa, N. Thompson,
J. Stinnett-Donnelly, P. Znojkiewicz, N. Habel, J. G.
Muller, J. H. Bates, J. S. Buzas, and P. S. Spector.
Electrogram fractionation. Circ Arrhythm

Electrophysiol, 55:909 – 916, Dec 2011.
[8] G. Frehse, C. L. Guernic, A. Donze, S. Cotton,
R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang,
and O. Maler. Spaceex: Scalable verification of hybrid
systems. In Proceedings of the 23d CAV, 2011.
[9] M. R. Gold et al. Prospective comparison of
discrimination algorithms to prevent inappropriate
ICD therapy: Primary results of the Rhythm ID
Going Head to Head Trial . Heart Rhythm, 9(3):370 –
377, 2012.
[10] R. Grosu, S. A. Smolka, F. Corradini, A. Wasilewska,
E. Entcheva, and E. Bartocci. Learning and detecting
emergent behavior in networks of cardiac myocytes.
Commun. ACM, 52(3):97–105, Mar. 2009.
[11] R. Hood. The EP Lab. Accessed 10/20/2015.
[12] Z. Huang, C. Fan, A. Mereacre, S. Mitra, and
M. Kwiatkowska. Invariant verification of nonlinear
hybrid automata networks of cardiac cells. In A. Biere
and R. Bloem, editors, CAV. 2014.
[13] M. A. Islam, A. Murthy, A. Girard, S. A. Smolka, and
R. Grosu. Compositionality results for cardiac cell
dynamics. HSCC, 2014.
[14] Z. Jiang, M. Pajic, S. Moarref, R. Alur, and
R. Mangharam. Modeling and Verification of a Dual
Chamber Implantable Pacemaker. Tools and
Algorithms for the Construction and Analysis of
Systems, 7214:188–203, 2012.
[15] R. Klabunde. Cardiovascular electrophysiology
concepts. Lippincott-Williams, 2 edition, 2011.
[16] S. Kong, S. Gao, W. Chen, and E. Clarke. dreach:
delta-reachability analysis for hybrid systems. In
C. Baier and C. Tinelli, editors, TACAS, volume 9035
of Lecture Notes in Computer Science. 2015.
[17] G. Lafferriere, G. J. Pappas, and S. Sastry. O-minimal
hybrid systems. Mathematics of Control, Signals and
Systems, 13(1):1–21, 2000.
[18] D. Mery and N. K. Singh. Pacemaker’s Functional
Behaviors in Event-B. Research report, INRIA, 2009.
[19] A. J. Moss et al. Reduction in inappropriate therapy
and mortality through icd programming. New England
Journal of Medicine, 367(24):2275–2283, 2012.
[20] M. Pajic, Z. Jiang, I. Lee, O. Sokolsky, and
R. Mangharam. Safety-critical medical device
development using the upp2sf model translation tool.
ACM Trans. Embed. Comput. Syst., 13(4), 2014.
[21] P. Prabhakar, V. Vladimerou, M. Viswanathan, and
G. E. Dullerud. Verifying tolerant systems using
polynomial approximations. In RTSS, 2009.
[22] M. Rosenqvist, T. Beyer, M. Block, K. Dulk,
J. Minten, and F. Lindemans. Adverse Events with
Transvenous Implantable Cardioverter-Defibrillators:
A Prospective Multi-center Study. Circulation, 1998.
[23] P. S. Spector. Visible EP. Accessed 10/20/2015.
[24] P. S. Spector, N. Habel, B. E. Sobel, and J. H. Bates.
Emergence of complex behavior: An interactive model
of cardiac excitation provides a powerful tool for
understanding electric propagation. Circulation:
Arrhythmia and Electrophysiology, 4(4):586–591, 2011.
[25] P. Tabuada. Verification and Control of Hybrid
Systems . Springer, 2008.
[26] K. Ten Tusscher, R. Hren, and A. V. Panfilov.
Organization of ventricular fibrillation in the human
heart. Circulation Research, 100(12):87–101, 2007.

[27] V. Vladimerou, P. Prabhakar, M. Viswanathan, and
G. Dullerud. Stormed hybrid systems. In Automata,
Languages and Programming. 2008.

APPENDIX
Proof of Lemma 5.3.
Proof. We show the resets are monotonic - the other properties are immediate. The state is x = (t, L2 , L1 , κ, σ2 )T .
The self-transition ACCUMULATE → ACCUMULATE is
initiated by VEvent (ventricular peak). At reset time, 0 ≤
t ≤ DL, we have that φ · (0 − t, t2 , t, 1, 0)T ≥ −φ1 DL +
W ant

φ4

≥ ζ.

The transition ACCUMULATE → FINALIZE, initiated at
the end of Duration, saves the value of the variance in σ2 .
This reset produces the constraint φ5 ((L2 − L21 /κ)/κ) ≥
ε|((L2 − L21 /κ)/κ)|. But the quantity in absolute value is
itself a variance and so is positive, therefore the constraint is
simply φ5 ≥ ε, compatible with the previous inequality.

Proof of Lemma 7.1.
Proof. This follows the lines of the elegant proof of [4] as
formulated in [25] and generalizes it to set-valued maps.
(The fact that using an approximate P ost operator yields
a simulation is a special case of a more general result on
transition systems but we prove it here for completeness.
Also note that this result holds for o-minimal systems [17]
generally, not just STORMED systems).
First observe that using approximate reachability on a system H is tantamount to replacing H with a system Hε whose
flows and reset maps are set-valued ε over-approximations
of the flows and resets of H (but is otherwise unchanged).
Therefore define the dynamical system Dε with state space
n
X and whose flow Θ : R × Rn → 2R is a set-valued ε overapproximation of θ` : Θ(t; x) = {y ∈ Rn | ||y − θ(t; x)||2 ≤
2 }. Let P := X/ ∼ be the partition induced by ∼. It follows
from the definability of θ and ||·||2 that Θ is definable. Given
P ∈ P, let Z(P ) = Θ−1 (P ) := {(x, t) | Θ(x, t) ∩ P 6= ∅}.
Then Z(P ) is definable because P and Θ are definable. Let
Zx (P ) = {t | (x, t) ∈ Z(P )} ⊂ R be the fiber of Z over x.
The number of connected components of Zx (P ) equals the
number of times that Θ(x, t) intersects P . Now it follows
from [25] Thm.7.11 that there exists a uniform upper bound
on the number of connected components of Zx (P ), independent of x. Let that bound be VP . Thus Θ(x, t) visits P at
the most VP times, regardless of x. Since there is a finite
number of blocks P ∈ P, then Θ(x, t) visits any block P a
maximum of V := maxP (VP ) times.
Thus we can associate to each x ∈ X a finite number of finite
strings q(x) = (`1 , `2 , . . . , `i−1 , `bi , `i+1 , . . . , `s ), where `i , `bi ∈
P. Each q(x) gives the sequence of blocks that Θ(x, t) visits
(with repetition), and in which `bi is the block containing
x. There may be more than one such string because the
set Θ(x, t) might intersect more than one block of P at a
time. The length of q(x) is thus uniformly upper-bounded by
V ·|P|, so there’s a finite number of different strings q(x). Let
Q(x) be the set of such strings associated to x, and let Q =
∪x Q(x). Then Q is the state space of the finite transition
system K = (Q, {∗}, −
→, Q0 ) whose transition relation is

∗
• `1 . . . `bi . . . `s −
→ `1 . . . `bi+1 . . . `s
∗
→ `1 . . . `s−1 `bs
• `1 . . . `s−1 `bs −

It is clear that K is non-deterministic and simulates D but
is not a bisimulation because of the over-approximation produced by Θ.

