Vol. 282, Issue 6, H2296-H2308, June 2002
A computationally efficient electrophysiological model of
human ventricular cells
O.
Bernus1,
R.
Wilders2,
C. W.
Zemlin3,
H.
Verschelde1, and
A. V.
Panfilov4
1 Department of Mathematical Physics and Astronomy,
Gent University, 9000 Gent, Belgium;
2 Department of Physiology, Academic Medical Center,
University of Amsterdam, 1105 AZ Amsterdam; and Department of Medical
Physiology, University Medical Center Utrecht, 3584 CG Utrecht;
3 Institute for Theoretical Biology, Humboldt
University, 10115 Berlin, Germany; and
4 Department of Theoretical Biology, Utrecht University,
3584 CG Utrecht, The Netherlands
 |
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
| |
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.
| |
APPENDIX. Model Equations |
Equilibrium Potentials
where R is the universal gas constant, T is the
absolute temperature, F is the Faraday constant, and
[Na+]e, [Ca2+]e,
and [K+]e are the extracellular
Na+, Ca2+, and K+ concentrations, respectively.
Inward Currents
Fast Na+ current.
|
(13)
|
![&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)
|
(14)
|
|
(15)
|
![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)
|
(16)
|
![&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)
|
(17)
|
Slow Ca2+ current.
|
(18)
|
|
(19)
|
|
(20)
|
|
(21)
|
![&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)
|
(22)
|
![&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)
|
(23)
|
![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)
|
(24)
|
Outward Currents
Transient outward current.
|
(25)
|
|
(26)
|
![&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)
|
(27)
|
|
(28)
|
|
(29)
|
|
(30)
|
![&tgr;<SUB>to</SUB>(V<SUB>m</SUB>)<IT>=</IT>1<IT>/</IT>[<IT>p&agr;</IT><SUB>to</SUB>(<IT>V</IT><SUB>m</SUB>)<IT>+p&bgr;</IT><SUB>to</SUB>(<IT>V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img049.gif)
|
(31)
|
|
(32)
|
Delayed rectifier K+ current.
|
(33)
|
For endocardial and epicardial cells
![X<SUB>∞</SUB>=0.988/[1+exp(−0.861<IT>−</IT>0.0620<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img053.gif)
|
(34)
|
![&tgr;<SUB>X</SUB>=240 · exp[−(25.5<IT>+V</IT><SUB>m</SUB>)<SUP>2</SUP><IT>/</IT>156]](/content/vol282/issue6/fulltext/H2296/img054.gif)
|
(35)
|
![&tgr;′<SUB>X</SUB>=40 · [1−tanh(160<IT>+</IT>2<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img056.gif)
|
(36)
|
For M cells
![X<SUB>∞</SUB>=0.972/[1+exp(−2.036<IT>−</IT>0.0834<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img057.gif)
|
(37)
|
![&tgr;<SUB>X</SUB>=380 · exp[−(25.5<IT>+V</IT><SUB>m</SUB>)<SUP>2</SUP><IT>/</IT>156]](/content/vol282/issue6/fulltext/H2296/img058.gif)
|
(38)
|
Inward rectifier K+ current.
|
(39)
|
|
(40)
|
![&agr;<SUB>K1</SUB><IT>=</IT>0.1<IT>/</IT>{1<IT>+</IT>exp[0.06<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>−E</IT><SUB>K</SUB><IT>−</IT>200)]}](/content/vol282/issue6/fulltext/H2296/img062.gif)
|
(41)
|
![&bgr;<SUB>K1</SUB><IT>=</IT>{3<IT> · </IT>exp[2<IT> · </IT>10<SUP><IT>−</IT>4</SUP><IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>−E</IT><SUB>K</SUB><IT>+</IT>100)]<IT>+</IT>exp[0.1](/content/vol282/issue6/fulltext/H2296/img063.gif)
|
(42)
|
Background Currents
Ca2+ background current.
|
(43)
|
Na+ background current.
|
(44)
|
Pump and Exchanger Currents
Na+-K+
pump.
|
(45)
|
|
(46)
|
![f′<SUB>NaK</SUB><IT>=</IT>{1<IT>/</IT>[1<IT>+</IT>(10<IT>/</IT>[Na<SUP>+</SUP>]<SUB>i</SUB>)<SUP>1.5</SUP>]}<IT> · </IT>[[K<SUP>+</SUP>]<SUB>e</SUB><IT>/</IT>([K<SUP>+</SUP>]<SUB>e</SUB><IT>+</IT>1.5)]](/content/vol282/issue6/fulltext/H2296/img070.gif)
|
(47)
|
![&sfgr;=0.1428 · [exp([Na<SUP>+</SUP>]<SUB>e</SUB><IT>/</IT>67.3)<IT>−</IT>1]](/content/vol282/issue6/fulltext/H2296/img071.gif)
|
(48)
|
Na+/Ca2+
exchanger.
|
(49)
|
![<IT>×</IT>exp(−0.024<IT> · V</IT><SUB>m</SUB>)}<SUP>−1</SUP><IT> · </IT>{[Na<SUP>+</SUP>]<SUP>3</SUP><SUB>i</SUB><IT> · </IT>[Ca<SUP>2+</SUP>]<SUB>e</SUB>](/content/vol282/issue6/fulltext/H2296/img074.gif)
|
(50)
|
| |
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.
| |
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
Received 14 August 2001; accepted in final form 21 January 2002.
| |
REFERENCES |
1.
Antzelevitch, C,
Nesterenko VV,
Muzikant AL,
Rice JJ,
Chen G,
and
Colatsky T.
Influence of transmural repolarization gradients on the electrophysiology and pharmacology of ventricular myocardium. Cellular basis for the Brugada and long QT-syndromes.
Phil Trans R Soc London A
359:
1201-1216,
2001.
2.
Arce, H,
Xu A,
González H,
and
Guevara MR.
Alternans and higher-order rhythms in an ionic model of a sheet of ischemic ventricular muscle.
Chaos
10:
411-426,
2000[Web of Science][Medline].
3.
Ashamalla, SM,
Navarro D,
and
Ward CA.
Gradient of sodium current across the left ventricular wall of adult rat hearts.
J Physiol (Lond)
536:
439-443,
2001.
4.
Beuckelmann, DJ,
Näbauer M,
and
Erdmann E.
Intracellular calcium handling in isolated ventricular myocytes from patients with terminal heart failure.
Circulation
85:
1046-1055,
1992.
5.
Beuckelmann, DJ,
Näbauer M,
and
Erdmann E.
Alterations of K+ currents in isolated human ventricular myocytes from patients with terminal heart failure.
Circ Res
73:
379-385,
1993[Abstract/Free Full Text].
6.
Cherry, EM,
Greenside HS,
and
Henriquez CS.
A space-time adaptive method for simulating complex cardiac dynamics.
Phys Rev Lett
84:
1343-1346,
2000.
7.
Clayton, RH,
Murray A,
and
Campbell RW.
Objective features of the surface electrocardiogram during ventricular tachyarrhythmias.
Eur Heart J
16:
1115-1119,
1995[Abstract/Free Full Text].
8.
Courtemanche, M,
Ramirez RJ,
and
Nattel S.
Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model.
Am J Physiol Heart Circ Physiol
275:
H301-H321,
1998[Abstract/Free Full Text].
9.
Davidenko, JM,
Pertsov AM,
Salomonsz R,
Baxter W,
and
Jalife J.
Stationary and drifting spiral waves of excitation in isolated cardiac muscle.
Nature
355:
349-351,
1991.
10.
Efimov, IR,
Krinsky VI,
and
Jalife J.
Dynamics of rotating vortices in the Beeler-Reuter model of cardiac tissue.
Chaos, Solitons and Fractals
5:
513-526,
1995.
11.
Efimov, IR,
Sidorov VY,
Cheng Y,
and
Wollenzier B.
Evidence of 3D scroll waves with ribbon-shaped filament as a mechanism of ventricular tachycardia in the isolated rabbit heart.
J Cardiovasc Electrophysiol
10:
1452-1462,
1999[Web of Science][Medline].
12.
Endresen, LP,
and
Skarland N.
Limit cycle oscillations in pacemaker cells.
IEEE Trans Biomed Eng
47:
1134-1137,
2000[Web of Science][Medline].
13.
Fast, VG,
Efimov IR,
and
Krinsky VI.
Transition from circular to linear cores in excitable media.
Phys Lett A
151:
157-161,
1990.
14.
Fenton, F,
and
Karma A.
Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation.
Chaos
8:
20-47,
1998[Web of Science][Medline].
15.
Franz, MR,
Schaeffer J,
Schottler M,
Seed WA,
and
Noble MIM
Electrical and mechanical restitution of the human heart at different rates of stimulation.
Circ Res
53:
815-822,
1983[Abstract/Free Full Text].
16.
Frazier, DW,
Wolf PD,
Wharton JM,
Tang ASL,
Smith WM,
and
Ideker RE.
Stimulus-induced critical point. Mechanism for electrical induction of reentry in normal canine myocardium.
J Clin Invest
83:
1039-1052,
1989[Web of Science][Medline].
17.
Glass, L.
Dynamics of cardiac arrhythmias.
Phys Today
49:
40-45,
1996.
18.
Gray, RA,
Pertsov AM,
and
Jalife J.
Spatial and temporal organization during cardiac fibrillation.
Nature
392:
75-78,
1998[Medline].
19.
Hodgkin, AL,
and
Huxley AF.
A quantitative description of membrane current and its application to conduction and excitation in nerve.
J Physiol (Lond)
117:
500-544,
1952.
20.
Jafri, MS,
Rice JJ,
and
Winslow RL.
Cardiac Ca2+ dynamics: the roles of rayodine receptor adaptation and sarcoplasmic reticulum load.
Biophys J
74:
1149-1168,
1998.
21.
Jalife, J,
Gray RA,
Morley GE,
and
Davidenko JM.
Self-organization and the dynamical nature of ventricular fibrillation.
Chaos
8:
79-93,
1998[Web of Science][Medline].
22.
Jongsma, HJ,
and
Wilders R.
Gap junctions in cardiovascular disease.
Circ Res
86:
1193-1197,
2000[Abstract/Free Full Text].
23.
Koumi, S,
Backer CL,
and
Arentzen CE.
Characterization of inwardly rectifying K+ channel in human cardiac myocytes. Alterations in channel behavior in myocytes isolated from patients with idiopathic dilated cardiomyopathy.
Circulation
92:
164-174,
1995.
24.
Li, GR,
Feng J,
Yue L,
and
Carrier M.
Transmural heterogeneity of action potentials and Ito1 in myocytes isolated from the human right ventricle.
Am J Physiol Heart Circ Physiol
275:
H369-H377,
1998[Abstract/Free Full Text].
25.
Liu, DW,
Gintant GA,
and
Antzelevitch C.
Ionic bases for electrophysiological distinctions among epicardial, midmyocardial, and endocardial myocytes from the free wall of the canine left ventricle.
Circ Res
72:
671-687,
1993[Abstract/Free Full Text].
26.
Luo, CH,
and
Rudy Y.
A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes.
Circ Res
74:
1071-1096,
1994[Abstract/Free Full Text].
27.
Main, MC,
Bryant SM,
and
Hart G.
Regional differences in action potential characteristics and membrane currents of guinea-pig left ventricular myocytes.
Exp Physiol
83:
747-761,
1998[Abstract].
28.
Morgan, JM,
Cunningham D,
and
Rowland E.
Dispersion of monophasic action potential duration: demonstrable in humans after premature ventricular extrastimulation but not in steady state.
J Am Coll Cardiol
19:
1244-1253,
1992[Abstract].
29.
Morgan, JM,
Cunningham D,
and
Rowland E.
Electrical restitution in the endocardium of the intact human right ventricle.
Br Heart J
67:
42-46,
1992[Abstract/Free Full Text].
30.
Noble, D.
A modification of the Hodgkin-Huxley equations applicable to Purkinje fibres action and pacemaker potential.
J Physiol (Lond)
160:
317-352,
1962.
31.
Noble, D,
Varghese A,
Kohl P,
and
Noble P.
Improved guinea-pig ventricular cell model incorporating a diadic space, IKr and IKs, and length- and tension-dependent processes.
Can J Cardiol
14:
123-134,
1998[Web of Science][Medline].
32.
Nygren, A,
Fiset C,
Firek L,
Clark JW,
Lindblad DS,
Clark RB,
and
Giles WR.
Mathematical model of an adult human atrial cell. The role of K+ currents in repolarization.
Circ Res
82:
63-81,
1998[Abstract/Free Full Text].
33.
Plonsey, R,
and
Barr RC.
Biolectricity: A Quantitative Approach. New York: Plenum, 1988, p. 213-216.
34.
Press, WH,
Teukolsky SA,
Vetterling WT,
and
Flannery BP.
Numerical Recipes in C: The Art of Scientific Computing. New York: Cambridge University Press, 1992, p. 870-ff.
35.
Priebe, L,
and
Beuckelmann DJ.
Simulation study of cellular electric properties in heart failure.
Circ Res
82:
1206-1223,
1998[Abstract/Free Full Text].
36.
Qu, Z,
and
Garfinkel A.
An advanced algorithm for solving partial differential equation in cardiac conduction.
IEEE Trans Biomed Eng
46:
1166-1168,
1999[Web of Science][Medline].
37.
Qu, Z,
Weiss J,
and
Garfinkel A.
Cardiac electrical restitution properties and stability of reentrant spiral waves: a simulation study.
Am J Physiol Heart Circ Physiol
276:
H269-H283,
1999[Abstract/Free Full Text].
38.
Rice, JJ,
Jafri MS,
and
Winslow RL.
Modeling gain and gradedness of Ca2+ release in the functional unit of the cardiac diadic space.
Biophys J
77:
1871-1884,
1999.
39.
Rush, S,
and
Larsen H.
A practical algorithm for solving dynamic membrane equations.
IEEE Trans Biomed Eng
25:
389-392,
1978[Web of Science][Medline].
40.
Saumarez, RC,
Camm AJ,
Panagos A,
Gill JS,
Stewart JT,
de Belder MA,
Simpson IA,
and
McKenna WJ.
Ventricular fibrillation in hypertrophic cardiomyopathy is associated with increased fractionation of paced right ventricular electrograms.
Circulation
86:
467-474,
1992.
41.
Saumarez, RC,
Heald S,
Gill J,
Slade AKB,
de Belder MA,
Walczak F,
Rowland E,
Ward DE,
and
Camm AJ.
Primary ventricular fibrillation is associated with increased paced right ventricular electrogram fractionation.
Circulation
92:
2565-2571,
1995.
42.
Volders, PG,
Sipido KR,
Carmeliet E,
Spätjens RL,
Wellens HJ,
and
Vos MA.
Repolarizing K+-currents Ito1 and IKs are larger in right than left canine ventricular midmyocardium.
Circulation
99:
206-210,
1999.
43.
Witkowski, FX,
Leon LJ,
Penkoske PA,
Giles WR,
Spano ML,
Ditto WL,
and
Winfree AT.
Spatiotemporal evolution of ventricular fibrillation.
Nature
392:
78-82,
1998[Medline].
44.
Zemlin, CW,
and
Panfilov AV.
Spiral waves in excitable media with negative restitution.
Phys Rev E
63:
041912,
2001.
45.
Zygmunt, AC,
Eddlestone GT,
Thomas GP,
Nesterenko VV,
and
Antzelevitch C.
Larger late sodium conductance in M cells contributes to electrical heterogeneity in canine ventricle.
Am J Physiol Heart Circ Physiol
281:
H689-H697,
2001[Abstract/Free Full Text].
46.
Zygmunt, AC,
Goodrow RJ,
and
Antzelevitch C.
INaCa contributes to electrical heterogeneity within the canine ventricle.
Am J Physiol Heart Circ Physiol
278:
H1671-H1678,
2000[Abstract/Free Full Text].
Am J Physiol Heart Circ Physiol 282(6):H2296-H2308
0363-6135/02 $5.00
Copyright © 2002 the American Physiological Society
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)
![{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)
![{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)
![{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)
![<IT>×</IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>34.2531)]}<IT>/</IT>{1<IT>+</IT>exp[−0.1732<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>34.2531)]}](/content/vol282/issue6/fulltext/H2296/img046.gif)
![<IT>×</IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>34.0235)]}<IT>/</IT>{1<IT>+</IT>exp[−0.1604<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>+</IT>34.0235)]}](/content/vol282/issue6/fulltext/H2296/img048.gif)
![<IT>+&bgr;</IT><SUB>to</SUB>(<IT>V</IT><SUB>m</SUB><IT>−V</IT><SUB>shift</SUB>)]](/content/vol282/issue6/fulltext/H2296/img051.gif)
![<IT>+</IT>182<IT> · </IT>[1<IT>+</IT>tanh(0.154<IT>+</IT>0.0116<IT> · V</IT><SUB>m</SUB>)]<IT>+&tgr;</IT><SUP><IT>′</IT></SUP><SUB><IT>X</IT></SUB>](/content/vol282/issue6/fulltext/H2296/img055.gif)
![<IT>+</IT>166<IT> · </IT>[1<IT>+</IT>tanh(0.558<IT>+</IT>0.0169<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img059.gif)
![<IT>×</IT>(<IT>V</IT><SUB>m</SUB><IT>−E</IT><SUB>K1</SUB><IT>−</IT>10)]}<IT>/</IT>{1<IT>+</IT>exp[−0.5<IT> · </IT>(<IT>V</IT><SUB>m</SUB><IT>−E</IT><SUB>K</SUB>)]}](/content/vol282/issue6/fulltext/H2296/img064.gif)
![<IT>+</IT>0.0365<IT> · &sfgr; · </IT>exp(−0.037<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img069.gif)
![f<SUB>NaCa</SUB><IT>=</IT>(87.5<SUP>3</SUP><IT>+</IT>[Na<SUP>+</SUP>]<SUP>3</SUP><SUB>e</SUB>)<SUP>−1</SUP><IT> · </IT>(1.38<IT>+</IT>[Ca<SUP>2+</SUP>]<SUB>e</SUB>)<SUP>−1</SUP><IT> · </IT>{1<IT>+</IT>0.1](/content/vol282/issue6/fulltext/H2296/img073.gif)
![<IT>×</IT>exp(0.013<IT> · V</IT><SUB>m</SUB>)<IT>−</IT>[Na<SUP>+</SUP>]<SUP>3</SUP><SUB>e</SUB><IT> · </IT>[Ca<SUP>2+</SUP>]<SUB>i</SUB><IT> · </IT>exp(−0.024<IT> · V</IT><SUB>m</SUB>)]](/content/vol282/issue6/fulltext/H2296/img075.gif)
