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
 |
INTRODUCTION |
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 (IKr) and slow (IKs)
components of the delayed rectifier K+ current, the L-type
Ca2+ current (ICa), the transient
outward K+ current (Ito), and the
inward rectifier K+ current (IK1),
are based on experimental data obtained on human myocytes. In addition,
parameters of the intracellular Ca2+ concentration
([Ca2+]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, [Ca2+]i undergoes substantial
changes during the AP. We did not intend to study the details of
Ca2+ dynamics but were mainly interested in the overall
properties of ICa 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
Ca2+-induced Ca2+ release as an all-or-none
event and reproduces Ca2+ transients with realistic
amplitudes but is phenomenological in construction and is not
sufficient to accommodate predictions of intracellular Ca2+
dynamics under more complex experimental conditions. Therefore, we
decided to assign a constant value of 400 nM to the free
[Ca2+]i appearing in the Ca2+
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
ICa 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 (INa) is described by the very fast
activation variable m and the fast and slow inactivation
variables h and j, respectively
|
(1)
|
where the superscript o denotes an equation used in the original
PB model, Vm is the membrane potential,
gNa is the maximum sodium conductance, and
ENa 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 INa 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
INa 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
|
(2)
|
The time course of v is determined by
![d<IT>v/</IT>d<IT>t=</IT>[<IT>v<SUB>∞</SUB></IT>(<IT>V</IT><SUB>m</SUB>)<IT>−v</IT>]<IT>/&tgr;<SUB>v</SUB></IT>(<IT>V</IT><SUB>m</SUB>)](/content/vol282/issue6/fulltext/H2296/img004.gif)
|
(3)
|
where 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 INa.
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
(Vm), which is
readily found from Eqs. 1 and 2
![v<SUB>∞</SUB>(V<SUB>m</SUB>)<IT>=</IT>[<IT>h<SUB>∞</SUB></IT>(<IT>V</IT><SUB>m</SUB>)<IT> · j<SUB>∞</SUB></IT>(<IT>V</IT><SUB>m</SUB>)]<SUP>½</SUP>](/content/vol282/issue6/fulltext/H2296/img005.gif)
|
(4)
|
Next, we determine
v(Vm) by simulating a
voltage-clamp experiment and investigating the time course of
INa at fixed membrane potentials over a time
interval [0, Te]. To obtain a reasonable fit
for all t, we used a least-squares method by minimizing the
function
Na[v(Vm)],
given by
![&Dgr;<SUB>Na</SUB>[<IT>&tgr;<SUB>v</SUB></IT>(<IT>V</IT><SUB>m</SUB>)]<IT>=</IT><LIM><OP>∫</OP><LL>0</LL><UL><IT>T</IT><SUB>e</SUB></UL></LIM> [<IT>I</IT><SUP><IT>o</IT></SUP><SUB>Na</SUB>(<IT>t</IT>)<IT>−I</IT><SUB>Na</SUB>(<IT>t</IT>)]<SUP>2</SUP> d<IT>t</IT>](/content/vol282/issue6/fulltext/H2296/img006.gif)
|
(5)
|
Because the Na+ current is a rapid current, we
performed a voltage clamp during Te = 100 ms for 10,000 equidistant values of membrane potential within
the interval [
100 mV, 90 mV].
The obtained functions,
v(Vm) and
v(Vm), 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
APPENDIX).
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 1.
A and B: steady-state value
v (A) and the time constant
v (B) of the inactivation gating
variable v of the fast Na+ current
(INa) in the reformulated Priebe-Beuckelmann
(PB) model. C and D: comparison of
INa in the reduced and the original PB model.
Activation of INa in response to a voltage-clamp
step from a holding potential of 90 mV to a test potential of 15 mV
(C) and +15 mV (D) is shown.
|
|
Figure 1, C and D, shows simulated voltage-clamp
records of the original and reformulated INa,
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).
CA2+ CURRENT.
Priebe and Beuckelmann modeled ICa using the
Hodgkin-Huxley-type (19) activation and inactivation
variables d and f. In addition, they introduced
[Ca2+]i-dependent inactivation by the
Vm-independent factor fCa
![I<SUP><IT>o</IT></SUP><SUB>Ca</SUB><IT>=g</IT><SUB>Ca</SUB><IT> · d · f · f</IT><SUB>Ca</SUB>([Ca<SUP>2+</SUP>]<SUB>i</SUB>)<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>−E</IT><SUB>Ca</SUB>)](/content/vol282/issue6/fulltext/H2296/img007.gif)
|
(6)
|
The activation of ICa is a fast process.
Hence, we eliminated d adiabatically and obtained
|
(7)
|
in which fCa has a constant value of 0.6, because [Ca2+]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
Ito is as follows
|
(8)
|
where r is a very fast activation variable and
to is an inactivation variable. Eliminating r
adiabatically yields our reduced current
|
(9)
|
where gto was increased by 33% (see below).
DELAYED RECTIFIER K+ CURRENT.
In the PB model, K+ current (IK)
consists of two components, IKr and
IKs, each with a single activation variable,
Xr and Xs, respectively
![I<SUP>o</SUP><SUB>K</SUB><IT>=I<SUB>Kr</SUB>+I</IT><SUB>Ks</SUB><IT>=g</IT><SUB>Kr</SUB><IT> · X</IT><SUB>r</SUB><IT> · </IT>(1<IT>/</IT>{1<IT>+</IT>exp[(<IT>V</IT><SUB>m</SUB><IT>+</IT>26)<IT>/</IT>23]})](/content/vol282/issue6/fulltext/H2296/img012.gif)
|
(10)
|
Using the same approach as for INa
but with Te set to 300 ms, we combined
IKr and IKs into a single
IK with a single activation variable
X and reversal potential EK for given
gKr and gKs
|
(11)
|
For 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
APPENDIX as Eqs. 34-36. The value of
gK 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 APPENDIX). This minor change had no marked influence on the shape of the AP or on
the time course of IK (data not shown). Note
that using a single IK 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 gKr and gKs
(see RESULTS).
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 2.
A and B: steady-state value
X (A) and time constant
X (B) of the activation variable
X of the delayed rectifier K+ current
(IK) in the reformulated PB model. C
and D: comparison of IK in the
reduced and the original PB model. Activation of
IK in response to a voltage-clamp step from a
holding potential of 90 mV to a series of test potentials of 20,
10, 0, +10, +20, and + 30 mV is shown. C:
reformulated model; D: original PB model.
|
|
Figure 2, C and D, shows simulated
voltage-clamp records of the original (sum of
IKr and IKs) and
reformulated IK. In the reformulated model,
IK 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
IK1, we used the original equations from the PB
model, which appear in the APPENDIX as Eqs.
43-50. For IK1, we used the original
Eqs. 39-42 (see the APPENDIX) but increased
its conductance (gK1) 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
|
(12)
|
where Cm is the cell capacitance per
surface unit, S is the surface-to-volume ratio, and 
is
the cellular resistivity. Iion is the net
membrane current, which consists of the aforementioned INa, ICa,
Ito, IK, and
IK1 as well as background Na+ and
Ca2+ currents (INa,b and
ICa,b, respectively) and the net currents generated by the electrogenic Na+-K+ pump and
Na+/Ca2+ exchanger (INaK
and INaCa, respectively). For
Cm = 2.0 µF/cm2 and
S = 0.2 µm1 (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, tmin = 0.02 ms and
tmax = 0.1 ms. In this scheme, for each
point of the medium, the ODE is solved with
tmax; if
dVm/dt 
1 V/s, we move on to the
next point; otherwise, we integrate the ODE at that point five times with tmin 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). Table
2 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 Vm 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. 3A and Table 3). We observed
an almost perfect match between INa in reduced
and original models (Fig. 3B). Peak
ICa increased by 40%, but the global time
course was preserved (Fig. 3C). The time course and
amplitude of IK were similar in the reduced and original models (Fig. 3D). In the reduced model,
Ito shows a significantly larger peak and an
earlier inactivation compared with the original model (Fig.
3E). 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
Ito (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).
View larger version (29K):
[in this window]
[in a new window]
|
Fig. 3.
Action potential (AP) and associated membrane currents in
the reformulated PB model and in the original PB model. The AP
shown is the fifth from a cell paced at 1 Hz with a 2-ms stimulus
current of two times threshold. A: membrane potential
(Vm); B: INa;
C: L-type Ca2+ current
(ICa); D: IK;
E: transient outward K+ current
(Ito).
|
|
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.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 4.
Restitution of the AP duration (APD). A: APD
restitution curves in the reformulated model, original PB model, and
experimental data obtained by Morgan et al. (29).
B and C: dependence of AP shape on diastolic
interval (DI). DI was gradually decreased from 500 to 50 ms, as
indicated. B: reformulated model; C: original PB
model.
|
|
Restitution of CV
The activation properties of the original PB model and our reduced
model are compared in Fig. 5. Figure
5A shows the stimulus strength-duration curves of both
models, whereas Fig. 5B shows the CV restitution curves of
the PB model. The decrease in CV with decreasing DI is due to the
decreasing amount of INa 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(Vm).
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 5.
Activation properties of the reformulated model and
original PB model. A: stimulus strength-duration curves;
B: conduction velocity restitution curves.
|
|
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 IKs (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 Ito.
From epicardium to endocardium, Ito shows a
decrease in amplitude, a negative shift in half-inactivation voltage
(Vshift), 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
Ito and IK (through its
slow component IKs). To create an endocardial
cell model, we only made changes to Ito.
The amplitude of Ito was changed by altering
gto. We increased the inactivation time of
Ito by multiplying the gating functions 
to(Vm) and
to(Vm) 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 Vshift. The resulting
equations for to(Vm) and
to(Vm) are listed in the
APPENDIX as Eqs. 31 and 32. To
prolong the APD of the M cells, we decreased gKs
by 50% in the original PB model and performed new voltage-clamp
simulations on the sum of IKr and
IKs to obtain X,
X, and gK for the M
cell. The thus-obtained tabulated values of X
and X are well approximated by the functions
listed in the APPENDIX as Eqs. 37 and 38. The new gK 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 6D 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.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 6.
AP
and APD restitution in 3 different configurations of the reformulated
model, representing the different cell types that are present in the
ventricular wall. A-C: AP of an epicardial cell
(A), a midmyocardial cell (M cell) (B), and an
endocardial cell (C). D: APD restitution curves
for each of the cell types.
|
|
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. 7A), the spiral tip drifts along the refractory tail of the wave (Fig. 7B) before curling back
(clockwise) and stabilizing approximately at the middle of the medium
(Fig. 7C). Figure 7D 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).
View larger version (74K):
[in this window]
[in a new window]
|
Fig. 7.
Spiral wave activity in a 2-dimensional (2-D) representation of the
reformulated model. A-C: formation of spiral wave by an
S1-S2 protocol. D: trajectory of spiral tip during 2 s
after stabilization. See text for details.
|
|
Figure 8A shows a record of
membrane potential in a point far from the spiral core (asterisk in
Fig. 7C) 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). Figure
8B 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).
View larger version (42K):
[in this window]
[in a new window]
|
Fig. 8.
Membrane potential and electrocardiogram (ECG) during
spiral wave activity shown in Fig. 7. Note the different time scales.
A: Vm record at point (650, 650) during 10 s after S1; B: nonscaled ECG
(time 0 at S1 delivery).
|
|
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.
9B), whereas the mean
[K+]i decreased nonmonotonically from 140 to
139.5 mM (Fig. 9C). 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. 9A).
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 9.
Beat-to-beat changes in APD (A), mean
intracellular Na+ concentration
([Na+]i; B), and mean
intracellular K+ concentration
([K+]i; C) in the original PB
model paced at 1 Hz.
|
|
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 [Ca2+]i is
observed (Fig. 10B) and it
takes >1 min of model time for the cell to reach a new steady state
(Fig. 10A).
View larger version (22K):
[in this window]
[in a new window]
|
Fig. 10.
Beat-to-beat changes in APD (A) and mean
intracellular Ca2+ concentration
([Ca2+]i; B) in the original PB
model ([Na+]i and
[K+]i fixed at 10 and 140 mM, respectively)
upon a step increase in stimulation frequency from 1 to 2 Hz at AP 100. C: beat-to-beat changes in APD in the reformulated model
using the same stimulation protocol.
|
|
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. 10C), 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 gK1 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 IK1 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
gK1, the maximum of IK1
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 Ito and ICa are, however,
larger than in the PB model, due to the adiabatic elimination of the
activation variables r and d, respectively. The
amplitude of Ito is furthermore enhanced by a
33% increase in gto. 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.
4A). 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., gCa)
(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.
5B). 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
Ito and IKs (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 IKs 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,
IK1 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 Ito and
IKs, 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, gto and gK1. 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
Vm and the gating variables of
INa (m and v),
ICa (f), Ito
(to), and IK (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 Ca2+ 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
Ca2+ channels and Ca2+-induced Ca2+
release from ryanodine receptor channels exhibiting adaptation. Whereas
this model reproduces important aspects of cardiac Ca2+
cycling, it also exhibits all-or-none rather than graded
Ca2+ 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.
TOP
ABSTRACT
INTRODUCTION
![<IT>f</IT><SUB>Ca</SUB>([Ca<SUP>2+</SUP>]<SUB><IT>i</IT></SUB>)<IT>=</IT>1<IT>/</IT>(1<IT>+</IT>[Ca<SUP>2+</SUP>]<SUB>i</SUB><IT>/</IT>0.0006)](/content/vol282/issue6/fulltext/H2296/img008.gif)
![E<SUB>Na</SUB><IT>=</IT>(<IT>R</IT>T<IT>/F</IT>) ln ([Na<SUP>+</SUP>]<SUB>e</SUB><IT>/</IT>[Na<SUP>+</SUP>]<SUB>i</SUB>)](/content/vol282/issue6/fulltext/H2296/img017.gif)
![E<SUB>Ca</SUB><IT>=</IT>(<IT>R</IT>T<IT>/</IT>2<IT>F</IT>) ln ([Ca<SUP>2+</SUP>]<SUB>e</SUB><IT>/</IT>[Ca<SUP>2+</SUP>]<SUB>i</SUB>)](/content/vol282/issue6/fulltext/H2296/img018.gif)
![E<SUB>to</SUB><IT>=</IT>(<IT>R</IT>T<IT>/F</IT>) ln [(0.043<IT> · </IT>[Na<SUP>+</SUP>]<SUB>e</SUB>](/content/vol282/issue6/fulltext/H2296/img019.gif)
![<IT>+</IT>[K<SUP>+</SUP>]<SUB><IT>e</IT></SUB>)<IT>/</IT>(0.043<IT> · </IT>[Na<SUP>+</SUP>]<SUB>i</SUB><IT>+</IT>[K<SUP>+</SUP>]<SUB>i</SUB>)]](/content/vol282/issue6/fulltext/H2296/img020.gif)
![E<SUB>K</SUB><IT>=</IT>(<IT>R</IT>T<IT>/F</IT>) ln ([K<SUP>+</SUP>]<SUB>e</SUB><IT>/</IT>[K<SUP>+</SUP>]<SUB>i</SUB>)](/content/vol282/issue6/fulltext/H2296/img021.gif)
![&agr;<SUB>m</SUB>=[0.32 · (V<SUB>m</SUB><IT>+</IT>47.13)]<IT>/</IT>{1<IT>−</IT>exp[−0.1<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>47.13)]}](/content/vol282/issue6/fulltext/H2296/img023.gif)
![v<SUB>∞</SUB>=0.5 · [1−tanh(7.74<IT>+</IT>0.12<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img025.gif)
![&tgr;<SUB>v</SUB>=0.25+2.24 · ([1−tanh(7.74<IT>+</IT>0.12<IT> · V</IT><SUB>m</SUB>)]<IT>/</IT>](/content/vol282/issue6/fulltext/H2296/img026.gif)
![{1<IT>−</IT>tanh[0.07<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>92.4)]})](/content/vol282/issue6/fulltext/H2296/img027.gif)
![<IT>×</IT>[(<IT>V</IT><SUB>m</SUB><IT>−</IT>22.36)<IT>/</IT>16.68]<SUP>2</SUP>})<IT>/</IT><FENCE>16.68<IT> · </IT><RAD><RCD>(2<IT>&pgr;</IT>)</RCD></RAD></FENCE>](/content/vol282/issue6/fulltext/H2296/img031.gif)
![<IT>×</IT>[(<IT>V</IT><SUB>m</SUB><IT>−</IT>6.27)<IT>/</IT>14.93]<SUP>2</SUP>})<IT>/</IT><FENCE>14.93<IT> · </IT><RAD><RCD>(2<IT>&pgr;</IT>)</RCD></RAD></FENCE>](/content/vol282/issue6/fulltext/H2296/img033.gif)
![&agr;<SUB>f</SUB>=(6.87 · 10<SUP>−3</SUP>)<IT>/</IT>{1<IT>+</IT>exp[−(6.1546<IT>−V</IT><SUB>m</SUB>)<IT>/</IT>6.12]}](/content/vol282/issue6/fulltext/H2296/img034.gif)
![&bgr;<SUB>f</SUB>={0.069 · exp[−0.11<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>9.825)]<IT>+</IT>0.011}<IT>/</IT>](/content/vol282/issue6/fulltext/H2296/img035.gif)
![{1<IT>+</IT>exp[−0.278<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>9.825)]}<IT>+</IT>5.75<IT> · </IT>10<SUP>−4</SUP>](/content/vol282/issue6/fulltext/H2296/img036.gif)
![f<SUB>Ca</SUB><IT>=</IT>1<IT>/</IT>(1<IT>+</IT>[Ca<SUP>2<IT>+</IT></SUP>]<SUB>i</SUB><IT>/</IT>0.0006)](/content/vol282/issue6/fulltext/H2296/img037.gif)
![&agr;<SUB>r</SUB>={0.5266 · exp[−0.0166<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>−</IT>42.2912)]}<IT>/</IT>](/content/vol282/issue6/fulltext/H2296/img040.gif)
![{1<IT>+</IT>exp[−0.0943<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>−</IT>42.2912)]}](/content/vol282/issue6/fulltext/H2296/img041.gif)
![<IT>×</IT>(<IT>V</IT><SUB>m</SUB><IT>−</IT>5.0027)]}<IT>/</IT>{1<IT>+</IT>exp[−0.1348](/content/vol282/issue6/fulltext/H2296/img043.gif)
![<IT>×</IT>(<IT>V</IT><SUB>m</SUB><IT>−</IT>5.186<IT> · </IT>10<SUP>−5</SUP>)]}](/content/vol282/issue6/fulltext/H2296/img044.gif)
