Keywords: Simulation, Cellular Automata, Hodgkin-Huxley, Luo-Rudy, ECG, Arrhythmia 

Introduction

Cellular automaton models of the electrical activity, in contrast, are computationally very efficient. The rules they are based on, however, usually are cartoon-like representations of the properties of excitable media or, at most, are applicable only in a very specific situation.

In this paper, the hybrid method is presented: This approach tries to combine the precision and flexibility of the membrane models with the computational efficiency of CAs. For this purpose a (quasi)one-dimensional tissue of cells, based on the membrane model of Luo and Rudy [2], is studied under a variety of situations. The results of these simulations are translated into situation-dependent rules for the CA. In this way, the CA receives a physiologically justified foundation and allows - within the limited range of situations - the realistic simulation of the electrical activity on larger scales; even, as shown here, the electrical behaviour of the whole human heart - with its activation and repolarisation patterns, surface potentials, ECGs - is accessible to simulation. In addition, the effects of drugs that modify the dynamics of membrane channels can easily be implemented in this scheme. The hybrid method allows to study the variation of the ECG under the influence of e.g. class I antiarrhythmia.

The Hybrid Method

Figure 1: Scheme of the hybrid method

The Linear Chain

For the membrane model, the Luo and Rudy model of the electrical activity of single ventricular cells [2] is chosen. Beside the fast inward sodium current I_Na, whose formulation was in the main focus of this model, a secondary inward, a slow calcium current and four outward potassium currents are represented. The equations of the different ionic currents are of the Hodgkin-Huxley [3] type.
To include drug effects, e.g. the effects of class I antiarrhythmia, the formulation of the fast sodium inward current has to be modified. Here, the "guarded receptor model" by Starmer et al. [4] is incorporated. With this, an additional equation for the relative number b of blocked sodium channels has to be solved: 

db dt = k D · (1-h) · (1-b)-L · b

(1)

D denotes the concentration of the drug, k the association rate of the drug to the sodium channel and L the analogous dissociation rate. h describes the inactivation gate of the sodium channel. The sodium current in the presence of an antiarrhythmic drug

IBNa = (1 - b) · INa

(2)

shows a use-dependent behaviour: The relative number of blocked channels b increases during the course of an action potential and decays exponentially thereafter. A higher frequency of excitation - in other words a lower basic cycle length BCL - leads to an increased number of blocked sodium channels.

To investigate the behaviour of collective phenomena, a small tissue sample has to be built with the Rudy and Luo cells. According to physiology the cell-to-cell coupling is ohmic. The simplest choice is a (quasi)one-dimensional medium: the linear chain. To avoid boundary effects, a special protocol is used for the computer experiments:

Figure 2: Coupling procedure used by the computer experiments with the linear chain.

We split the configuration in three chains. In the beginning chain A and B are coupled (symbolized by K_ab). Cell #1 is stimulated N times with a fixed basic cycle rate BCL=1000 ms. (N has to be large enough to get "stationary conditions" for this value of BCL.) This leads to N waves of excitation in the direction of cell #270. The Nth stimulation of cell #1 is associated with t=0. At a certain time t=t^* cell #271 is stimulated leading to a wave of excitation in the direction of #330. If the transmembrane potential of #300 is greater than that of cell #31, then chain A and B are uncoupled and C and B are linked together (K_cb), so that - instead of #301 - cell #31 is adjacent to #300.
The results refer to the N+1st stimulation of the central part (cells #60 to #240) of chain B. The behaviour of these cells strongly depends on the individual last cycle length t_hist, that is the actual individual time difference to the last excitation. The cell's response to the excitation is in general an action potential (AP) which propagates with a certain velocity v. The shape of the AP may well be characterized by the following parameters: _peak, the peak value of the AP, t_peak the duration of the fast upstroke (phase 0) of the AP, _plateau, the value of the transmembrane potential _M during the plateau phase and t_ref, the duration of the AP.

For a BCL=1000 ms the parameters concerning the shape of the AP and the velocity v are displayed as functions of t_hist in figs. 3 and 4 respectively. In these figures, the results for three different situations are plotted: Squares label the simulation with normal model parameters, that means the original Luo and Rudy model coupled to a one-dimensional medium. Triangles mark the case where the maximum sodium conductivity is reduced by 66%; this is approximately the sodium conductivity, which is effective in the case of the class I antiarrhythmic drug (BCL=1000 ms, kD=4.0 Hz, L=0.5 Hz) at t_hist=1000 ms. The simulations with this class I antiarrhythmic drug are symbolized by circles.

Figure 3: The influence of the cycle length t_hist on the characteristic parameters of the AP. Squares label the simulation with normal parameters. Circles represent the results with the class I drug (kD=4.0 Hz, L=0.5 Hz) and triangles the case with the reduced sodium conductivity (by 66%).

Figure 4: The influence of the cycle length t_hist on the (normalized) propagation velocity of the excitation wave. Squares label the simulation with normal conditions. Circles represent the velocities with the class I drug and triangles the case with the reduced sodium conductivity.

As expected the curves for the class I drug and of the reduced sodium conductivity cross at about t_hist=1000 ms. All of the class-I-drug-associated curves show the use dependent behaviour, the feature of recovery for longer cycle lengths. A simple reduction in sodium conductance is inappropriate to mimic the effects of class I drugs.

The Cellular Automaton (CA)

In this section, a short description of the components of the CA (geometry, neighborhood, states, rules) will be given.
A data set of the geometry of a real human heart was made available by Dr. Meinzer from the German Cancer Research Center (DKFZ). According to these data, approximately 90000 model-cells are arranged on a cubic array with a grid spacing g=0.15 cm. To account for the parameter variations of the heart cells (typical parameters), the different regions of the heart (atria, ventricles) and specific structures (sinoatrial node, atrioventricular node, specific conducting system, isolating layer between atria and ventricles) had to be identified manually. The condition for the neighborhood is a distance smaller or equal to .
In the CA, the cells take one of the following three states: During the quiescent state (I), the transmembrane potential is equal to the resting membrane potential, the active state (II) marks the starting point of the AP and the refractory state (III) covers the whole course of the AP: upstroke, plateau phase and the complete recovery.

A graphical representation of the rules of the CA is given in figure 5.

Figure 5: Flow diagram for the state evolution of a single cell in the cellular automaton. Assignments are symbolized by rectangles, ovals mean queries. If one condition is met, the vertical line is the path to follow. See text.

A more detailed description of the rules can be found elsewhere [1].
quiescent state (I): If a cell is in the quiescent state, the value of the counter history, this corresponds to the cycle length, will be increased by one. An excitation, a change from the quiescent state to the active state, is only possible if a certain time since the last AP has elapsed and the membrane potential of the neighboring cells exceeds a certain threshold. The value of Err_max determines the number of timesteps the cell has to stay in state I; this controls the velocity of excitation.

active state (II): The cell changes in the refractory state.

refractory state (III): This state is divided in two parts, in the absolute refractory state and the relative refractory state, governed by the number of timesteps (Refrak) elapsed in this state (relative to Refrak_max). Only in the second time period a cell's state could be changed to the active one. For this, conditions similar to transition I to II have to be met.

The Determination of the Actual Parameters

For every model cell, the number of timesteps since the last excitation (history in figure 7) determines the t_hist-value and using the simulations with the linear chain the corresponding characteristic parameters. These parameters are modified depending on the location of the heart (as mentioned above).
These actual parameters are then used to calculate the membrane potential depending on the value of Refrak, which is given by the CA.

So, the CA uses the results of the linear chain to compute for every timestep and cell the value of the transmembrane potential. With this and certain assumptions (quasistationarity, bidomain model, infinite torso model), the extracardial potential can be related to the Laplacian of the transmembrane potential (see review of Gulrajani [5]). The positions of the electrodes on the "surface" of the infinite torso are chosen according to Thiry and Rosenberg [6].
As all cells of the model heart are assumed to behave like ventricular cells, there is no automaticity built in in this model: A "heartbeat" has to be initiated "manually" by changing the state of the cells representing the sinoatrial node to the active state.

Figure 6: Simulation of the standard-12 ECG with normal electric parameters.

A comparison of the duration of the simulated P-Wave with typical values yields a time resolution of =0.15 ms.
The simulated standard-12 ECG (e.g. right axis deviation) resembles ECGs of enlarged right chambers. These were the reasons to check again the geometric properties of the heart model: It seems, that the geometric data were derived from a right hypertrophic human heart.

Another point of interest is the simulated sequence of excitation and repolarisation. In figures 7 and 8, these electrical events are shown color coded on slices of the model heart. The upper left (lower right) slice corresponds to the most (least) distant portion of the heart (frontal perspective).

Figure 7: Isochrones of activation for normal parameters from a frontal perspective. Lead V₃ is displayed at the bottom.

Figure 8: Isochrones of repolarisation for normal parameters from a frontal perspective. Lead V₅ is displayed at the bottom.

The Isochrones of activation are defined by the time (relative to SAN), when the cells are in the active state. Similarly, the isochrones of repolarisation are defined by the transition refractory to quiescent state. The sequence of excitation (figure 7) complies well with the description given by Durrer et al. [7]. In addition the repolarisation patterns are in satisfactory agreement with the qualitative behaviour given in the literature [8,9]: The ventricular cells are activated first, which are located on the endocardial layers of the heart exhibit the longest refractory periods and repolarise at the very end of the repolarisation sequence and vice versa.
Atrial cells have more or less a constant refractory period; therefore, their isochrones for activation and repolarisation look quite similar.

Now, the influence of the cycle length (t_hist) as well as drug effects on the electrical activity of the whole heart model, expressed by the ECG is demonstrated. To get the influence of the cycle length on the different components (waves, intervals, segments) these structures have to be identified. This is done by visual inspection of the lead V₆ and its time derivative.
As examples, the PQ segment, the QRS complex and the QT interval are selected: The duration of each of these structures as a function of the cycle length t_hist is displayed. This is done for three cases (analogous to the simulations with the linear chain): The results with normal model parameters are labeled with squares. Triangles mark the case where the maximum sodium conductivity is reduced by 66% (for motivation see text above). The simulations with the class I antiarrhythmic drug are symbolized by circles.

Figure 9: Duration of the PQ segment t_PQ as a function of the cycle length t_hist. For the explanation of the symbols, see figure 3 or 4.

Figure 10: Duration of the QRS complex t_QRS as a function of the cycle length t_hist. For the explanation of the symbols, see figure 3 or 4.

Figure 11: Duration of the QT interval t_QT as a function of the cycle length t_hist. The parallel lines mark the curve t_QT^{normal, max} = (t_hist)^0.5 · 0.42 s (with t_hist in seconds), the maximal value of normal QT intervals for fixed cycle lengths, that means steady state conditions. For the explanation of the symbols, see figure 3 or 4.

In presence of the antiarrhythmic drug, the PQ interval is enlarged by 51%; this is mainly the effect of the increased PQ segment (79%). The very slight enlargement of the QT interval (+4%) stems from the increased QRS duration (+17%). The JT interval keeps its value constant.
These results with the hypothetical class I drug, at least for steady state conditions (cycle length t_hist=1000 ms), are consistent with "clinical properties" of class Ic drugs: Class Ic medicaments increase markedly the P-R interval [10,11] and (even at low concentration) the duration of the QRS complex [10,11,12] but the duration of the J-T interval is almost not changed [10,12]. 

Conclusions

I just want to add a few comments to this model:

Acknowledgments

References

[2] Luo, C.-H. and Rudy, Y., "A model of the ventricular cardiac action potential depolarization, repolarization, and their interaction", Circulation Research, vol. 68, no. 6, pp. 1501-1526, 1991.

[3] Hodgkin, A. L. and Huxley, A. F., "A Quantitative Description of Membrane Current and its Application to Conduction and Excitation in Nerve", Journal of Physiology, vol. 117, pp.500-544, 1952.

[4] Starmer, C. F., Lastra, A. A., Nesterenko, V. V. and Grant, A. O., "Proarrhythmic response to sodium channel blockade - theoretical model and numerical experiments", Circulation, vol. 84, pp. 1364-1377, 1991.

[5] Gulrajani, R. M., "Models of the electrical activity of the heart and computer simulation of the electrocardiogram", CRC Critical Reviews in Biomedical Engineering, vol. 16, no. 1, pp. 1-66, 1988.

[6] Thiry, P. S. and Rosenberg, R. M., "On electrophysiological activity of the normal heart", Journal of the Franklin Institute, vol. 297, pp. 377-396, May 1974.

[7] Durrer, D., van Dam, R. T., Freud, G. E., Janse, M. J., Meijler, F. L. and Arzbaecher, R. C., "Total excitation of the isolated human heart", Circulation, vol. 41, pp. 899-912, June 1970.

[8] Burgess, M. J., Green, L. S., Millar, K., Wyatt, R. and Abildskov, J. A., "The sequence of normal ventricular recovery", American Heart Journal, vol. 84, no. 5, pp. 660-669, Nov. 1972.

[9] van Dam, R. T.,  and Durrer, D., "Experimental study on the intramural distribution of the excitability cycle and on the form of the epicardial T wave in the dog heart in situ", American Heart Journal, vol. 61, pp. 537-542, 1961.

[10] Somberg, J. C., "Clinical use of class Ic antiarrhythmic drugs", in Antiarrhythmic Drugs, E. M. Vaughan Williams and T. J. Campbell, Eds., vol. 89 of Handbook of Experimental Pharmacology, chapter 11, pp. 235-277. Springer Verlag, 1989.

[11] Gülker, H., Haverkamp, W. and Hindricks, G., "Leitfaden zur Therapie der Herzrhythmusstörungen", de Gruyter, 2 edition, 1992.

[12] Vaughan Williams, E. M., "Classification of antiarrhythmic actions", in Antiarrhythmic Drugs, Vaughan Williams, E. M. and Campbell, T. J., Eds., vol. 89 of Handbook of Experimental Pharmacology, chapter 2, pp. 45-67. Springer Verlag, 1989.