## Abstract

Recent experimental and theoretical results have stressed the importance of modeling studies of reentrant arrhythmias in cardiac tissue and at the whole heart level. We introduce a six-variable model obtained by a reformulation of the Priebe-Beuckelmann model of a single human ventricular cell. The reformulated model is 4.9 times faster for numerical computations and it is more stable than the original model. It retains the action potential shape at various frequencies, restitution of action potential duration, and restitution of conduction velocity. We were able to reproduce the main properties of epicardial, endocardial, and M cells by modifying selected ionic currents. We performed a simulation study of spiral wave behavior in a two-dimensional sheet of human ventricular tissue and showed that spiral waves have a frequency of 3.3 Hz and a linear core of ∼50-mm diameter that rotates with an average frequency of 0.62 rad/s. Simulation results agreed with experimental data. In conclusion, the proposed model is suitable for efficient and accurate studies of reentrant phenomena in human ventricular tissue.

- action potential
- computer simulation
- mathematical model
- reentrant arrhythmia
- spiral wave

the history of modeling biological excitable media such as nerve and heart tissue started 50 years ago with the Hodgkin-Huxley model of the giant squid axon (19). The first model of cardiac tissue (Purkinje fibers) was proposed by Noble in 1962 (30) and consisted of four variables. During the following decades, the experimental techniques for studying the properties of the cell membrane were improved continuously, leading to new cardiac tissue models of increasing accuracy, e.g., the phase-2 Luo-Rudy (LR) model of a single guinea pig ventricular cell (26), which has 12 variables, and the models developed by Noble et al., consisting of 40–60 equations (see, for example, Ref. 31). Recently, comprehensive models have been proposed to account for specific properties of human cardiac tissue: Nygren et al. (32) and Courtemanche et al. (8) developed models of single human atrial cells, whereas Priebe and Beuckelmann (35) proposed a model of a single human ventricular cell.

The Priebe-Beuckelmann (PB) model is based on the phase-2 LR model. However, five major ionic currents, including the fast (*I*
_{Kr}) and slow (*I*
_{Ks}) components of the delayed rectifier K^{+} current, the L-type Ca^{2+} current (*I*
_{Ca}), the transient outward K^{+} current (*I*
_{to}), and the inward rectifier K^{+} current (*I*
_{K1}), are based on experimental data obtained on human myocytes. In addition, parameters of the intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}) handling were changed in such a way that simulated transients are comparable to observed experimental data on human myocytes (4). The remaining currents have been adjusted from the LR model, with their amplitude scaled to fit human cell data.

Priebe and Beuckelmann developed their model to compare the electrophysiological properties of failing and nonfailing ventricular myocytes. It can be used for accurate simulations of the ionic currents and concentrations in a single cell during electrical activity. Our objective, however, is to simulate reentrant sources of arrhythmias in two (2-D) and three dimensions (3-D), which are believed to underlie most ventricular tachycardias and ventricular fibrillation (18,21, 43). Although the PB model is based on experimental measurements in human heart tissue, we refrained from using it for the study of reentrant arrhythmias for several reasons. First, the complexity and the high number of variables make the PB model inefficient for extensive computer simulations of large pieces of tissue. Second, the PB model is a “second-generation model” in which, besides membrane potential and gating variables, ion concentrations vary in time. Second-generation models made up of ordinary differential equations describing the time dependence of membrane potential, gating variables, and ion concentrations, are inherently unstable, showing complications involving long-term drifts of ion concentrations and degeneracy of equilibria, as recently emphasized by both Arce et al. (2) and Endresen and Skarland (12).

A different approach to study reentrant sources of arrhythmias in 2-D and 3-D is using low-dimensional, FitzHugh-Nagumo-type models, e.g., the Fenton-Karma model (14), which allow one to fit some overall tissue properties using three state variables. This approach, however, does not reproduce the shape of the action potential (AP), which is believed to be important for many problems, e.g., in electrocardiography. It also does not permit the study of effects of individual ionic currents on cardiac activation.

The main issue of this paper is the development of a model of an intermediate class. It should be computationally efficient and less complex than the “second-generation” ionic models and free of their inherent instabilities. At the same time, the model should retain important properties of human ventricular tissue such as restitution properties, AP shape, and change of AP shape under variation of major ionic currents. This aim was achieved by reformulating the PB model as a six-variable model, which keeps the major ionic currents but discards some single cell processes that do not substantially affect the desired tissue properties. It should be noted that the reformulation procedure does not depend on the PB model and can be applied to any ionic model.

## MATERIALS AND METHODS

The several stages in the reformulation of the PB model are discussed for each ionic current of the reduced model. We also pay attention to the numerical accuracy and stability of the reduced model and we discuss numerical schemes used in the simulations. Throughout this paper, time is in ms, membrane potential is in mV, ion concentrations are in mM, rate constants are in ms^{−1}, current is in pA/pF, and conductance is in nS/pF. Single cell membrane capacitance, cell length, and cell diameter are 153.4 pF, 100 μm, and 22 μm, respectively (same as in the original PB model).

### Reduction of the PB Model

#### Intracellular ion concentrations.

The intracellular ion concentration handling in second-generation models can be a source of instabilities (2, 12). This also holds for the PB model, which shows long term drifts of intracellular K^{+} concentration ([K^{+}]_{i}) and intracellular Na^{+} concentration ([Na^{+}]_{i}) when paced at different frequencies (see *Stability, Ionic Drift and Computational Efficiency*). Therefore, noting that the variations of [K^{+}]_{i} and [Na^{+}]_{i}over the course of an AP are extremely small, we set [K^{+}]_{i} and [Na^{+}]_{i} to constant values of 140 and 10 mM, respectively, i.e., the initial values in the PB model.

In contrast, [Ca^{2+}]_{i} undergoes substantial changes during the AP. We did not intend to study the details of Ca^{2+} dynamics but were mainly interested in the overall properties of *I*
_{Ca} and its influence on the AP shape and restitution properties of cardiac tissue. Moreover, in the PB model, the calcium subsystem is represented by a model that mimics Ca^{2+}-induced Ca^{2+} release as an all-or-none event and reproduces Ca^{2+} transients with realistic amplitudes but is phenomenological in construction and is not sufficient to accommodate predictions of intracellular Ca^{2+}dynamics under more complex experimental conditions. Therefore, we decided to assign a constant value of 400 nM to the free [Ca^{2+}]_{i} appearing in the Ca^{2+}current equations of the original PB model, as an average of its value during the course of an AP. It should be noted that we do not claim that the intracellular calcium handling is irrelevant for cardiac activity but that the simplified equations for*I*
_{Ca} listed below are sufficient for a correct representation of AP shape and restitution of AP duration (APD).

#### Ionic currents.

##### fast na^{+} current.

In the PB model, the kinetic behavior of the fast Na^{+}current (*I*
_{Na}) is described by the very fast activation variable *m* and the fast and slow inactivation variables *h* and *j*, respectively
Equation 1where the superscript o denotes an equation used in the original PB model, *V*
_{m} is the membrane potential,*g*
_{Na} is the maximum sodium conductance, and*E*
_{Na} is the Nernst potential (similar notations are used in the equations for the other currents; i.e., *g* is conductance and *E* is potential).

The activation of *I*
_{Na} is the fastest process during an AP. Hence it is tempting to eliminate *m*adiabatically. However, such reduction induces an unrealistically steep upstroke of the AP and drastically decreases the numerical accuracy of the model (see below). Therefore, we refrained from changing*I*
_{Na} activation.

Inactivation of *I*
is modeled by the product of two variables, *h* and *j*. For our purposes, a single variable *v* of the Hodgkin-Huxley type (19) suffices. Retaining the quadratic dependence of Na^{+} current inactivation, we obtain
Equation 2The time course of *v* is determined by
Equation 3where *t* is time and τ is the time constant of relaxation. The constraint we put on *v* is that the voltage-clamp behavior of *I*
should be well reproduced by the above-formulated *I*
_{Na}. Thus our aim was to obtain a reasonable approximation of*I*
(*t*) at all times for each clamped potential. The kinetics of the gating variable *v* can be determined in two steps. First, we determine*v*
_{∞}(*V*
_{m}), which is readily found from *Eqs. 1
* and *
2
*
Equation 4Next, we determine τ_{v}(*V*
_{m}) by simulating a voltage-clamp experiment and investigating the time course of*I*
_{Na} at fixed membrane potentials over a time interval [0, *T*
_{e}]. To obtain a reasonable fit for all *t*, we used a least-squares method by minimizing the function Δ_{Na}[τ_{v}(*V*
_{m})], given by
Equation 5Because the Na^{+} current is a rapid current, we performed a voltage clamp during *T*
_{e} = 100 ms for 10,000 equidistant values of membrane potential within the interval [−100 mV, 90 mV].

The obtained functions,*v*
_{∞}(*V*
_{m}) and τ_{v}(*V*
_{m}), are shown in Fig. 1, *A* and *B*, respectively. In subsequent simulations, we used values directly tabulated from the voltage-clamp simulations, but they can be well approximated by *Eqs. 16
* and *
17
* (see
).

Figure 1, *C* and *D*, shows simulated voltage-clamp records of the original and reformulated *I*
_{Na}, illustrating that at both test potentials the two closely match, both in amplitude and time. Similar results were obtained at other test potentials between −70 and +50 mV (data not shown).

##### CA^{2+} CURRENT.

Priebe and Beuckelmann modeled *I*
_{Ca} using the Hodgkin-Huxley-type (19) activation and inactivation variables *d* and *f*. In addition, they introduced [Ca^{2+}]_{i}-dependent inactivation by the*V*
_{m}-independent factor *f*
_{Ca}
Equation 6
The activation of *I*
_{Ca} is a fast process. Hence, we eliminated *d* adiabatically and obtained
Equation 7in which *f*
_{Ca} has a constant value of 0.6, because [Ca^{2+}]_{i} was set to 400 nM in the reduced model, as set out above.

##### TRANSIENT OUTWARD K^{+} CURRENT.

The equation that Priebe and Beuckelmann used for*I*
_{to} is as follows
Equation 8where *r* is a very fast activation variable and*to* is an inactivation variable. Eliminating *r*adiabatically yields our reduced current
Equation 9where *g*
_{to} was increased by 33% (see below).

##### DELAYED RECTIFIER K^{+} CURRENT.

In the PB model, K^{+} current (*I*
_{K}) consists of two components, *I*
_{Kr} and*I*
_{Ks}, each with a single activation variable,*X*
_{r} and *X*
_{s}, respectively
Equation 10
Using the same approach as for *I*
_{Na}but with *T*
_{e} set to 300 ms, we combined*I*
_{Kr} and *I*
_{Ks} into a single*I*
_{K} with a single activation variable*X* and reversal potential *E*
_{K} for given*g*
_{Kr} and *g*
_{Ks}
Equation 11For the parameter values of the original PB model, the resulting*X*
_{∞} and τ_{X} are plotted in Fig. 2, *A* and *B*, respectively. Continuous approximations are listed in the
as *Eqs. 34-36
*. The value of*g*
_{K} was determined by the constraint that*X*
_{∞} should range between 0 and 1. To obtain good restitution of APD, we increased τ_{X}slightly around the resting potential by adding a sigmoidal function τ′_{X} (see the
). This minor change had no marked influence on the shape of the AP or on the time course of *I*
_{K} (data not shown). Note that using a single *I*
_{K} does not prevent our model to be applied for, e.g., studies of drug action on the separate rapid and slow components of the K^{+} current. Therefore, one should just reapply the above fitting procedure with different initial values of *g*
_{Kr} and *g*
_{Ks}(see results).

Figure 2, *C* and *D*, shows simulated voltage-clamp records of the original (sum of*I*
_{Kr} and *I*
_{Ks}) and reformulated *I*
_{K}. In the reformulated model,*I*
_{K} shows slightly faster kinetics than in the original PB model, resulting in a somewhat larger current at the end of the 300-ms pulse.

##### TIME-INDEPENDENT CURRENTS.

For each of the time-independent currents except*I*
_{K1}, we used the original equations from the PB model, which appear in the
as *Eqs.43-50
*. For *I*
_{K1}, we used the original*Eqs. 39-42
* (see the
) but increased its conductance (*g*
_{K1}) by 56% (see below).

##### PARAMETER VALUES.

The differences between the parameter values of the original PB model and our reduced model have all been listed above. We did not change any of the other parameters. The values of all parameters of the two models are compared in Table 1.

### Numerical Approach

#### Mathematical model of wave propagation.

The propagation of an AP is modeled by the (2-D) cable equation as follows
Equation 12
where *C*
_{m} is the cell capacitance per surface unit, *S* is the surface-to-volume ratio, and ρ is the cellular resistivity. *I*
_{ion} is the net membrane current, which consists of the aforementioned*I*
_{Na}, *I*
_{Ca},*I*
_{to}, *I*
_{K}, and*I*
_{K1} as well as background Na^{+} and Ca^{2+} currents (*I*
_{Na,b} and*I*
_{Ca,b}, respectively) and the net currents generated by the electrogenic Na^{+}-K^{+} pump and Na^{+}/Ca^{2+} exchanger (*I*
_{NaK}and *I*
_{NaCa}, respectively). For*C*
_{m} = 2.0 μF/cm^{2} and*S* = 0.2 μm^{−1} (same value as in the PB model) and ρ_{x} = ρ_{y} = 180 Ω · cm (isotropic case), we obtained a plane wave conduction velocity (CV) of 70 cm/s. The value was the same as that reported by Jongsma and Wilders (22) for longitudinal CV in a linear cable of PB model cells at ρ = 181 Ω · cm.

In the one-dimensional case (no derivative with respect to*y*) we used the Crank-Nicholson as well as forward Euler methods to integrate (*Eq. 12
*). In the 2-D case, we used the operator-splitting method to split *Eq. 12
* into an ordinary differential equation (ODE) for reaction and a partial differential equation (PDE) for diffusion. The PDE was solved using an alternating-direction implicit scheme (34), whereas for the ODE we used a time-adaptive forward Euler scheme with two possible time steps, Δ*t*
_{min} = 0.02 ms and Δ*t*
_{max} = 0.1 ms. In this scheme, for each point of the medium, the ODE is solved with Δ*t*
_{max}; if d*V*
_{m}/d*t* ≤ 1 V/s, we move on to the next point; otherwise, we integrate the ODE at that point five times with Δ*t*
_{min} from its initial value. The relaxation equations for the gating variables were integrated using techniques presented by Rush and Larsen (39). All simulations were coded in C++ and run either on a PC with an Intel Pentium III 500-MHz CPU or on a Compaq AlphaServer DS20E with a 666-MHz Alpha 21264 CPU.

#### Numerical accuracy.

We tested the numerical accuracy of our model in a cable by changing the time and space step and measuring the CV of a propagating AP. We used a Crank-Nicholson method as well as a forward Euler method, yielding similar results (data not shown). Table2 shows the results we obtained. In both the original and reduced model, an increase of the time step from 0.01 to 0.02 ms caused a change of <1%, whereas an increase in Δ*x* from 0.01 to 0.03 cm reduced CV by 10% [comparable to results obtained by Qu et al. (37)]. Table 2 also show the accuracy of the reduced model upon adiabatical elimination of*m*. The >50% decrease in CV that occurred when Δ*x* was increased from 0.01 to 0.03 cm readily explains why we refrained from eliminating *m* adiabatically.

#### Restitution of APD and CV.

Restitution of APD was obtained using an S1-S2 protocol on a single cell, consisting of 10 S1 stimuli at 1-Hz frequency and an extrastimulus delivered at some diastolic interval (DI) after the last S1 AP. The restitution curve was then obtained by plotting the APD of the extra AP versus DI. Stimulus current lasted 2 ms, and its strength was two times threshold. The APD was defined using a threshold potential of −76 mV, corresponding to 90% repolarization.

Restitution of CV was measured by periodically pacing one end of a one-dimensional cable of cells. The CV was measured at a point in the middle, and, by increasing the pacing rate, we were able to measure CV at various DIs. Simulations were performed using Crank-Nicholson integration in a cable of 4-cm length with Δ*x* = 0.10 mm and Δ*t* = 0.02 ms.

#### Spiral waves and electrograms.

We obtained spiral waves in a 2-D sheet of ventricular tissue using an S1-S2 protocol. We first paced one side of the tissue (S1), producing a plane wave propagating in one direction. When the refractory tail of this wave crossed the middle of the medium, we moved our pacing electrode to that site and applied a single stimulus (S2) parallel to the S1 wavefront but not over the whole width of the medium. Stimulus currents lasted for 2 (S1) and 5 ms (S2), and their strength was two times threshold for both. The S2 wave front curled around its free end and produced a spiral wave. We calculated the trajectory of the spiral tip using an algorithm presented by Fenton and Karma (14) along an isopotential line of −60 mV.

We simulated electrograms of the spiral wave activity in 2-D by calculating the dipole source density of the membrane potential*V*
_{m} in each element, assuming an infinite volume conductor (33). We recorded the electrogram with a single electrode located 10 cm above the center of the sheet of ventricular tissue.

## RESULTS

### AP and Main Ionic Currents

Figure 3 shows a steady-state AP and the main (gated) ionic currents in the reduced and original PB models (solid and dotted lines, respectively) obtained by pacing a single cell at 1 Hz with a 2-ms, 3-nA depolarizing current. The AP in the reduced model is almost identical to the AP in the original model (Fig.3
*A* and Table 3). We observed an almost perfect match between *I*
_{Na} in reduced and original models (Fig. 3
*B*). Peak*I*
_{Ca} increased by 40%, but the global time course was preserved (Fig. 3
*C*). The time course and amplitude of *I*
_{K} were similar in the reduced and original models (Fig. 3
*D*). In the reduced model,*I*
_{to} shows a significantly larger peak and an earlier inactivation compared with the original model (Fig.3
*E*). However, the net charge carried by this current in the reduced model only increased by ∼40% compared with the original model (data not shown). Except for the large amplitude of*I*
_{to} (see below), the AP and ionic currents in the reduced model agreed well with those in the PB model. Moreover, we investigated the relative influence of each gated current on the AP shape by gradually decreasing the conductances and obtained results that closely match those obtained by Priebe and Beuckelmann (35) for the full model (data not shown).

### Restitution of APD

The APD restitution is known to play a crucial role in the formation and stability of reentrant wave patterns like spiral waves (17, 37). Figure 4A shows the APD restitution curves of the reduced model and the original model as well as the experimental data obtained from the human right ventricle by Morgan et al. (29). The two model curves are quite similar; the reduced model matching the experimental data even slightly better than the original model. Figure 4, *B* and*C*, illustrates the changes in shape of the AP with decreasing DI. At short DI, we observed a lowering of the plateau in both models. We conclude that the restitution properties of the PB model are well reproduced in our reduced model, both with respect to AP shape and APD, and that the reformulated model nicely reproduces available experimental data.

### Restitution of CV

The activation properties of the original PB model and our reduced model are compared in Fig. 5. Figure5
*A* shows the stimulus strength-duration curves of both models, whereas Fig. 5
*B* shows the CV restitution curves of the PB model. The decrease in CV with decreasing DI is due to the decreasing amount of *I*
_{Na} that is available for (re)activation. In the full model, this amount is determined by the fast and slow inactivation variables *h* and *j*, whereas in our model it is solely determined by *v*. The curves behave similarly at long DI, but at short DI the reduced model shows a slightly steeper dependence than the original model, caused mainly by τ_{v} being slightly “too small.” We did not try to get a perfect match with the CV restitution curve in the full model, because there exist no reliable experimental data on CV restitution in normal human ventricular tissue. Note, however, that the CV restitution properties can easily be adapted to new experimental data by minor modifications of τ_{v}(*V*
_{m}).

### AP Heterogeneity

The transmural heterogeneity of the ventricular AP that has been observed in animal and human myocardium (1, 24, 25, 27) plays an important role in the electrocardiogram (ECG) transcription, especially with respect to the J and T wave. As we intend to investigate and simulate ECGs of ventricular arrhythmias, it is important to take this heterogeneity into account. Three main AP morphologies are observed through the ventricular wall: the epicardial AP with the typical notch and spike-and-dome shape, the M cell AP with a less pronounced notch and prolonged repolarization phase, and the endocardial AP with a more triangular shape (1, 24, 25). The prolonged APD of the M cells is associated with a reduced*I*
_{Ks} (25), whereas the change from a spike-and-dome to a more triangular shape is due to a transmural gradient in the amplitude and kinetics of *I*
_{to}. From epicardium to endocardium, *I*
_{to} shows a decrease in amplitude, a negative shift in half-inactivation voltage (*V*
_{shift}), and an increase in reactivation time (24).

As a test on the role of selected ionic currents, we investigated how the reformulated model could reproduce the observed AP heterogeneity. The full PB model and our reduced version thereof show a clear notch and spike-and-dome morphology and were assumed to resemble an epicardial cell. To create an M cell model, we changed both*I*
_{to} and *I*
_{K} (through its slow component *I*
_{Ks}). To create an endocardial cell model, we only made changes to *I*
_{to}.

The amplitude of *I*
_{to} was changed by altering*g*
_{to}. We increased the inactivation time of*I*
_{to} by multiplying the gating functions α_{to}(*V*
_{m}) and β_{to}(*V*
_{m}) with a constant factor of *p* < 1, thus increasing τ_{to}while preserving the function to_{∞}. The half-inactivation voltage was changed through a shift in the gating functions of*to* by *V*
_{shift}. The resulting equations for τ_{to}(*V*
_{m}) and to_{∞}(*V*
_{m}) are listed in the
as *Eqs. 31
* and *
32
*. To prolong the APD of the M cells, we decreased *g*
_{Ks}by 50% in the original PB model and performed new voltage-clamp simulations on the sum of *I*
_{Kr} and*I*
_{Ks} to obtain *X*
_{∞}, τ_{X}, and *g*
_{K} for the M cell. The thus-obtained tabulated values of *X*
_{∞}and τ_{X} are well approximated by the functions listed in the
as *Eqs. 37
* and *
38
*. The new *g*
_{K} was found to equal 0.013 nS/pF.

The AP heterogeneity was obtained by using the parameter values listed in Table 4. Note that all relative changes are well within the experimentally observed range (24). The simulated APs for epicardium, M cells, and endocardium are shown in Fig. 6,*A–C*, and have a duration of 360, 362, and 400 ms, respectively, and are similar to the APs observed experimentally (1, 24). Figure 6
*D* shows the APD restitution curves for the different cell types. Like those reported by Antzelevitch et al. (1), the restitution curves for epicardium and endocardium almost coincide, whereas the restitution curve for the M cells is considerably steeper.

### Spiral Waves

The 2-D simulations were performed on a 700 × 700 square lattice with Δ*x* = 0.25 mm. We obtained spiral waves using an S1-S2 stimulation protocol, as described in materials and methods. The results of this computation are shown in Fig. 7. After application of the S2 stimulus (Fig. 7
*A*), the spiral tip drifts along the refractory tail of the wave (Fig. 7
*B*) before curling back (clockwise) and stabilizing approximately at the middle of the medium (Fig. 7
*C*). Figure 7
*D* shows the trajectory of the tip after stabilization. The tip motion occurs along lines that are ended by a quick turn through ∼180°. The linear part of the trajectory rotates counterclockwise in space at an angular velocity of ∼0.62 rad/s. This tip motion is characteristic of the spiral waves with a linear core observed in experiments on thin slices of myocardium (9) and in the whole heart (11, 16) as well as in mathematical models of cardiac tissue (10, 13).

Figure 8
*A* shows a record of membrane potential in a point far from the spiral core (asterisk in Fig. 7
*C*) during 10 s after S1. Due to restitution, the APD changed from ∼350 to ∼280 ms. The average period of the spiral wave was 304.48 ± 3.09 ms (means ± SD). Figure8
*B* shows the ECG associated with the spiral wave. It is similar to ECGs obtained during polymorphic ventricular tachycardia. The Fourier transform of the ECG shows a dominant frequency around 3.29 Hz (not shown), which is comparable to the value of 3.65 ± 0.46 Hz obtained in clinical studies (7).

### Stability, Ionic Drift, and Computational Efficiency

We compared the speed of computations by pacing a cell at 1 Hz during 1,000 s of model time. All computations were performed on the 666-MHz Alpha 21264 CPU of a Compaq AlphaServer DS20E using Euler-type integration with a time step of 10 μs. For the original model, the simulation took 386.5 s of computer time, whereas for the reduced model, the simulation finished in 79.5 s. During this simulation, the original model showed a drift in both [Na^{+}]_{i} and [K^{+}]_{i}. Over 1,000 s of model time, mean [Na^{+}]_{i}decreased from 10 to 8.1 mM (Fig.9
*B*), whereas the mean [K^{+}]_{i} decreased nonmonotonically from 140 to 139.5 mM (Fig. 9
*C*). As a consequence, the PB model did not reach a steady state, and, after 1,000 s of model time, the APD ran up to 425 ms (Fig. 9
*A*).

Because [Na^{+}]_{i} and [K^{+}]_{i} change very little over the course of one AP, it is tempting to eliminate their drifts by putting them to a constant value. If the above simulation of 1,000 s of model time is repeated with [Na^{+}]_{i} and [K^{+}]_{i} set to constant values of 140 and 10 mM, respectively, the PB model reaches a steady state with an APD of ∼360 ms and the simulation takes 305 s of computer time (data not shown). However, if the cycle length is abruptly changed, e.g., to 2 Hz, a very slow adaptation of [Ca^{2+}]_{i} is observed (Fig. 10
*B*) and it takes >1 min of model time for the cell to reach a new steady state (Fig. 10
*A*).

Both these observations are unrealistic because paced isolated cells adapt in a few beats to a given cycle length. They were the main reason why we refrained from using the original PB model, even with [Na^{+}]_{i} and [K^{+}]_{i}set to constant values, for extensive simulations of a 2-D spiral wave. In those cases, it would take many rotations before all transient drifts would have disappeared. In the reformulated model, the cell adapts after a few beats to a given cycle length, even after an abrupt change (Fig. 10
*C*), just like real isolated cells. Given the increase in computational speed obtained with the reformulated model (an increase by a factor of 4.9 compared with the full PB model and a factor of 3.8 compared with the PB model with constant [Na^{+}]_{i} and [K^{+}]_{i}) as well as the much faster and more realistic adaptation to changes in cycle length, it is clear that the use of the reduced model saves a huge amount of computational time when investigating reentrant arrhythmias. Also, with the number of model variables reduced from 15 to 6, the demand for computer memory is reduced by a factor of 2.5.

## DISCUSSION

### Action Potential

The results of our simulations show that the overall properties of the AP in the original PB model are well reproduced in the reduced model. The shape of the AP is similar in both models (Fig. 3) and AP parameters are preserved (Table 3). To have comparable resting membrane potential and APD, we had to increase *g*
_{K1} by 56%. With a maximum outward current in the negative voltage range of 0.36 pA/pF at −60 mV, Beuckelmann et al. (5) found quite low values of *I*
_{K1} in human ventricular myocytes, whereas higher values, up to 8 pA/pF (23), have been reported in other studies. In our reduced model, with increased*g*
_{K1}, the maximum of *I*
_{K1}in the negative voltage range is 1.33 pA/pF, which is well within the range of experimental data.

### Main Gated Currents

The overall time course of the individual ionic currents resembled that in the original model (Fig. 3). The amplitudes of*I*
_{to} and *I*
_{Ca} are, however, larger than in the PB model, due to the adiabatic elimination of the activation variables *r* and *d*, respectively. The amplitude of *I*
_{to} is furthermore enhanced by a 33% increase in *g*
_{to}. This increase was necessary to preserve the typical notch that follows the AP upstroke. Although these amplitudes may seem less realistic, we found that they had no influence on the properties of the AP we studied.

### Restitution Properties

The restitution properties of the original PB model were obtained far from steady states, where large ionic drifts occur (see Fig. 9). Therefore, we did not solely compare our model with the original PB model but also with available experimental data. The reformulated model closely fits the data obtained by Morgan et al. (29) (Fig.4
*A*). Experimental studies of restitution in human myocardium have also shown biphasic restitution curves (15, 28). This phenomenon is still poorly understood and, so far, no ionic models have been developed to reproduce these nonmonotonic curves. Biphasic restitution curves can easily be obtained by introducing a DI dependency of a model parameter (e.g., *g*
_{Ca}) (44). Such approach is, however, detached from the underlying electrophysiology. Therefore, we decided to wait until more experimental data on this property become available before trying to modify the APD restitution curve.

The CV restitution is steeper than in the full model (Fig.5
*B*). Because there are no experimental data on CV restitution in human myocardium, we did not attempt to flatten the CV curve in the reduced model. Moreover, when Saumarez et al. (40,41) measured the latency of propagating pulses in human ventricular tissue, they found that normal ventricular tissue (less susceptible to cardiac arrhythmias) showed a steeper increase in latency at short DI than the tissue of patients who had cardiac arrhythmias. Therefore, the steeper curve in the reduced model may be an even more accurate model of normal ventricular tissue.

### AP Heterogeneity

On the basis of experimental findings on*I*
_{to} and *I*
_{Ks} (24,25), we were able to reproduce the observed AP heterogeneity in ventricular muscle. We obtained an APD prolongation in the M cells through a 50% block of *I*
_{Ks} like reported by Liu et al. (25). The simulated APD difference between M cells and epicardium and endocardium is, however, smaller than the difference observed by Li et al. (24) in myocytes isolated from the human right ventricle. This could indicate that other currents are involved in the APD prolongation of M cells. Indeed,*I*
_{K1} seems to be somewhat smaller in M cells from the canine ventricle (25), and transmural inhomogeneities in other currents are starting to be reported as well (3, 45,46). On the other hand, quite small APD differences were measured in guinea pig ventricular tissue by Main et al. (27). These observations reflect species differences, and, more recently, differences between right and left ventricles have also been observed (42). For the moment, we feel that our simulations of the observed transmural AP heterogeneity, based on differences in *I*
_{to} and*I*
_{Ks}, are realistic enough, but we are aware that we may have to include transmural inhomogeneities in other currents as well, once these are demonstrated experimentally in the human ventricle.

### Numerical Efficiency

We were able to reformulate the PB model with five gating variables by altering only two electrophysiological parameters,*g*
_{to} and *g*
_{K1}. The total number of variables has been reduced from 15 in the original model to 6 in the reformulated model, i.e., the membrane potential*V*
_{m} and the gating variables of*I*
_{Na} (*m* and *v*),*I*
_{Ca} (*f*), *I*
_{to}(to), and *I*
_{K} (*X*). The computational gain at the cellular level was of the order of a factor of four. The computational speed can be increased further by the choice of sophisticated numerical schemes, such as adaptive space masks (6) or adaptive time step schemes (36). We used a scheme similar to the latter in 2-D simulations of reentrant arrhythmias. Furthermore, there was a considerable gain in stability and rate adaptation. Moreover, in the reduced model, it was less difficult to relate changes in parameter values to changes in the AP. Taken together, the reformulated model is very efficient for extensive 2-D and 3-D simulations of reentrant arrhythmias.

### Spiral Waves in 2-D

The spiral waves generated with the reduced model show meandering dynamics, preceded by a short transient drift, and have a period of ∼300 ms. The ECG of the 2-D spiral waves is polymorphic and supports the idea that meandering reentrant arrhythmias underlie polymorphic ventricular tachycardias (18, 21, 43). We are aware of the quantitative limitations of 2-D simulations of spiral waves and associated ECGs. The geometry and anisotropy of the ventricles can have an important influence on the dynamics of vortexes and the electrical properties of the human body have to be taken into account. The AP heterogeneity addressed above is also important for a correct description of heart repolarization, e.g., in modeling the J and T waves of the ECG. The 2-D spatial domain in which we study spiral waves is quite large compared with the size of the human heart (17.5 × 17.5 cm). However, our computations show that once a spiral is formed, it can rotate in a medium as small as 8.2 × 8.2 cm. In the whole heart, spiral waves may be formed because of geometry and anisotropy.

### Limitations of the Model

The inherent limitations of the full PB model have been discussed in the paper by Priebe and Beuckelmann (35). The main limitation of our reformulated model is the lack of ionic concentration handling. Giving constant values to all intracellular concentrations is a drastic way of reducing the number of equations and puts some constraints on the applicability of the model; the dynamics of calcium overload, sodium overload, etc., can no longer be studied. It should, however, be noted that the equations representing the calcium subsystem are still evolving and even the latest comprehensive models fail to adequately represent fundamental properties of calcium handling. Jafri et al. (20) incorporated an entirely new description of intracellular Ca^{2+} handling into the phase-2 LR model, resulting in a set of as many as 30 ordinary differential equations that provide an improved description of the L-type Ca^{2+} channels and Ca^{2+}-induced Ca^{2+}release from ryanodine receptor channels exhibiting adaptation. Whereas this model reproduces important aspects of cardiac Ca^{2+}cycling, it also exhibits all-or-none rather than graded Ca^{2+} release (38).

In summary, we developed an efficient electrophysiological model of single human ventricular cells that may facilitate and improve large scale computational studies of reentrant arrhythmias.

## We are thankful to Kirsten H. W. J. ten Tusscher for help in 2-D computations and Dr. S. McNab for editorial advice.

This research was supported by a specialization grant from the Flemish Institute for the Promotion of Scientific and Technological Research in Industry and by the Netherlands Organization for Scientific Research via the Research Council for Earth and Life Sciences.

## Model Equations

### Equilibrium Potentials

where *R* is the universal gas constant, T is the absolute temperature, *F* is the Faraday constant, and [Na^{+}]_{e}, [Ca^{2+}]_{e}, and [K^{+}]_{e} are the extracellular Na^{+}, Ca^{2+}, and K^{+} concentrations, respectively.

### Inward Currents

#### Fast Na^{+} current.

Equation 13 Equation 14 Equation 15 Equation 16 Equation 17

#### Slow Ca^{2+}* current.*

Equation 18 Equation 19 Equation 20 Equation 21 Equation 22 Equation 23 Equation 24

### Outward Currents

#### Transient outward current.

Equation 25 Equation 26 Equation 27 Equation 28 Equation 29 Equation 30 Equation 31 Equation 32

#### Delayed rectifier K^{+}* current.*

Equation 33For endocardial and epicardial cells Equation 34 Equation 35 Equation 36For M cells Equation 37 Equation 38

#### Inward rectifier K^{+} current.

Equation 39 Equation 40 Equation 41 Equation 42

### Background Currents

#### Ca^{2+}* background current.*

Equation 43

#### Na^{+}* background current.*

Equation 44

### Pump and Exchanger Currents

#### Na^{+}*-K*^{+}pump.

^{+}pump.

Equation 45 Equation 46 Equation 47 Equation 48

#### Na^{+}*/Ca*^{2+}exchanger.

^{2+}exchanger.

Equation 49 Equation 50

## Footnotes

Address for reprint requests and other correspondence: O. Bernus, Dept. of Mathematical Physics and Astronomy, Krijgslaan 281 (S9), 9000 Gent, Belgium (E-mail:olivier.bernus{at}rug.ac.be).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.First published January 24, 2002;10.1152/ajpheart.00731.2001

- Copyright © 2002 the American Physiological Society