A Probabilistic Model of Cardiac Electrical Activity Based on a Cellular Automata System

INTRODUCTION

Computer simulation models of myocardial electrical activation have been widely used in cardiology. In the field of electrophysiology, several types of models have provided a deeper understanding of the mechanisms implicated in arrhythmias,^1,2 the properties of sinoatrial node pacemaker activity,³ conduction,^4,5 and reentry circuits.^6-8 More complex phenomena such as fibrillatory activity have also been investigated with these techniques,^9,10 and it has even been possible to predict some of the effects of antiarrhythmic agents on ventricular fibrillation^11,12 and other arrhythmias.^13,14 The main activation models can be classified into 2 major groups. The first is based on the cellular automata concept, which considers the tissue to be a set of discrete elements connected to each other. Each element, or automaton, can have a finite number of allowable states (e.g., active or inactive), each of which varies as a function of the preceding state and the state of the neighboring elements according to a set of predefined rules. The transitions between individual states govern the evolution and behavior of the system as a whole and are deterministic, since the system evolves in an identical manner under the same initial conditions. Cellular automata are easy to program and allow fast simulations with a moderate computational burden; nevertheless, they have serious limitations when attempting to reproduce various phenomena of significant interest in electrophysiology, such as the effects of curvature on the activation fronts.¹⁵ The second group of models, known as reaction-diffusion models, reproduce membrane dynamics by simulating the flow of ionic currents through the ion channels and the associated membrane potential.^16-18 The connection between the cells is modeled by electrical resistors and the propagation of the impulse is calculated by solving the equations used to describe electrical circuits.^19,20 These models are much more realistic and more accurately reproduce complex situations, but require an extremely high computational burden and often need many calculation hours in multiprocessors to simulate 1 or 2 seconds of electrical activity.²¹ The purpose of this study was to design a computer model of electric activation based on a cellular automata system that simulates complex electrophysiological phenomena and does not require the computational burden needed for reaction-diffusion models.

MATERIALS AND METHODS

Model Design

The cardiac tissue was modeled as a grid of discrete elements, or cells, that represent groups of cells with an average intrinsic behavior, and interact with the neighboring cells according to a probabilistic rule of electrical impulse transmission. The behavior of each cell resembles a cellular automaton that can have 3 states: resting (relaxed and excitable), refractory 1 (excited, able to excite neighboring cells), and refractory 2 (excited, but not able to excite neighboring cells). The refractory 1 period is maintained for a fraction F (around 10%) of the duration of the action potential. During the rest of the action potential, the cell remains in refractory 2 period, then changes to the resting value, which corresponds to the diastolic interval. The transition between states is governed by three laws: partial repolarization (transition from refractory 1 to refractory 2), total repolarization (transition from refractory 2 to rest), and depolarization (transition from rest to refractory 1). Both partial repolarization and total repolarization take place in a deterministic manner, when both the instant of depolarization and the action potential duration are known. In contrast, depolarization is defined in probabilistic terms and is based on two factors: a) the cell excitability (E), which increases with the time the cell remains at rest, and b) the excitation quantity around each cell (Q), in which the greater the excitation quantity, the greater the probability that the element is excited. If P_j^exc denotes the probability that cell j is excited, then these two factors can be combined in the following formula:

where i is an adjacent element, A_i is the binary excitation state (1 if refractory, 0 if not), and d_ij is the distance between the midpoints of elements I, and j. Note that the formula can be applied to both two- and three-dimensional substrates, by extending the summatory function to the neighboring elements on a surface (flat or curved) or to a three-dimensional space. Both factors (excitability and excitation quantity) were estimated using 2 macroscopic variables, action potential duration and conduction velocity, taking into account the electrical restitution properties of cardiac tissue. According to these properties, the action potential duration and conduction velocity depend on the frequency of the stimulation to which the tissue is exposed. A high frequency reduces the diastolic interval, leading to a brief action potential duration and a slow conduction velocity, whereas a low frequency produces the opposite phenomenon. This behavior is shown in the restitution curves (Figure 1) used to adjust the model parameters.

Figure 1. Restitution curves for the action potential (left) and conduction velocity (right), based on the diastolic interval. APD indicates action potential duration; CV, conduction velocity.

Calculation of the Electrograms

In order to calculate the electrograms associated with the activation standard, the transmembrane current must be calculated from the action potential. In our case, we used a prototype action potential obtained from the model devised by Luo and Rudy,¹⁷ which associates the time for which each cell has been in its current state with its voltage level. Assuming that cardiac tissue can be modeled as a medium having a single, isotropic, resistive, and homogenous domain, the electrical potential obtained at a point near the extracellular surface would be calculated as the weighted sum of the transmembrane currents that originate in the tissue cells.^22,23 A unipolar electrogram is modeled as the recording of the extracellular potential measured by a single electrode of positive polarity, with the reference (zero potential) set to infinity. The distance between the electrode and the surface quantifies the zone influenced by the electrode, such that the nearer the tissue, the greater the field capture. A bipolar electrogram is modeled as the difference in potential between two points near each other, i.e., as the difference between 2 unipolar electrograms.

The simulation code and display of the experiments was done using the Matlab^® computational software package (The MathWorks, Inc.), using an Intel^® Pentium^® 4 PC computer (2.19 GHz, single-processor, 256 Mb RAM).

RESULTS

Using a two-dimensional tissue slice of 80x120 elements with homogeneous properties as a substrate, the initial, simultaneous activation of a line of elements generated a flat activation front with a uniform propagation velocity. Figure 2 contains the bipolar and unipolar electrograms obtained from various points in space. Focal stimulation on the same substrate produced a round front that moved away from the point of stimulation. The simulated unipolar electrogram presented a QS morphology with the recording electrode located on the stimulation area, unlike the rS morphology, which was obtained at some distance from the middle. The introduction of unexcitable areas and slow conduction zones in the model generated stable and unstable reentry circuits. Figure 3 shows the simulation of a reentry circuit in a figure 8 shape, similar to those observed in experimental models and in patients with sustained ventricular tachycardia in chronic myocardial infarction.^24-26 The gray rectangles indicate necrotic areas, with the corridor between the two an area of slow conduction. During stimulation from the far left portion (top figure), the electrode located on the slow conduction corridor records an initial potential (blue arrow) that corresponds to activation of the rapid conduction area, located outside the necrotic areas. Subsequently, a late, low-amplitude potential (red arrow) is recorded separated from the initial one by an isoelectric line. This potential corresponds to slower activation of the corridor between the two scars. Similar potentials have been described in patients with ventricular tachycardia and myocardial infarction.²⁷ During reentry (Figure 3, bottom), activation of the slow conduction area generates a mid-diastolic low-amplitude potential (red arrow). The probabilistic concept of excitation in our model is particularly appropriate for simulating the behavior of activation fronts in areas with a pronounced curvature. The convex border of the front encounters more surface area to activate in the direction of activation, and therefore the probability that neighboring elements in this direction will be excited decreases. Because conduction velocity depends essentially on this parameter, it is inversely proportional to the curvature of the front (Figure 4) and reproduces the actual behavior in this situation with no need for additional ad hoc hypotheses. The model also allows simulations to be performed in fibrillation processes. Using a tissue slice of homogeneous properties and 80x80 elements as substrate, rapid stimulation with a frequency above 8-10 Hz generates curvature and extraordinary fragmentation of the wavefronts in a few seconds, which is typical of fibrillatory conduction. Figure 5 contains an example of stimulation at 50 Hz. Note how the electrograms obtained in the neighborhood of the substrate reproduce those obtained during fibrillation. Finally, the model allows three-dimensional simulation of the activation process, providing depth to the grid that sustains the cellular automata system. The three-dimensional substrate can be defined with homogeneous properties, assuming the same behavior as activation on the surface and in depth, but can also generate a heterogeneous model by modifying the calculation of the distance according to the direction. This last aspect allows different conduction velocities to be simulated in transversal and longitudinal directions and in depth, reproducing the anisotropic phenomenon. By using a homogeneous substrate, successive stimulation of a cube of 80x80x10 elements with two flat fronts in the perpendicular direction generates an activation front with a conduction velocity that decreases as it nears its rotation center. The extreme curvature of the activation wave makes it turn around the core, which remains unexcited, generating a rotor similar to those observed during atrial and ventricular fibrillation in various experimental models^28-30 (Figure 6, video 1 available in the January issue at: www.revespcardiol.org/).

Figure 2. Simulation of a flat wavefront and the respective unipolar and bipolar electrograms. From left to right, 3 time points during activation in the tissue slice are shown (1: start, 2: middle, 3: end). The positions of the recording electrode are indicated by the letters a, b and c. The unipolar (U) and bipolar (B) electrograms obtained in each position are shown at the top of the tracing. The right bar represents the color map of the membrane voltage, expressed in mV.

Figure 3. Simulation of a reentry circuit in the shape of a figure "8." The red arrows on the tissue slice indicate the direction of activation. Slow conduction is depicted by a wavy arrow. On each slice, the bipolar electrograms obtained on the corridor between the scar areas are shown (see text). The rectangle on the right represents the color scale of membrane voltage, expressed in mV.

Figure 4. Effect of front curvature on conduction velocity, indicated by white arrows. Note the decrease in velocity with the increase in curvature. The voltage map is shown as in Figures 2 and 3.

Figure 5. Fibrillatory conduction in response to high-frequency (50 Hz) stimulation of the far left portion of the tissue slice. The upper box contains the electrogram obtained when the electrode is placed on the tissue surface. The rectangle on the right represents the voltage map, as in the previous figures.

Fig. 6. Simulation of a spiral activation wave (rotor). The various colors represent the activation time, and the midpoint where they converge corresponds to the core of the rotor (see text).

DISCUSSION

The cellular activation model described is able to reproduce many aspects of the electrical behavior of myocardial tissue, including flat activation fronts, more or less complex reentry circuits, the effects of the curvature front, fibrillatory conduction and the genesis of stable rotors described in experimental models of atrial and ventricular fibrillation. These phenomena can be simulated with a cellular automata model based on simple rules to reduce the computational burden. The probabilistic nature of excitation is fundamental to the model, as it naturally reproduces propagation behavior in curvilinear fronts and allows even stable rotors to be simulated in homogenous, two- or three-dimensional media.

The published cardiac impulse activation and propagation models differ in their implementation of the macroscopic and microscopic cardiac anatomy and their approach to the electrophysiological properties of the cell and the intercellular connections. Some models attempt to simulate cellular electrophysiology by reproducing the behavior of the ionic currents that generate the transmembrane action potential. These models are based on the original equations of the action potential model in the giant squid axon of Hodgkin and Huxley,³¹ and are adapted to the cardiac tissue in subsequent modifications.^16,17,32-35 This approach is the basis for the so-called reaction-diffusion models and requires complex mathematical calculations to describe the processes of excitation and impulse propagation. These types of models have been able to reproduce the action potential of cardiac cells under normal conditions as well as in pathological situations such as ischemia or heart failure.³⁶ They are also highly suitable for assessing the electrophysiological effects of drugs that act on the ion channels.^12,37,38 The effects of curvature on the fronts simulated with these models are quite realistic, with the conduction velocity decreasing as the degree of curvature decreases. Rotors and vortices have also been simulated with these models, in both 2 and 3 dimensions.^9,15 The main drawback of reaction-diffusion models is that very long computation times in high-capacity computers are needed. Even in simplified models,^39,40 the computational needs are very high, making them suitable for simulations of small tissue specimens for a few seconds. An important advantage of our model with respect to those described is that it is simpler to program, has a much lower computational burden and can therefore simulate complex and/or long-lasting phenomena without the need for a sophisticated computing infrastructure. In fact, our evaluation was done on readily available commercial computers in reasonable time frames for complex simulations, such as those involved in fibrillation phenomena. With high-capacity computers, such as those currently used to run reaction-diffusion models, the resolution of the model would increase considerably, achieving results closer to the experimental preparation and enabling the phenomena occurring in an entire chamber of the heart to be simulated in realistic anatomic detail. The second group of activation models includes cellular automata systems and is based on representing the cardiac tissue as a grid of discrete elements, that can adopt different states as a function of their preceding state and the state of neighboring cells.⁴¹ In published models, the communication between one element and its neighbors occurs according to a deterministic rule, with each element obliged to change its state or to retain it, depending on the situation of the surrounding cells. This methodology has been able to reproduce one-directional functional blocks with induction of stable reentries around an obstacle. Spiral waves that rotate around a fixed center have also been simulated, with fragmentation achieved by introducing anisotropy in the substrate.⁴¹ In addition, the inherent calculation speed of these models means they can be used for large surfaces or volumes, simulating the behavior of complete heart chambers. Nevertheless, the reproduction of some key phenomena in the dynamics of stimulus propagation is rather unrealistic or requires the incorporation of additional, hard-to-justify hypotheses in the model. This occurs, for instance, with the decrease in conduction velocity produced when the frequency of stimulation is increased and particularly, with the behavior of curvilinear fronts, in which the conduction velocity varies inversely with the degree of curvature. This fact forms the basis of many electrophysiological observations of interest and is hard to reproduce with the available models. The introduction of the probabilistic factor in our model, whereby the excitation probability of each element depends on the time since its last repolarization and on the number of active cells nearby, reproduces both phenomena naturally, without the need for additional hypotheses and without increasing the calculation time. Another phenomenon we were able to reproduce directly and that has not been previously reported with cellular automata systems is fibrillatory conduction in the substrate when the stimulation frequency is very high. Since unipolar and bipolar virtual electrograms can be obtained at any position in space and at any angle with respect to the activation surface, the model is more useful and also establishes a relationship with the electrical tracings directly obtained in experimental models or the electrophysiology laboratory.

Limitations

The model developed is based on macroscopic properties of activation, and therefore does not directly take into account the behavior of the cellular ion channels and currents. As a result, the model in its current state would be inadequate for simulating the effects of drugs or genetic diseases that modify these currents. This limitation, common to all cellular automata models, could be at least partially overcome by introducing parameters that indirectly reflect the most relevant aspects of channel dynamics. Secondly, the simulations have been done on a standard computer with a limited calculation capacity. The use of more powerful computers would allow better resolution and probably confer a more realistic aspect to the simulations. In addition, our model was evaluated on a simple, homogeneous, two- or three-dimensional virtual substrate. For the time being, the anatomy of the chambers, the integration of anatomical obstacles and the orientation of the muscle fibers have not been included. Finally, although the model is capable of qualitatively reproducing many interesting phenomena in basic electrophysiology, we have not performed a quantitative validation in experimental preparations. This validation would allow a better fit of the model parameters and an assessment of their predictive capacity.

Funded in part by a basic research grant from Guidant Spain.

See Editorial on Pages 6-9

---

Correspondence: Dr. A. García Alberola.
Pza. del Roble, 36. 30150 La Alberca. Murcia. España.
E-mail: [e-mail].