Introduction

Membrane models have turned out to be the relevant ingredients for the simulation of electrical activity of cardiac tissue on small scales. Based on a set of specific membrane channels and associated ionic currents with a complicated, nonlinear dynamics, they give a reliable, in a wide range of experimental situations well tested microscopic description of electrical activity. Of course, membrane models could, in principle, be the starting point for the simulation of e.g. a complete ventricle; the enormous amount of computer power required, however, renders this approach rather impractical.

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

In this section the hybrid method to simulate the electrical activity of the heart is demonstrated. For this purpose the following situation is studied as an example: It is assumed that the sinoatrial node (SAN) stimulates periodically and that the heart has reached a "stationary" periodic state characterized by the stimulation period, the basic cycle length (BCL). Now, the electrical response of the heart to an "aperiodic" stimulus i.e. a stimulation by the SAN with a cycle length t_hist differing from BCL, is investigated. Obviously, this is a "one-parameter-family of situations" labeled by t_hist. To be precise, since the periodic and aperiodic excitations propagate with different velocities, t_hist will in general be a locally varying quantity.

Appropriate simulations with a linear chain of membrane model cells (details see below) and a variety of t_hist values are conducted and the results obtained are processed suitably: It turns out that the electrical response, the action potential (AP) and its velocity of excitation can satisfactorily be described by a general shape formula with five t_hist-dependent parameters (normalized results).

There arises a slight complication, which in addition is accounted for within the hybrid method. Typical electrical parameters of the AP (normal situation: no drugs, basic cycle length BCL=1000 ms, that means a heart frequency of F=1 Hz.) depends on the location in the heart. However, the corresponding adjustments of the Luo and Rudy membrane model are not known or at least not published.

Therefore, it is assumed here that the local behaviour may be described by a t_hist-independent deformation of the linear chain results. These deformations are defined by rescaling the normalized parameters such that at t_hist=BCL they are equal to those typical parameters given in textbooks.

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<(p>

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.

Results

A simulation under normal conditions ( t_hist=1000 ms, no drug, normal electrical parameters) shows the basic properties of this computer model.

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

In this paper, a new method is presented which fills the gap between membrane models and cellular automata approaches to the simulation of electrical activity. With the hybrid method the rules of the CA get a quantitative justification. So effects, which could be simulated with membrane models on small tissue samples are "translated" with the hybrid method to the electrical activity of whole heart models.
Here, the effects of the cycle length and of a (hypothetical) class I drug on the electrical activity of the whole heart are simulated. But this method is not restricted to class I medicaments: any effect which can be implemented in membrane models is a possible candidate for the hybrid method.

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

This membrane model is only valid for ventricular cells. As soon as membrane models for other types of cells (atrial cells, nodal cells, cells of the SCS, etc.) are available this restraint could be overcome.

All the results are produced for a basic cycle length BCL=1000 ms. Only for this value of the cycle length ("stationary conditions"), medical textbooks present sufficient data for comparison. To cover a wider range of frequencies, all the simulations (membrane model, CA) have to be conducted with other values of the BCL.

The real human heart, the geometric data set was derived from, seems to exhibit right hypertrophic properties.

Acknowledgments

Financial support by the Max-Planck-Gesellschaft is gratefully acknowledged.

References