- Email: [email protected]

AND

Computer

BIOMEDICAL

Simulation

MAREKMALIK,*

16,454-468

RESEARCH

(1983)

of the Cardiac

THOMAS

COCHRANE,~

Conduction

System

AND A. JOHN CAMMS

Computing Unit for Medicul Sciences and Department of CardioioLqys St. Bartholomew’s Hospital. 38 Little Britain, London ECIA 7BE, Great Britain Received March 4. 1983

A discrete computer model simulating the operation of the cardiac impulse transmission apparatus and the electropotential changes of the heart musculature has been developed on a NORD-100 minicomputer. The model is written in NORD-FORTRAN and allows description of practically all basic pathologies of the transmission apparatus. A simulated ECG curve is produced in each case. The current version of the model is especially suited for studies of heart rhythms. Pacemaker descriptions (of most important pacemaker types) have been introduced into the model and the computer computations make it possible to assess the effects of different pacemaker modes under different pathological situations. The principle on which the model is based is described and the paper shows concisely several examples in the form of simulated ECG records.

1. INTR~DuCTI~N In general, the methods of computer modeling and computer simulation can be used when investigating any complex system, the logic of which and the relationships between different parts of which are comparatively well known. If the modeled system is known exactly enough, its computer description can be exact too. Results obtained by experimenting with such an exact model can be easily applied to the original. On the other hand, in the case where the computer representation of the system under investigation contains some important simplifications, the results of the model are only approximate and their application to the original cannot be used for essential conclusions. This fact is probably the most important limitation to the use of computer modeling for investigation of biomedical systems. The relationships within biomedical systems are usually only known roughly and our knowledge about many biomedical systems is restricted to the description of their external behavior. The description of such external manifestations cannot be used to * Address correspondence and reprint requests to Dr. Marek Malik, Dept. of Computer Science, Charles University, Malostranske namesti 25, 118 00 Praha 1, Czechoslovakia. t From Dept. of Medical Electronics, St. Bartholomew’s Hospital, London. $ From Dept. of Cardiology, St. Bartholomew’s Hospital, London. 454 0010-4809/83 $3.00 Copyright 0 1983 by Academic Press, Inc. All tights of reproduction in any form reserved.

SIMULATION

OF THE

CARDIAC

CONDUCTION

SYSTEM

455

prepare the specification of a model. Thus the method of computer simulation can be successfully used to study only a few biomedical systems. The cardiac conduction system and the electropotential changes in heart muscle form one biomedical system the logical operation of which is known comparatively accurately. Many properties of this system and many relations between its different parts have been described in an exact mathematical way (1). Action-potential changes of muscle fibers have been modeled using mathematical description based on differential equations (2, 3). This paper describes a complex computer model of the cardiac conduction system which has been created using the methods of discrete event simulation (4). The paper deals with the principles on which the model is based, describes its current realization, and shows some examples of its results in the form of simulated ECG records. 2. MODEL

DESCRIPTION

The aim of the whole project has been the creation of a model which would enable simulation of the principal known abnormalities and changes of the cardiac conduction system. This means that the structure of the model had to consist of several parts corresponding to the more important parts (from the point of view of the cardiac electropotential activity) of the heart and its conduction system. The heart musculature, conducting fibers, both nodes, and some pathological structures (such as parasystolical centers, preexcitation fibers bypassing the AV node, etc.) had to be represented in the model both in data structures and in algorithms. We shall describe the principal data structures of the model and its algorithms separately for each part of the model. 2.1 Model Realization

of the Heart Muscle

The common data structures used in the realization of the model can be divided into two parts. The first group of data structures corresponds to the simulation calendar, i.e., contains variables serving for synchronization of all simulation processes. The other data group contains the parameters and local variables of the different simulation processes. For the computer model, the heart musculature has been considered to consist of a great number of hexagonal elements, which are ordered into small plane triangular-shaped structures. Each of these triangles contains 45 basic elements which are numbered (identically within each triangle) and the model contains a special table giving the numbers of neighbors to each element in the triangle (Fig. 1). These triangular structures have been used as construction elements to create the three-dimensional representation of the heart muscle (i.e., the three-dimensional image is formed from two-dimensional planes only). The triangular form of these elements enables them to be put together in different ways and model different three-dimensional structures (the triangular form has been chosen for its inherent simplicity).

456

MALIK,

COCHRANE,

23 ' 25

FIG. 1. Six-sided “musculature” The central area (elements 40.

AND CAMM

3

elements are arranged into triangular-shaped planar structures. . 45) may be contacted by a Purkinje fiber.

To assemble triangular structures, it is necessary to describe how they are connected. It is supposed that each triangle can connect with one or two others at each side. (Two different neighbors are used to describe how the septa contact both ventricles or both atria.) With regard to the numbering of triangle elements, there are 12 different ways to link one triangle to one side of another. The model contains another common table describing these possibilities. All triangles are numbered and the connection between two of them is described by the numbers of both linked structures and by a pointer into the table containing the linkage description, Each small musculature element is determined by the number of its triangle and by its position within it. This data organization enables the neighbors to each musculature element (some of them may not be in the same triangle) to be searched rapidly. Each musculature element represents independent instances of a simulation process which consists of three events. ’ These events describe the depolarization and repolarization of an element in the following simplified way: if an element receives an action signal, it waits for a certain time, then sends the signal to its neighbors and eventually to other structures (for instance, to conduction fibres-in those pathological situations when retrograde conduction occurs) and waits again for a longer time when signal transmission to this element is impossible (the repolarization period). The way in which this action influences the simulated ECG signal will be described later. The individual process instances of all musculature elements are planned to be activated at a certain value of the simulated system clock. The activations of all these processes have to be synchronized and this synchronization means (in effect) that the minimum of all plan values has to be searched at each step of the I The technical term simulation process is used for the algorithmical prototype describing a part of modeled actions. (For instance, the action of muscle elements is described by a fixed algorithm-the musculature process.) During the model calculation, several dynamic copies of this prototype may be created. These copies are called the process instances.

SIMULATION

OF THE CARDIAC

CONDUCTION

SYSTEM

457

FIG. 2. Diagram of the data structure containing the planes of musculature process instances. The instances of elements belonging to one triangle are ordered in “local” trees and these trees are ordered in a “global” tree. (The open circles show the leftmost elements of the local trees.)

computation. Since it can be supposed that the whole heart muscle will be modeled by a great number of such small elements, a more complicated data structure making it possible to search the minimal plan quickly has been used. The plan values of process instances belonging to all elements in one triangular block are ordered in a binary tree and the common data on each triangle contain two pointers to the root of the corresponding tree and to the leftmost element (containing the minimal plan value). The triangular blocks themselves are ordered in another global binary tree after the values of their local minimal plans (Fig. 2). These local minima have to be calculated only if the leftmost element of the local tree is changed; otherwise they are accessible using the pointers to the leftmost elements. In this hierarchical structure, the searching and the plan changes (canceling and new planning of events) can be done with low algorithmical complexity (5). The model algorithms relating to the representation of the musculature realize the simulation process (repeating the three events discussed) and handle the hierarchical tree structure. 2.2 Realization

of the Conduction

System

The transmission system has been supposed to consist of separate linear sections. Each of these sections has again been realized by a process instance. The sections can be connected together and a picture of the conduction tree can

458

MALIK,

COCHRANE,

AND CAMM

be built from them. In addition, the sections are oriented (the beginning and end of a section are distinguished) which means that the differences between anterograde and retrograde transmission can be implemented. As with the musculature blocks, the conduction sections are numbered and a common data structure saves the numbers of front and back neighbors for each section. The sections having no front neighbors represent the Purkinje fibers and each of these contacts the central area of one musculature triangle and can transmit the depolarization signal to the musculature elements in this area. The simulation process realization is different for the transmission fibers in atria and in ventricles since it is supposed that the ventricular fibres can radiate their own impulses under special circumstances (corresponding to their own spontaneous repetition frequencies). The state diagram of the process realizing one ventricular section is shown in Fig. 3. The signal transmission consists of three separate events: the section waits for a certain time and sends the signal in the anterograde direction, then waits again for a certain time and sends the signal in the retrograde direction, and finally waits for the repolarization period during which time no signal can be elaborated (i.e., transmitted to this section and sent in the anterograde and/or retrograde way). The organization of the process realizing operation of atria1 sections is simpler because no spontaneous activity is supposed. The time plans of process instances of all transmission sections are saved in two ordered binary trees, one of which contains the plan values of atria1 sections, the other plan values of ventricular sections. The algorithmical realization of the conduction system is composed partly of routines computing the process actions and partly of routines handling these calendar trees.

FIG. 3. State diagram of the simulation process realizing one ventricular conduction section. SE,. SE2, and SE3 (= Signal Elaboration 1, 2, 3) are the three separate events of signal transmission described in the text. WS (= Wait Start) is the initial state in which the process waits before starting to radiate its own impulses. WP (= Wait Period) is the state in which the process waits for its own period. The “s” state changes correspond to the explicit changes when a signal has been received; the other state changes are the implicit changes made after time plans.

SIMULATION

2.3 Realization

OF THE

CARDIAC

CONDUCTION

SYSTEM

459

of the Nodes

Special data structures are used to introduce the sinus and atrioventricular node into the model. Both structures are unique and their simulation processes each exist in only one instance. Consequently, there are no problems with their synchronization. The process realization of the sinus node is quite simple. The node sends the depolarization signal periodically to the atrial transmission fibers and waits for its repolarization time. Transmission fibers conducting from the node are able to depolarize it retrogradely (for instance, if an atria1 parasystolic center is considered). On the contrary, the process realization of the AV node is the most complicated process in the model. Similar to the ventricular transmission sections, the node is able to radiate its own impulses and, in addition, it contains an inner counter which sums the signals obtained from the atria1 conduction system. The signal is sent to the ventricles only when the value of this counter exceeds a radiation threshold. The model enables either of the nodes to be protected in such a way that no signal can be transmitted to it. Such a node radiates its own impulses regularly and no pathological situation can depolarize it. Other special data structures are used to introduce pathological parasystolic centers into the model. Such a center is connected to the central area of a musculature triangle and periodically activates the musculature elements in this area. (Such activation is possible only if the elements in the connected area are not depolarized.) 2.4 Result Collection Since the model has to serve clinical cardiological research, it is desirable that the output be available in an easily interpreted form. For this we have chosen to simulate the ECG signal in three orthogonal XYZ projections. The orthogonal projections can be modeled in a somewhat simpler way than the full 1Zchannel record and for the initial experimenting with the model, the orthogonal projections have been fully sufficient. On the other hand, to extend the model and to simulate the full standard ECG record would be comparatively easy. (We shah revert to this point in the last section of the paper.) The signal is produced by the model computation in digital form, filtered and smoothed (using the method of moving averages), and finally plotted on a graphical terminal. To produce the raw digital form, the polarization changes from the whole heart muscle are simulated. The output contributions are taken only from the musculature elements. Depolarization and repolarization of each musculature element is described by a common sampled function (the same for all elements). The sampled values are read as input data allowing the function to be changed easily.

460

MALIK,

COCHRANE,

AND

CAMM

The way in which the depolarization and repolarization of each musculature element influences the changes of the polarization vector of the whole heart depends on the topological position of this element. It is further supposed that the topological contribution of all musculature elements belonging to one triangular block is the same. The data for each musculature triangle describe this topological position by three projection constants. (As we have mentioned, the modeled polarization vector is projected into three directions corresponding approximately to the orthogonal XYZ channels and each projection constant refers to one channel.) If a muscle element is depolarized, the common output signal is multiplied by the projection constants of the corresponding block and the products obtained for all muscle elements are summed. The global output signal F, for the channel w is (in the course of one systole) given by the formula

where M is the number of musculature blocks, pi,*, is the projection constant of the kth triangle for channel W, fis the elementary output signal, and qk,] is the time when the jth element of the kth triangle is depolarized. The projection constants correspond not only to the topological position of the triangle, but also to the thickness of the musculature area which the triangle represents. The difference between the ECG signal contribution of the atria and of left and right ventricles is realized by differences in projection constants for the different triangles. 2.5 Pacemakers

For special rhythm studies, the possibility to model different pacemakers was requested. In addition, the model had to allow different pacemaker modes to be introduced into the simulated system. The common data structures contain descriptions of atria1 and ventricular pacing areas in the form of sets of musculature triangles. The atria1 or ventricular pacing then corresponds to the activation of muscle elements belonging to one of the blocks in the area considered. The routines realizing atria1 or ventricular stimulation contribute to the output data. Pacing spikes are modeled in this way. Similarly, common data structures contain descriptions of atrial and ventricular sensing areas and two logical variables enable initiation of sensing and inhibition of sensing during the pacemaker refractory period. If the sensing key is open and a musculature element in the sensing area is depolarized, a signal is sent to the pacemaker, the sensing key is closed, and the simulation process containing the pacing algorithm is activated. For the pacemaker algorithms, the model is open in the sense that different routines realizing different pacemaker modes but using the prescribed common data structures can be linked together with other parts of the model.

SIMULATION

OF THE CARDIAC CONDUCTION

3. MODELREALIZATION

SYSTEM

461

AND EXPERIMENTS

The model has been developed on a NORD-100 minicomputer in the Computing Unit for Medical Sciences at St. Bartholomew’s Hospital, London. Currently, it is in its sixth version, The program is written in NORD-FORTRAN (an extension of the FORTRAN-IV standard) and has approximately 7000 lines. Result calculation requires approximately 120 kbytes. The calculation speed depends on the complexity of the situation which is simulated; to model 6 set of real time (the case for results presented here) needs between 5 (simple cases) and 15 (most complicated cases) min CPU time. 3.1 fnput Data

Atria1 musculature, including the atria1 septum, has been modeled by 54 and the ventricular musculature, including the septum, by 88 triangular blocks. The topological form of the heart muscle has been simplified and represented by an ellipsoid with a central plane representing the septum (Fig. 4). The common signal describing the action-potential changes of one musculature element has been taken in the form of a sampled approximation of the first derivative of the action-potential curve. It consists of two nonzero parts; the first corresponds to the depolarization change, the second to the repolarization phase. Modeled P waves and QRS complexes originate as the summation of depolarization signals; the T and Tp waves originate as the sum of repolarization outputs. To describe the conduction system, a detailed map of the fiber tree has been prepared. The map of atria1 part of this system is shown in Fig. 5 (linear conductive sections and contacted musculature blocks are numbered). A similar map has been prepared to describe the ventricular system. In atria 65, and in ventricles 181, linear fibers have been considered (Fig. 6).

FIG. 4. The three-dimensional image of the modeled heart muscle is built from planar triangles. The topological form is simplified and only an ellipsoid with a central plane (the septum) is considered.

462

MALIK,

COCHRANE,

AND CAMM

FIG. 5. Diagram of the atria1 conduction system. The numbers in circles represent the numbering of transmission fibers; the numbers in squares designate the conducted musculature triangles.

3.2 Experiments

The fact that the musculature is modeled by point elements is a quite imporPolarity changes of ECG signal when depolarization and repolarization radiate through the musculature in a direction opposite to normal systole cannot be modeled satisfactorily. This influences the results in cases where such a specific pathological situation is modeled. For some important morphological changes, however, the model gives correct results. Figure 7 shows the simulated ECG when a normal conduction cycle was supposed (case a) and the development of right bundle branch block (case b). For rhythm studies, the current version of the model seems to be adequate. It permits descriptions of individual rhythm pathologies and simulation of the situation when all of them are present simultaneously. Figure 8 contains the simulation result for sinus tachycardia with second degree AV block (2 : I). In this case, the AV node threshold for the impulse radiation to the ventricles has been increased in such a way that a QRS complex occurs after every second P wave. Figure 9 shows the simulated ECG of the combination of normal sinus rhythm with a parasystolic center in the right ventricle. Figure IO shows a modeled paroxysm of AV reentrant tachycardia caused by connecting a musculature block in the superior part of left ventricle to another one in the inferior part of left atrium. This connection has been done by a simple change in the tant simplification.

SIMULATION

OF THE CARDIAC CONDUCTION

SYSTEM

463

FIG. 6. The global organisation of the conduction system. Black dots are the connections between neighboring fibres; the open circles correspond to the contacts between transmission apparatus and musculature.

input data describing the topological organization of the musculature triangles and the contacts between them. The current version of the model contains six special routines modeling basic pacemaker types VOO, DOO, VVI, DVI, VDD, and DDD (6) and it is possible to describe further routines. Figure 11 shows the simulated effect of a DVI pacemaker in the case of underlying sinus bradycardia with protected sinus node. The pacemaker paces both atria and ventricles periodically with a constant AV delay. A normal QRS complex which occurs during this AV delay inhibits ventricular pacing (first complex in the figure) and each normal QRS complex starts a new pacing period (this can be seen after the second normal complex in the figure). Results

464

MALIK,

COCHRANE,

AND CAMM

FIG. 7. Simulated ECG of (a) normal conduction cycle, and (b) the development of right bundle branch block.

CZI

V"

Y"

1"

I FIG.

1

f I

8. Simulated sinus tachycardia with second degree AV block.

SIMULATION

OF THE CARDIAC CONDUCTION

AA

CZI

’ ”

-”

”

V”

V”

465

SYSTEM

AA

AA

1”

Y”

FIG. 9. A simulated combination of normal sinus rhythm with a pathological center in the right ventricle which generates its own impulses.

presented in Fig. 12 compare the effects of VDD (case a) and DDD (case b) pacemakers in the same situation of sinus bradycardia with a complete AV block and with sinus and AV nodes both protected. The VDD pacemaker senses but does not stimulate the atria: some complexes have no P wave. Finally Fig. 13 shows a simulation of a DDD pacemaker causing endless-loop tachycardia in presence of AV block and intact retrograde VA conduction. 4. PERSPECTIVES

AND DISCUSSION

The experiments with the current version of the model give several suggestions how the model might be developed in future versions.

I

I

FIG. 10. A model of AV reentrant tachycardia caused by connecting together two muscular elements from left ventricle and left atrium.

FIG. I I. Simulated effect of a DVI pacemaker and underlying sinus bradycardia with protected sinus node (i.e., entrance block).

b FIG. 12. Simulated VDD (a) and DDD (b) pacemaker in sinus bradycardia, complete AV block. Both nodes are protected. 466

SIMULATION

OF THE CARDIAC CONDUCTION

SYSTEM

467

FIG. 13. Simulation of a DDD pacemaker causing tachycardia in the case of anterograde AV block and intact retrograde VA conduction. The pacemaker stimulates the ventricles; the depolarization is radiated from ventricles to atria where it is sensed by the pacemaker. In response to the sensing of atrial depolarization, the pacemaker paces the ventricles and this sequence is repeated resulting in so-called endless-loop AV tachycardia.

The most important limitation of the model has been mentioned. In the next version the pointlike description of the musculature will be changed and a fibershaped structure will be introduced for the muscle elements. For detailed morphological studies, it seems to be necessary to introduce a truly three-dimensional image of the musculature and not only a three-dimensional shell built from two-dimensional planar blocks. Another limitation is the presupposition that the elementary output signal is the same for all musculature elements. This prevents consideration of different properties of muscular elements in different areas (such as, for instance, result from regional ischemia, etc.). A more realistic approach would be not to realize the output contribution by a common discrete signal but to describe the polarization changes of each muscle element by a differential equation, realize the algorithms solving these equations in the form of simulation processes, and synchronize them in the way in which the depolarization and repolarization is radiated through heart muscle. In realizing this idea, the computation time would be increased several times. It is, therefore, thought expedient to have two different versions of the model. One of these could be simpler and useful in initial stages of cardiological research projects (and, perhaps, for educational purposes); the other would be more complex and more realistic, but its computation would be slower. For the more realistic version, parametric descriptions of ECG channels could be introduced in a quite simple way. This would enable the collected channels to be different for different cases. It would be possible to model special channels (such as an intracardiac or an esophageal channel). In addition, the next versions of the model must consider the cycle rate dependence of signal transmission velocity and duration of repolarization. This

468

MALIK,

COCHRANE,

AND CAMM

dependence is very important especially for pacemaker studies such as the simulation of tachycardia terminating pacemakers, etc. The new version of the model will be developed in the near future and it should prove helpful for evaluation of new and special pacemaker types and in other cardiac research projects. ACKNOWLEDGMENTS The authors thank Mr. A. W. Dunlop and Mr. A. Sakal from the Computing Unit for Medical Sciences at St. Bartholomew’s Hospital, London, for their helpful suggestions and discussions. During the course of this work, Dr. Malik was supported by the British Council. Dr. Cochrane acknowledges the help of the Joint Research Board of St. Bartholomew’s Hospital. Dr. Camm is a Wellcome senior lecturer.

REFERENCES C. V., AND GESELOWITZ, D. B. (Eds.), “The Theoretical Basis of Electrocardiology.” Oxford Univ. Press, Oxford, 1976. DROUHARD. J. P., AND ROBERGE, F. A., The simulation of repohuization events of the cardiac Purkinje fiber action potential, IEEE Trans. Bio. Med. Electron. BME29(7), 481 (1982). DROUHARD, J. P., AND ROBERGE, F. A., A simulation study of the ventricular myocardial action potential, IEEE Trans. Bio. Med. Electron. BME29(7), 494 (1982). ZEIGLER, B. P., “Theory of Modelling and Simulation,” Wiley, New York, 1976. Vol. 3, Addison-Wesley, New York, KNUTH, D. E.. “The Art of Computer Programming,”

1. NELSON, 2. 3. 4. 5.

1973. 6. PARSONNET,

V., FURMAN. S., AND SMYTH, N. P. D., “A Revised Code for Pacemaker Identification, PACE,” Vol. 4, pp. 400-403, July-August 1981.