A Cellular Automata Model of the Heart and Its Coupling with a Qualitative Model

COMPUTERS AND BIOMEDICAL RESEARCH ARTICLE NO.

29, 222–246 (1996)

0017

A Cellular Automata Model of the Heart and Its Coupling with a Qualitative Model P. SIREGAR, J. P. SINTEFF, M. CHAHINE,

AND

P. LEBEUX

De´partement d’Information Me´dicale, Universite´ de Rennes I, 35033 Rennes Cedex, France

Received March 15, 1995

Cellular Automata (CA) models offer a good compromise between computational complexity and biological plausibility while qualitative models have expressive power for explicitly describing dynamic processes. In this paper we present a 2D CA model and its coupling with a qualitative model. The CA model includes elements characterizing muscle, nodal tissue, and bypass conduction. Each element exhibits adaptive properties to cycle length and to the prematurity of incoming impulses. A crude electrocardiogram is also simulated via an equivalent source formulation. Arrhythmias such as the Wenckebach phenomenon, atrial flutter, or extrasystole-triggered tachyarrhythmias can be simulated using relatively simple models when they incorporate the fast conduction system with muscle tissue and when the model elements exhibit adaptive properties. We then illustrate how a CA model can be coupled to a qualitative model to produce a system that combines the fine grained description of CA models with the high level interpretative role of qualitative models.  1996 Academic Press, Inc.

I. INTRODUCTION In a previous paper (27) we described a qualitative model of the heart and discussed its role as a simulator and as a multilevel descriptive tool to study cardiac rhythmic disorders. The model is part of the CARDIOLAB framework which, ultimately, will constitute an intelligent computational environment for teaching, research, diagnosis, and monitoring (25–28). We discussed the central role of qualification in different tasks such as qualitative simulation, explanation, and decision-making. In model-based diagnosis, it is essential that computational models of a target system should be quick in reproducing the range of behavior they are intended to simulate since the problem is to explain some observed or measured data. Qualitative models are particularly well suited to achieve this task (6, 23, 25). However, their generally crude level of detail can be a limiting factor for simulating certain kinds of rhythmic disorders such as fibrillation. More detailed models are necessary and the CA approach has thus far shown promising results. It has been adopted to model the specialized conduction system (9, 18), patches of muscle (4, 5, 21), as well as the ventricles (3, 17, 20). On the other 222 0010-4809/96 $18.00 Copyright  1996 by Academic Press, Inc. All rights of reproduction in any form reserved.

CELLULAR AUTOMATA MODEL OF THE HEART

223

hand, CA models cannot—as stand alone models—perform diagnosis, nor are they capable of describing in any intelligible way what they are computing. Hence our effort has been directed toward creating a multiscale computational model of cardiac electrical activity with both quantitative and qualitative features. In this paper, we describe a 25,000-elements 2D CA model and show how it can be coupled with a qualitative model that has been previously reported (27). The resulting system combines the fine-grained description of CA models with the high-level interpretative role of the qualitative model. II. CURRENT RESEARCH Four problems can be associated to modeling and simulation in ECG: (1) model the electrical conduction system, (2) model the electrical activity of muscle tissues, (3) model the effects of volume conduction, and (4) produce summaries (descriptions) of the salient events or states associated to a particular ECG simulation. To the authors’ knowledge, no single model tackles these four problems that, together, address both quantification and qualification. Indeed, in general, heart models adopt one of three representational schemes: (i) qualitative representations based on formalisms such as first-order logic (FOL) where processes are qualitatively described (6, 23, 27), (ii) CA models that mimic electrophysiologic processes as discrete time, space, and state simulations (3–5, 9, 17, 18, 20, 21), and (iii) networks of cells governed by the so-called Hodgkin–Huxley (H–H) type equations (14). The H–H type equations compute the electrodynamics of each cell in terms of nonlinear relationships between transmembrane currents and transmembrane potential (13, 15, 16, 29). While the first kind of model addresses the qualification problem, the second and third kinds address (broadly speaking for CA models) quantification. Models that simulate quantitative data (e.g, isochrones, depolarization wavefronts, propagating impulses) can be separated into two broad classes: models of the specialized conduction system (9, 18), and models of cardiac muscle tissue (3–5, 13, 15–18, 20, 21, 29). A complete model should explain arrhythmias in terms of conduction impairment and impulse formation disturbances of the conduction system and the myocardium. As mentioned in (9), the most versatile models should include realistic three-dimensional (3D) geometry in the specialized conduction system and the myocardium, while simulating an ECG through an equivalent source formulation (3, 20). For model-based diagnosis, we believe that to be truly versatile, a framework should also include different grain-sized models, and most importantly, it should also possess a metalevel allowing it to describe and ‘‘know’’ what it is computing. That is, the framework should produce concept symbols denoting high level descriptions from which decisions could be inferred, possible course of actions undertaken, and simulation success/failure summaries automatically produced. Moreover, the price to pay if a single fine-grained quantitative model is used, is the computational space and time complexities. Indeed, in modelbased diagnosis, the problem of model parametrization is a totally understated one. The adjustments required to tune the model-parameters so that the ensuing

224

SIREGAR ET AL.

simulation matches the data is simply intractable if complex models are used. For a given context, some a priori information about initial tissue states should either be known in advance or be derived by using more efficient means than brute force adjustments of the model. Thus, model-based diagnosis of arrhythmias may call for the use of interactive models of intermediate complexity. These models would, at the same time, offer some of the desired features of the most versatile models. In this paper we describe such a model. It is a 2D, CA-qualitative model of the specialized conduction system and muscle tissue. An ECG is simulated via an equivalent-source formulation and serves as a visual aid for analyzing the time-course of the model’s behavior. III. THE CELLULAR AUTOMATA MODEL CA models mimic physical laws by a series of simple rules that are easy to compute. In CA heart models, each element may represent hundreds or thousands of individual cardiac cells. For the myocardium, state space may be reduced to two states: active and refractory or resting and excitable (4, 5). The timevarying patterns obtained by these models are maps exhibiting the states of each model element at sampled time instants (3–5, 17, 20, 21). Electrical conductionsystem models are generally built as networks partitioned into segments (9, 18). Segment parameters may include the slow diastolic depolarization, the antero and retrograde conductions, and the refractory periods (12). Our model includes both nodal and muscle tissue. It is made up of approximately 25,000 square elements. Assuming characteristic size from base to apex neighboring 12 to 15 cm, each element corresponds to tissue patches whose surface is around 0.8 mm2. In the present model, nodal and tissue elements take one of four states (Fig. 1a): resting, depolarization, absolute refractory, and relative refractory, while in the specialized conduction tissue elements, the resting state is substituted with the slow-diastolic-depolarization. The elements form a two-dimensional array, and at discrete instants, a state matrix S(i, j ) defines the state of the (i, j )th element. Each element undergoes a natural state transition (e.g., absolute refractory to relative refractory), or a forced state-transition (e.g., relative refractory to depolarization). Adaptive properties have also been introduced since many arrhythmias cannot be reproduced if rate dependencies and slowed conduction—as in the case of a propagating impulse in partially recovered tissue—are ignored. They are given as follows: Frequency adaptation. This property is modeled by a function particularized for each model cell. The refractory period of a given cell adapts to the length of its preceding cycle (i.e., to its current rate of depolarization) (Figs. 1b and 2). When the cycle length increases, the refractory period increases and inversely. The rate dependency at a given instant, n, is given by a function adjusted from real data (2, 7, 10, 24) and has the following form:

tabs (n) 5 tabs,857 1 (Dtabs(Tn21) 2 Dtabs(857)),

[1]

CELLULAR AUTOMATA MODEL OF THE HEART

225

FIG. 1. (a) The states and state transitions of the model’s cells. Four states characterize muscle and nodal tissue cells. For muscle tissue cells, the transition from resting to depolarizing (or from relative refractory to depolarizing) is triggered when neighbor cells are in an active state. For nodal tissue cells, depolarization automatically follows slow diastolic depolarization. The transition from relative refractory state to depolarization (for muscle tissue) or from slow diastolic depolarization to depolarization (for nodal tissue) occurs when the give cell is discharged by an incoming excitation wavefront. (b) Theoretical adaptation curve of the absolute refractory period to the cycle length. Dtabs designates the increasing duration of the absolute refractory period with cycle length. This function is used to determine the instantaneous value of tabs as defined in Eq. [1]. This curve is based on data given in (2). (c) Theoretical adaptation curve to impulse prematurity. When a cell is discharged by an excitation wavefront while in a relative refractory state, the cell’s depolarization 1 duration increases by a factor given byt dep . When the prematurity P (see Eq. [2]) equals one, the impulse is blocked and no conduction occurs. When P equals zero, the impulse is normally conducted across the cell. This function is based on data given in (10).

226

SIREGAR ET AL.

FIG. 2. Adaptive properties of the model elements. Upon user request, the transmembrane potential, fm , of any given element can be displayed during the simulation. This figure depicts the time-course of fm of a right ventricle element at a normal sinus pace of 70 bpms (upper display) and during the ventricular fibrillation shown in plate 7. Observe the narrower action potential corresponding to the shorter cycle-length during fibrillation when compared to the normal sinus rhythm. Also, the influence of premature depolarizations on the conduction speed can be assessed in the lower display. The slope of the action potential decreases as the prematurity of the depolarization increases. This results in a slowed conduction across the given element.

where Tn21 is the preceding cycle’s length, tabs,857 is the duration of the refractory period for a normal 70 beats per minute (bpms) rate (857 ms cycle length), and Dtabs(.) is the function depicted in Fig. 1a. Conduction velocity adaptation. In accordance to known physiological data (7), conduction delays are modeled by adapting each element’s depolarization duration to the prematurity of the afferent impulses. An element in a partialrefractory phase will be depolarized by an incoming impulse but the rise-time will be increased by a factor proportional to the partial refractoriness. Conduction delays across a given cell are quantified according to the relation depicted in Fig. 1c: 1 t dep 5 f (tdep , P)

P512

Dt , trel

[2]

where P is the prematurity factor, Dt is the time lag between the beginning of the refractory-state and the instant the impulse reaches the structure, trel is the cell’s relative refractory period, and tdep is the normal depolarization duration (i.e., corresponding to a frequency of 70 bpms). To the authors’ knowledge, there is no exact expression for the function f. However, this function has been extrapolated from data obtained in the guinea pig (9). To summarize, an element’s absolute refractory period is reduced when the cycle length is reduced, and conduction velocity in a given element is also reduced

CELLULAR AUTOMATA MODEL OF THE HEART

227

if an afferent impulse reactivates the element in its relative refractory period. The reference cycle length is 70 bpms and the adjustments are based on data collected from the literature (2, 7, 9, 10, 24). Propagation of the excitation wavefront. Propagation in cardiac tissue corresponds to the action potential that satisfies the wave equation. The case where the velocity c is constant and isotropic with conductivity s is given by:

s=2fm 5

1 ­ 2fm , c 2 ­t 2

[3]

where fm is the potential across the cell membrane. CA models mimic this equation using heuristics such as the Huygens construction method for isotropic and anisotropic propagation (4, 5, 11, 21). At each discrete instant, propagation is modeled by generating circular or ellipsoidal propagation wavelets around excited cells (Fig. 3a). The size of the wavelet corresponds to the distance traveled by the wavelet between consecutive instants. Excitable cells inside the wavelet are depolarized and the global activation wavefront is obtained by repeating the process over the entire model at each time increment. Depending on the ratio of the local conduction vectors, local isotropic or anisotropic propagation can be achieved (Fig. 3b). IV. MODELING

THE

ECG

The relationship between the propagating wavefronts of activation and repolarization on the one hand, and the ECG on the other, can be analyzed by simulating the spatio-temporal potential gradients in the myocardium and the effect of volume conduction. This is the so-called forward problem (3, 8, 11, 17, 20). Volume and body surface potentials at a given point r can be predicted from the spatio-temporal source distribution Imv (r, t) that in turn follows from Poisson’s equation =.s=fe 5 2=.Js 2=.Js 5 Imv ,

[4]

where fe is the extra-cellular potential and Js( A/cm2) is the impressed or source current density. Using the bidomain theory (11, 20), Js can be obtained from the cells’ transmembrane potential Js 5 s=fm .

[5]

In the isotropic case, s is a scalar and can be expressed as a bulk conductivity that depends on the interstitial (so) and intracellular (s i) conductivities (17):

s5

si so . si 1 so

[6]

228

SIREGAR ET AL.

FIG. 3. (a) Huygen’s construction method. Propagation is modeled by generating circular or ellipsoidal propagation wavelets around excited cells. Let C 5 (cx , cy ) designate the conduction velocity across the two-dimensional medium, and let Dt designate the time interval following the onset of a wavelet initiation. Then the wavelet axes (rx , ry ) are given by rx 5 Cx Dt and ry 5 Cy Dt, for isotropic propagation. Excitable cells inside the wavelet are depolarized and the global activation wavefront is obtained by repeating the process over the entire model at each time increment. (b) Anisotropic propagation in the model. In this figure, the anisotropic propagation results in the difference of activation wavefront velocity across nodal and muscle tissues.

In the anisotropic case, the ratio sTan / sOrth can be assigned for each element, where sTan designates the conductivity along the direction tangent to the muscle fibers, and sOrth designates the orthogonal conductivity. As the model stands now, fm is not calculated in terms of transmembrane ion flow, but is given a priori as a finite support continuous piecewise linear function that mimics the action potential. Given an element e, the i th phase is associated to a linear function of the form:

CELLULAR AUTOMATA MODEL OF THE HEART

229

Ye, i (t) ; ae, i (t). t 1 be, i (t),

[7]

where the parameters ae, i (.) and be, i (.) are time-varying since they are context dependent, as depicted in Fig. 2. Hence for an element, e, fm is given by:

fm, e (t) 5 Ye, i (t) if at the time sample t, the element e is in state i.

[8]

For now, the context of a given cell e is reduced to three factors. Namely, its state, the eventual prematurity of an incoming impulse, and the preceding cycle’s length (the time between two consecutive depolarizations). As noted, the cell e will be depolarized if it is in a partially or totally recovered state, and its rate of depolarization will be dependent on whether it has partially recovered or not. As for the refractory periods, they are conditioned by the cell’s preceding cycle length. Now, substituting Eqs. [6] and [7] into [5] yields an equivalent dipole moment density that can be substituted into [4] to compute the current source Imv . The extracellular potential, fe , at a given point r, and hence the simulated ECG, was computed by solving Eq. [4] for an infinite homogeneous medium. Since the CA model is planar, one cannot expect to have a correct simulated ECG by using the computed Imv as is. We ‘‘solved’’ this problem by introducing a third dimension as follows: Since during the depolarization of the ventricles, the radial component of the impressed current densities plays a most important role, we eliminated the current sinks (negative Imv’s) as if they were behind the current sources with respect to the observation point r. The result yielded a qualitatively acceptable ECG as depicted in plates I to VII. Of course, the gross approximation excludes any quantitative analysis of the simulated signal. That is, the signal is essentially a visual aid from which basic information such as atrial and ventricular rates can be derived. V. COUPLING

THE

CA MODEL

WITH THE

QUALITATIVE MODEL

In model-based diagnosis, there are good reasons why a CA model should be coupled to a qualitative model. Among them we can include: (i) qualitative models can help determine the initial conditions of more detailed models (Fig. 4), and (ii) they can translate CA results into high-level concept symbols whose semantics are straightforward (Fig. 5). Qualitative models operate into the symbolic realm by explicitly representing high level concepts (6, 23, 27). In qualitative process simulation predictions of global behavior can be inferred by rules such as ‘‘if the AVN is in an absolute refractory state then it cannot be excited by an incoming impulse.’’ In CARDIOLAB’s qualitative model, both the specialized conduction system and the myocardium are represented (27). There are 21 elements belonging to muscle or nodal tissue and the level of detail corresponds to the one usually adopted by clinicians when they explain arrhythmias. Each structure exhibits the same phases as those characterizing the elements of the CA model. Autorhythmicity of the nodal tissue elements and frequency adaptation are part of the region’s properties. All representations are based on firstorder logic (FOL) and have been implemented in PROLOG. Dynamic events

230

SIREGAR ET AL.

FIG. 4. Instantiating the CA model using the qualitative model. This task is particularly important for explaining observed data such as in model-based diagnosis. Given a particular real system behavior, the goal is to explain the behavior in terms of initial model-parameter values. A grossly determined initial parameter vector can be established using the qualitative model. Depending on how well the simulated and desired behavior match, the qualitative model’s parameters can be further tuned or the CA model’s initial parameter vector computed by the Instantiate task. The CA parameters can then be fine-tuned in a hypothesized-and-test loop.

have finite temporal support and are represented within the framework of the so-called FOL temporal logics (1, 19). In Allen’s temporal logic (1), events are represented as terms of the ‘‘Holds’’ predicate which explicitly states the temporal occurrence of the event. For instance, the fact that ‘‘during the interval [t0 , t1[, the upper atria (UA), was in an absolute refractory state’’ can be written as Holds(Absolute-refractory-state(UA), During([t0 , t1[)). In CARDIOLAB, there are three levels of abstraction (Fig. 5): the PROCESS, EVENT, and DIAGNOSTIC levels. Events of the PROCESS level constitute the canonical set of events

CELLULAR AUTOMATA MODEL OF THE HEART

231

FIG. 5. Interpreting CA simulation results using the qualitative model. (a) Qualitative process simulation. Given user-defined initial conditions, the qualitative model generates canonical descriptions in the PROCESS level workspace. These descriptions are abstracted into EVENT and DIAGNOSTIC descriptions (see text and (27)). (b) It is the CA model that generates the PROCESS level concept symbols (see also Fig. 6). The qualitative model’s knowledge sources are then triggered to derive higher-level concept symbols in the EVENT and DIAGNOSTIC workspace.

from which all other events are derived. They are generated in the system’s blackboard via qualitative process simulation that corresponds to the deduction of events when the initial conditions (i.e., the initial theory) are given. PROCESS level facts basically describe the phase states and transitions (27). Conduction sequences are abstracted from the PROCESS level and written at the next EVENT level using FOL intentional definitions (27). Basically, they are temporal assertions about propagating wavefronts such as ‘‘a branch to branch reentry occurred between t0 and t1’’. Finally, EVENT level descriptions are further abstracted at the DIAGNOSTIC level where the heart is considered as a whole. It corresponds to the level of description usually employed in clinical diagnosis. Typical concepts include those of impulse formation, e.g., premature beat, capture beat, escape beat, and those of impulse conduction, e.g., normally and abnormally

232

SIREGAR ET AL.

conducted impulses. For instance, the system generates the following: ‘‘a normal sinus rhythm at 70 bpms occurred between t0 and t1.’’ The coupling between the CA model and the qualitative model is established by substituting qualitative process simulation by the CA model simulation (Figs. 4 and 6). In other words, PROCESS level events are produced in the blackboard by the CA model (Fig. 6). This is achieved by assigning semantic links between PROCESS level events and CA-derived relations. The general strategy is to assign to each elementary event of the PROCESS level, one or more discriminating relations concerning quantities simulated by the CA model. Such quantities include the spatio-temporal patterns produced by the tissue phases, scalar fields such as the transmembrane potential fm , or vector fields such as the current sources Js . Given a region and given a temporal interval, we have the following metarule: ;p;I add Holds(p, I ) in the blackboard if and only if Rp( f ) is true during I. The variable symbol p designates an elementary event, Rp( f ) designates the corresponding relation in the CA realm, and f is some function of the CAsimulated quantity. Because of the bijection, the elementary events generated in the blackboard are determined in function of which of the Rp’s hold (see below). For instance, suppose that the transmembrane potential fm (.) is chosen as the discriminating quantity: Let Eh.j be the time average over a given time interval [t0 , t1[, and we let r denote a region of muscle tissue, then the following equivalencies can be used: Qualitative model Holds(depolarization(r), [t0 , t1[)

CA model } 0 , ldep # E

Holds(absolute-refractory(r), [t0 , t1[) } lrel , E

HE

r

HE ­­tf

m

Holds(relative-refractory(r), [t0 , t1[)

} E

Holds(resting-state(r), [t0 , t1[)

} «#E

r

HE

HE

r

­ fm dr ­t

J

J

­ fm dr , « , 0 ­t

J

[9]

dr # lrel

J

­ fm dr , ldep r ­t

FIG. 6. Qualification of the CA phase patterns into PROCESS level concepts. While the CA model simulates data such as the spatio-temporal phase patterns, the spatio-temporal patterns of scalar fields such as the transmembrane potential fm , or the spatio-temporal patterns of vector fields such as the current sources Js , events of the PROCESS level are concurrently generated by the model. They are subsequently abstracted in the qualitative model’s blackboard by triggering the appropriate knowledge sources (see 27).

CELLULAR AUTOMATA MODEL OF THE HEART

233

234

SIREGAR ET AL. TABLE 1 REFERENCE VALUES

OF THE

MODEL’S ‘‘NORMAL’’ HEART

AT

70

BPMS

Regions

Conduction speed (mm/s)

Absolute refractory (ms)

Relative refractory (ms)

Atria AV node His Bundle Bundle branches Ventricles

800 150 1500 2000 400

230 230 120 120 190

100 100 120 120 150

Here, f 5 E her ­ fm / ­t drj and Rp ( f )’s are the inequalities derived from the CA model. Since the logical foundation of the current model is standard logic, the tissuestates of any given region over any given interval are mutually exclusive. Therefore, the right-hand relations concerning the scalar f should also be mutually exclusive. In fact, this is true for any relation concerning any given CA-simulated quantity: the logical status between the temporal logic assertions and the Rp relations should be the same. Another aspect that needs a thorough investigation is the computing cost which depends, to a great extent, on the quantity being used as the discriminating factor. This problem will be discussed in a following paper. VI. RESULTS Seven simulation sessions are illustrated in the plates below. There are no limits to the number of rhythmic disorders that can be produced since before each simulation session the user chooses interactively the initial conditions. An important point is that in many instances the CA and qualitative models agree. That is, the phase durations of the CA model element comprising one of the 21 main regions of the qualitative model can be searched first via the physiological window of the qualitative model (27) until some desired global behavior is reached and then fine-tuned via the physiological window of the CA model (Fig. 4). Within each of the main regions, subregions of arbitrary shapes and sizes can be delimited with the corresponding tissue parameter values specified independently from the surrounding tissue. By this means, differential conductivity and refractoriness can be introduced for any given patch of tissue. Since each of the model’s cells is adaptive, any differential may increase, decrease, or remain the same throughout a simulation session. In the model, the reference state is a normal heart at a rate of 70 bpms. Approximate parameter values associated to this state have been derived from the literature (2, 7, 22). Table 1 depicts the values adopted for the model. In the plates below, all modifications reflecting pathological predisposing factors of arrhythmias are given as percentage ratios of the reference state values. Finally, context dependent variations are computed

CELLULAR AUTOMATA MODEL OF THE HEART

235

TABLE 2 SIMULATION TIME

Phase fm Imn ECG

OF A ‘‘NORMAL’’ HEART CYCLE (70 IBM/PC DX2/66 COMPUTER

BPMS) ON

CPU time(s)

CPU 1 DISPLAY time(s)

45 160 420 480

160 270 660 510

Note. The user can interactively choose any combination of the above information. Notably, the phase patterns can be selected alone. This operates a significant reduction in CPU time if one considers the time needed to simulate an ECG (i.e., compute the phases, fm , J, Imn , and then fe ).

with respect to this reference state. This property introduces versatility into the model while keeping computation-time reasonable. Table 2 gives some figures concerning this aspect. Note that the most important information with respect to cardiac electrical activity, i.e., uses about 10 times less CPU-time than the simulation of the ECG. One of the reasons is that the extracellular potential fields are directly computed from the current-source distributions. Important time-reduction can be obtained by merging the sources in predefined regions and then computing the fields from the local equivalent source (20). In each plate, selected instants illustrate the state (or phase) of each model element at the displayed time sample. The extracellular potential selected at a point 5–6 cm above the heart yields the simulated ECG. One of the artifacts of the crude ECG approximation is the absence of negative P waves that normally can be observed when the atrium is depolarized retrogradely. Plate 1: Atrial flutter. This disorder has been obtained by simulating a macroreentry at the level of the patch in the left atria (T 5 130 ms, T 5 190 ms), leading to a stable spatio-temporal pattern during the first 4000 ms. From around 4000 ms on, a more disorganized self-sustained activity sets in. The ventricular rate is constant at around 60 bpms. The ECG shows pattern reminiscent of the Fwaves that characterize atrial flutter at a frequency of around 350 bpms. The predisposing factors are a twofold decrease in conduction velocity (again, with respect to the reference state) and a 50% decrease of the refractory periods in the above-mentioned patch. The trigger is the late repolarization of this area while a depolarization front propagates through the surrounding medium (e.g., T 5 130 ms). Plate 2: Atrial fibrillation and first-degree AVblock. The predisposing factor is 15 patches in the atria exhibiting reduced velocity and a diminished duration of the effective refractory periods. The mean velocity over the 15 patches reflects a fourfold decrease and the range of dispersion is around 300 mm/s. Likewise, the mean refractory period over these patches is 110 ms which corresponds to about 40% decrease with respect to the reference value, and the range of disper-

SIREGAR ET AL.

PLATE 1

236

237

PLATE 2

CELLULAR AUTOMATA MODEL OF THE HEART

238

SIREGAR ET AL.

sion is 50 ms. Hence in this case study, both inhomogeneous conduction and dispersion of the refractory periods characterize the atrial tissue. The trigger is two simultaneous left and right atrial extra-systoles at a period when the tissue is only partially repolarized. Multiple independent reentries (e.g., T 5 155–510 ms) result from these events. The interfering initially self-sustained excitation loci are followed by higher frequency depolarization wavefronts as can be seen in the icons corresponding to T 5 2440–5860 ms. From the second ventricular beat on, the ventricular rate is regular at around 60 bpms, reflecting impaired conduction in the AV node. Plate 3: A 4 : 3 Wenckebach phenomenon. The predisposing factor is a twofold increase of both absolute and relative refractory periods of the AVN. The sinus rate is constant at around 78 bpms. The three first sinus impulses are conducted to the ventricles, with a slowing of the AV node conduction for the second and third impulses as can be observed by the increasing PR intervals (e.g., T 5 60–1870 ms). The fourth sinus impulse is blocked at the AVN (T 5 2375–2560 ms), and a new cycle is initiated. The simulation of decremental conduction in the AV node is made possible by the adaptive response with respect to the prematurity of incoming impulses (Figs. 1 and 2). Plate 4: Trigeminy. The arrhythmia was obtained as a result of two predisposing conditions: (i) there is a differential in refractoriness between the left (twofold increase) and right (50% decrease) AV node pathways, and (ii) the upper bundle branch’s absolute refractory periods are abnormally long (400% increase). The differential in refractoriness leads to an echo originating in the AV node, while the bundle branch block causes an antegrade depolarization of the left ventricle followed by a retrograde activation of the right ventricle. The following cyclic conduction is observed: (1) The first sinus impulse is normally conducted in the ventricles (T 5 60–195 ms), (2) the second sinus impulse is conducted in the ventricles with a delay in the AV node, (3) while the ventricle are depolarizing, a retrograde reentry, at the level of the right AV node, is occurring (T 5 760 ms), leading to a retrograde depolarization of the atria (T 5 905–1010 ms), (4) the traveling impulse follows a path in a clockwise motion in the AV node and reactivates antegradely the ventricles in the presence of a phasic right bundle branch block (T 5 1195–1260 ms) (AV node zoom). The process repeats itself. Plate 5: Paroxysmal tachycardia, first degree block, and phasic right bundle branch block. This simulation was obtained by a 50% reduction of the conduction velocity in the AV node. Otherwise, the initial conditions are identical to the preceding trigeminy case. The main difference is a decrease of conduction speed when the media is in a partially recovered state. The third ventricular beat corresponds to an echo originating in the AV node, similar to the one illustrated in plate 4. However, the second echo (sixth ventricular beat) initiates a tachycardia that results in a local lower right bundle branch self-sustained reentry (T 5 3175–3845 ms) that reactivates the ventricles. The tachycardia lasted for around 6 sec. Plate 6: Ventricular flutter. The predisposing factor is a 30% reduction of conduction velocity in the ventricles and a 50% reduction of the absolute refractory period in the ventricles. This condition characterizes ischemic tissue. The

239

PLATE 3

CELLULAR AUTOMATA MODEL OF THE HEART

SIREGAR ET AL.

PLATE 4

240

241

PLATE 5

CELLULAR AUTOMATA MODEL OF THE HEART

SIREGAR ET AL.

PLATE 6

242

CELLULAR AUTOMATA MODEL OF THE HEART

243

arrhythmia is triggered by two reentries in the right and left ventricles over patches of partially recovered tissue (e.g., T 5 160–265 ms). This results in a branch to branch activation between the right and left bundle branch (e.g., T 5 1380–1440 ms). The spatio-temporal patterns stabilize at around 4400 ms where the ventricular rate is stable and approximates 240 bpms. The ventricles control the atria as depicted by the retrograde atrial depolarization (e.g., T 5 780– 800 ms). Plate 7: Ventricular fibrillation. The predisposing factor is 10 patches in the ventricles exhibiting a dispersion of repolarization (range 5 120 ms), a 50% reduction and dispersed (range 5 90 mm/s) in the ventricles. The trigger is a lower-right bundle branch extra-systole (T 5 5 ms) initiating a reentrant circuit (two in the right ventricle and one in the left (T 5 275–800 ms)). The circuits interfere with each other but remain in existence throughout most of the simulation. Two or three vortices could be observed with spirals moving clockwise and sometimes counter-clockwise. The resulting condition is a desynchronization of the ventricles as can be seen in the simulated ECG. VII. DISCUSSION

AND

CONCLUSION

In this paper, we have described an interactive 2D CA model of cardiac electrical activity and gave an overview on how it can be coupled with a qualitative model previously reported (27). Both CA and qualitative models are designed to simulate propagation in the heart’s specialized conduction system and in the myocardium. They offer two complementary approaches to modeling: (i) they reflect two distinct levels of detail, and (ii) they separately produce quantitative data and qualitative descriptions. The model resulting from their coupling offers both features: The advantages of qualitative models with their ease of manipulation, computer-time efficiency and expressivity, and the usually more appropriate level of detail of finer-grained models. The main role of the CA model presented here is to reproduce arbitrarily complex and qualitatively plausible arrhythmias, while keeping computer-time within reasonable bounds. For diagnosis, quick assessment of global behavior can be obtained by tuning the qualitative-CA models’ initial parameter values until some desired behavior is obtained (Fig. 4). These initial conditions can then be used as a first approximation to a more correct quantitative analysis using 3D models in a multistage resolution process. Much work remains to be done in this area, however, and answers to the following questions should improve the efficiency of the coupled Qualitative-CA system: (1) What are the converging behavioral modes between the qualitative model and the CA model? (2) where do they diverge? and (3) what do we do about the discrepancies? Despite the obvious difficulties lying ahead, there are good reasons for integrating quantification and qualification in frameworks dedicated to the study and diagnosis of complex systems: While theoretical quantitative studies provide the soundest approach for explaining experimental data, qualification could be central for the development of intelligent systems via which detailed models could be analyzed, interpreted, and where natural language

SIREGAR ET AL.

PLATE 7

244

CELLULAR AUTOMATA MODEL OF THE HEART

245

summaries highlighting interesting points could be conveyed to the user or stored in a knowledge base. In classical mathematical models, all forms of reasoning about the processes are implicit in the equations. Describing dynamic properties of a model by interpreting the simulation results is done at a meta-level by a human agent and usually in a nonmathematical language where the equations and the results constitute the object level (26). One possible route then is to embody detailed computational models within a reasoning process to produce an intelligent computational framework that would play an active role in problemsolving tasks such as explaining real data, correlating events in a data base, explaining the salient points of a simulation, or simply asking the appropriate questions. To reach this goal, the qualitative and quantitative realms need some meeting point where strings of numbers and 2D or 3D numerical field patterns can be translated into semantically clear concept symbols and vice versa. The material presented in this paper could be one starting point since it describes a multiscale approach to modeling and an introductory work on the problem of translating simulated data into qualitative descriptions at different levels of abstraction. ACKNOWLEDGMENT This work was fully supported by the Conseil Re´gional de Bretagne.

REFERENCES 1. ALLEN, J. F., Towards a general theory of action and time, Artif. Intell. 23, 123 (1984). 2. ANTZELEVITCH, C., LITOVSKY, S. H., AND LUKAS, A. Epicardium versus endocardium: Electrophysiology and pharmacology. In ‘‘Cardiac Electrophysiology—From Cell to Bedside,’’ (D. P. Zipes and J. Jalife, Eds.), pp. 386–395. Saunders, Philadelphia, 1990. 3. AOKI, M., OKAMOTO, Y., MUSHA, T., AND HARUMI, K. Three-dimensional simulation of the ventricular depolarization and repolarization processes and body surface potentials: Normal heart and bundle branch block. IEEE Trans. Biomed. Eng. 34, 454 (1987). 4. AUGER, P., BARDOU, A., COULOMBE, A., GOVAERE, M. C., ROCHETTE, L., SCHREIBER, J. P., VON EUW, D., AND CHESNAIS, M. Computer simulation of ventricular fibrillation and of defibrillating electric shocks. Effects of Antiarrhythmic drugs. Math. Comput. Modelling 14, 576–581 (1990). 5. BARDOU, A., AUGER, P., DUMEE, P., FARAH, R., BOURNIER, P., AND GOVAERE, M. C. Cellular automata models applied to cardiac arrhythmias. Modelling, Measurement & Control, AMSE Press, 38, (No 4), 23–31 (1993). 6. BRATKO, I., MOZETIC, I., AND LAVRAC, N. Automatic Synthesis and Compression of Cardiological Knowledge. Machine Intelligence 11 (1986). 7. BRAUNWALD, E.‘‘ Heart Diseases, A Textbook of Cardiovascular Medicine,’’ Saunders, Philadelphia, 1984. 8. CUFFIN, B. N., AND GEZELOWITZ, D. B. Studies of the electrocardiogram using realistic cardiac and torso models. IEEE Trans. Biomed. Eng. 242, 24 (1977). 9. DASSEN, W., BRUGADA, P., RICHARDS, D., GREEN, M., HEDDLE, B., AND WELLEN, H. A mathematical model of the conduction system to study the mechanisms of cardiac arrhythmias. In ‘‘Computers in Cardiology’’ (K. L. Ripley, Ed.), IEEE Comp. Soc. Press, Siver Spring, MD, 1982. 10. DELMAR, M., AND JALIFE, J. Wenckebach Periodicity: From deductive electrocardiographic analysis to ionic mechanisms. In ‘‘Cardiac Electrophysiology—From Cell to Bedside’’ (D. P. Zipes and J. Jalife, Eds.), pp. 128–138, Saunders, Philadelphia, 1990.

246

SIREGAR ET AL.

11. GESELOWITZ, D. B., On the theory of the electrocardiogram. Proc. IEEE 77, 6 (1989). 12. GULRAJANI, R. M. Models of electrical activity of the heart and computer simulation of the electrocardiogram. CRC Crit. Rev. Biomed. Eng. 16 (Issue 1), 1(1988). 13. GUL’KO, F. B., AND PETROV, A. A. Mechanism of the formation of closed pathways of conduction in excitable media. Biophysics 17, 271 (1972). 14. HODGKIN, A. L., AND HUXLEY, A. F. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117 (1952). 15. JANSE, M. J., AND VAN CAPELLE, F. J. L. Electrotonic interactions across an inexcitable region as a cause of ectopic activity in acute regional myocardial ischemia. Circ. Res. 50, 527 (1982). 16. KEENER, J. P., A mathematical model for the vulnerable phase in the myocardium. Math. Biosci. 3–18 (1988). 17. LORANGE, M., AND GULRAJANI, R. M., Computer simulation of the Wolf-Parkinson-White preexcitation syndrome with a modified Miller-Gezelowitz heart model. IEEE Trans. Biomed. Eng. 33, 862 (1986). 18. MALIK, M., COCHRANE, T., AND CAMM, A. J. Computer simulation of the cardiac conduction system. Comput. Biomed. Res. 16, 454 (1983). 19. MCDERMOTT, D. A temporal logic for reasoning about process and plans. Cog. Sci. 101 (1982). 20. MILLER, W. T., AND GESELOWITZ, D. B. Simulation studies of the electrocardiogram, I. The normal heart. Circ. Res. 43, 301 (1978). 21. MOE, G. K., RHEINBOLT, W. C., AND ALBISCOV, J. A. A computer model of atrial fibrillation. Am. Heart J. 67, 200 (1964). 22. SHAMROTH, L. ‘‘The Disorder or Cardiac Rhythm,’’ Blackwell Sci., Oxford, 1975. 23. SHIBAHARA, T. On using knowledge to recognize vital signals: Knowledge-based interpretation of arrhythmias. IJCAI 307 (1985). 24. SHRIER, A., DUBARSKI, H., ROSENGARTEN, M., GUEVARA, M. R., NATTEL, S., AND GLASS, L., Prediction of Complex Atrio-ventricular Conduction rhythms in humans with use of the Atrioventricular nodal recovery curve. Circulation 76 1196–1205 (1987). 25. SIREGAR, P., MABO P., AND COATRIEUX, J. L. How can deep knowledge be used in CCU Monitoring? IEEE Eng. Med. Biol. 4 (N0 ? 12), 92–99 (1993). 26. SIREGAR, P. Theoretical Cardiology: From Mathematical to Qualitative Models. J. Biol. Sys. 4 (N0 ? 1), 1996. 27. SIREGAR, P., CHAHINE, M., LEMOULEC, F., AND LE BEUX, P. An Interactive Qualitative Model in Cardiology Comput. Biomed. Res. 28 (N0 ? 6), 443 (1995). 28. SIREGAR, P., CHAHINE, M., SINTEFF, J. P. AND LEBEUX, P. Deep knowledge and computer-assisted instruction in cardiology. MEDINFO, Vancouver, Canada, 1995. 29. VAN CAPELLE, F. J. L., AND DURRER, D. Computer simulation of arrhythmias in a network of coupled excitable elements. Circ. Res. 47, 454 (1980).