Poole, Matthew

English

Name not available

Simulation of Cardiac Action Potential Propagation Using Hybrid Models
MJ Poole
University of Portsmouth, Portsmouth, UK
Abstract
This paper investigates hybrid computational models as
an approach to efficient simulation of electrical activity in
cardiac tissue. In the hybrid approach, computational re-
quirements are reduced by modelling only a local region
using a detailed model, and embedding this within a sim-
pler model of a larger tissue. We argue that the validity of
this approach is dependent upon the method used to cou-
ple the different component models. We develop a simple
coupling method for creating hybrid models in which the
ionic Luo-Rudy phase 1 model is embedded within a tissue
that uses the 3-variable Fenton-Karma model. The hybrid
models are shown to preserve depolarisation fronts propa-
gating through tissue as well as accurate action potential
behaviour within embedded regions.
1. Introduction
Mathematical models of the cardiac action poten-
tial include complex biophysically-derived ionic models,
simpler generic ionic models, and phenomenologically-
derived models that caricature the action potential and its
propagation. Complex models are computationally de-
manding, often prohibitively so depending on the space
and time scales that need to be simulated. To investigate
many spatio-temporal behaviours (e.g. reentrant arrhyth-
mias), simpler models are often more computationally ef-
fective. These include reduced variable continuous models
or completely discrete models such as cellular automata.
However, these models are often insufficiently rich to rep-
resent many dynamics and pathological states of interest.
Hybrid, or multi-level, models have been proposed as
a way of gaining many of the benefits of detailed ionic
models but with significantly reduced computational costs
[1, 2]. In the hybrid approach, regions (or single cells) of
tissue modelled at a biophysically-detailed level are em-
bedded within a larger tissue that employs a simpler action
potential model. Thus points or regions of specific inter-
est (e.g. tissue affected by ischemia, or an area around
the tip of a spiral wave) can be realistically modelled, but
with computational requirements greatly reduced for the
surrounding tissue.
The validity of the hybrid approach requires that action
potentials propagate from one region to another without
significant changes in velocity, and with correct morphol-
ogy; in particular, the shape (and duration) of action poten-
tials within an ionic model-based embedded cell or region
should not be adversely affected by the different action po-
tential shape of a simpler model in the surrounding space.
Construction of hybrid models requires that computational
elements from the component models are coupled at the
interface between the regions, and it is this coupling that
determines the influence of one model on another.
Where both models are continuous, an obvious approach
to inter-model coupling is to allow element-element dif-
fusion to apply unaltered across the interface. In [1] it is
shown that, using this diffusive method of coupling, the ve-
locity of depolarisation wavefronts can be preserved across
models, but the action potentials within the embedded re-
gion can take an incorrect shape, influenced by those of the
other model.
In this paper, we consider hybrid simulations which use
two well-known action potential models, the ionic Luo-
Rudy phase 1 model [3] and the simplified Fenton-Karma
model [4]. We demonstrate shortcomings of the diffusive
coupling method when used with these component mod-
els. We then present an alternative method in which inter-
model coupling is dependent upon the phases of the action
potential at the interface of the two regions. Our method
is shown to perform well in simulations of propagation in
1D and 2D hybrid models.
2. Methods
Our method of coupling applies to component compu-
tational models which have a regular spatial structure and
share the following property: the state (i.e. variables) of
each computational element i at time t +∆t is dependent
upon the element’s voltage vt
i
and its other variables (e.g.
gating variables) at time t, and on a diffusion term dt
i
. This
diffusion term involves voltages of neighbouring elements
at time t; for example, for a 1D finite difference model,
dt
i
= D(vt
i−1 + v
t
i+1 − 2v
t
i
)/∆x2, where D is a diffusion
coefficient and ∆x the numerical space-step between ele-
ments. Commonly, the model we wish to embed will have
ISSN 0276−6574 633 Computing in Cardiology 2010;37:633−636.
a finer spatial and temporal structure and each element will
contain a greater number of variables; it will be (possibly
vastly) more computationally intensive. The two models
may also be of different mathematical type (e.g. numer-
ical approximation to a continuous model versus cellular
automaton). For simplicity, in this paper we consider two
specific continuous models, and assume that the compu-
tational implementations share both space-step and time-
step. Specifically, we embed elements or regions of tissue
modelled with the 8-variable ionic Luo-Rudy I cell dynam-
ics (LR) within a larger space which utilises the compu-
tationally less demanding 3-variable Fenton-Karma (FK)
model. Figure 1 gives the structure of such hybrid models.
LR elements
FK elements
Figure 1. Structure of 1D and 2D FK/LR hybrid mod-
els. Single elements or regions of the LR model (shaded)
are embedded within a tissue composed of FK elements.
Interface elements are denoted by • (LR) and ◦ (FK).
We divide the original LR model time constants τ d and
τ f in [3] by two in order to give speeded up calcium ki-
netics. We use the parameter values given in [4] for the
FK model to approximate the restitution properties of the
modified LR model. Thus, the models’ (i) action potential
duration versus diastolic interval (DI) behaviour, and (ii)
conduction velocity versus DI behaviour, both match ap-
proximately. This is clearly an important prerequisite for
a hybrid model to demonstrate valid behaviour. Figure 2
shows action potentials of each of the LR and FK models
at a point on a 1D tissue strand.
-70
-30
 10
 0  200  400  600  800
v 
(m
V)
time (ms)
 0  2  4
time (ms)
Figure 2. Left: Action potentials generated by three stim-
uli at one end of a 1D cable for the LR model (solid line)
and the fitted FK model (dotted line). Right: Upstrokes of
LR and FK action potentials.
We now consider how to couple elements at an FK/LR
interface; i.e., what the diffusion terms d should be at in-
terface elements (those that include one or more elements
from the other component model within their neighbour-
hood). The simplest approach, taken in [1, 2], is to use
a diffusion term equivalent to that used for all other non-
boundary elements. This method has been shown to ade-
quately preserve depolarisation fronts passing through the
interface regions. However, as demonstrated in Figure 3,
action potentials within embedded LR region(s) are af-
fected by the different morphology of those in the sur-
rounding FK space. There are two problems. Firstly, the
plateau and recovery phases at LR interface elements tend
to take the shape of an FK action potential, with this effect
only slowly diminishing as we move towards the centre
of an embedded region (experiments suggest that an LR
‘buffer’ region of approx. 5mm is needed for the shape
and duration to closely approximate those of the LR-only
model). Second, an LR interface element fails to depo-
larise correctly to peak voltage because the upstroke of its
neighbouring FK elements are not sufficiently steep (vari-
ations of the effect illustrated occur for left and right inter-
face elements of a larger embedded LR region).
-70
-30
 10
v 
(m
V)
-70
-30
10
0 100 200
v 
(m
V)
time (ms)
0 2 4
time (ms)
Figure 3. Top: Action potential and upstroke (red) for a
single LR element embedded within a 1D FK model using
diffusive coupling, shown with normal LR (dashed) and
FK (dotted) action potentials. Bottom: Action potential at
the centre of a 1.1mm embedded LR region.
The new coupling method aims to address both prob-
lems, but particularly the former since valid and accurate
behaviour within an embedded region is what motivates
the use of hybrid models. The poor LR upstroke at the in-
terface appears not to affect LR elements further away, and
so this issue needs to be resolved mainly if we wish to use
very small (or single element) embedded regions.
For a typical action potential model the diffusion term
d is significant only during and just after depolarisation.
Figure 4 shows, for the LR model, that d returns rapidly to-
wards zero immediately after depolarisation, and remains
insignificantly small for the remainder of the action poten-
tial. At the interface elements of Figure 3, it is a significant
634
-200
-100
 0
 100
 0  100  200
d 
(m
V/
ms
)
time (ms)
 0  2  4
time (ms)
Figure 4. Diffusion term d (solid) shown alongside an
action potential (dashed) for an LR element on a 1D strand.
value of d that leads to an incorrect action potential shape.
In our coupling method, the diffusion term at an in-
terface element is only used at certain phases of its ac-
tion potential; at other times diffusion is ignored (i.e. as-
sumed to be 0). To determine whether the diffusion term
is used, each interface element is assigned a discrete ‘cou-
pling variable’ which at any time is in exactly one of five
possible states, as illustrated in Figure 5. An element’s
coupling variable changes to a new state only when the
element’s voltage v and diffusion term d satisfy the condi-
tion given on the corresponding transition. For example, if
an element is in the ‘Stimulated’ state, and its voltage be-
comes greater than a threshold vth, then its coupling vari-
able changes to the ‘Upstroke’ state.



	

ﬀ



	

6
6



	

?



	

ﬀ
-



	

Early Re-
polarisationRecovery
Plateau & d > dr
‘off’ states ‘on’ states
Upstroke
Resting Stimulated
and ≥ 1 neighbour in an ‘on’ state
d > dth
v ≤ vth and d ≤ 0
v > vth
dv
dt
≤ 0
v ≤ vr
Figure 5. States and transitions for interface elements.
The three states on the right-hand side of Figure 5 are
denoted ‘on’ states—in these states the diffusion term is
used normally. For the ‘off’ states on the left-hand side,
the diffusion term is ignored. Note that for the transi-
tion from ‘Resting’ to ‘Stimulated’, not only does diffu-
sion need to exceed a threshold dth, but at least one neigh-
bouring element must be in an ‘on’ state. This prevents
an element from becoming prematurely stimulated by a
recovering neighbouring element of the other component
model which has a higher voltage (e.g. due to differ-
ences in action potential shape or restitution). Because
of this, all non-interface neighbours of interface elements
also need to include a coupling variable, although the dif-
fusion terms for such elements are always used normally.
For our FK/LR coupling, we use the following parameter
values: vth = −60, vr = −70; for LR elements, dth = 1.5,
dr = −70; for FK elements, dth = 0.01, dr = −2.5.
To improve upon the upstroke behaviour illustrated in
Figure 3, an LR interface element requires that neighbour-
ing FK elements give a better approximation to an LR up-
stroke. To achieve this, we transform each FK element’s
voltage before it is used in the LR element’s diffusion term.
We note the following from Figure 2: (i) the initial phase
of the models’ upstrokes (upto approx. −60mV) is simi-
lar; (ii) after this, the LR upstroke is steeper; (iii) the LR
upstroke takes a more complex form (with peak); (iv) the
height of the LR upstroke is dependent on the diastolic in-
terval. We introduce a simple piecewise linear map φ to
convert FK voltages, by increasing upstroke steepness:
φ(v) =
{
v if v ≤ −60
1.54(v + 60)− 60 if v > −60.
We do not attempt to reproduce the DI-dependent peak
and early-repolarisation of the LR model when converting
an FK voltage. Rather, the current diastolic interval of each
LR interface element is recorded and used to determine an
upper limit for φ-converted FK voltages used in its diffu-
sion term. We use a pre-computed lookup table to store
these limits, which range from 9mV for DI < 36ms, upto
21mV for DI > 150ms. Thus, the upstrokes of neighbour-
ing FK elements appear to an LR interface element quite
similar to those of square waves of variable amplitude.
To convert a neighbouring LR voltage for an FK inter-
face element, we take an even simpler approach. Voltages
upto −60mV are left unaltered, as for φ. Above −60mV
we use a lookup table to trace a generic FK action potential
independently of the LR element’s voltage.
3. Results
Computer simulations were performed for isotropic
1D and 2D tissues, with diffusion coefficient D =
0.1mm2/ms, using a finite difference method with time-
step∆t = 0.01ms and space-step∆x = 0.1mm. Diastolic
intervals were determined using a threshold of −60mV.
For the 1D simulations, we used two cables of length
45mm, one using the LR model, and one using the FK
model with a single LR element embedded at 30mm. The
first 1mm of each cable was stimulated at different rates
to generate action potentials along the strands. Figure 6
compares the action potentials of the LR model with the
hybrid FK/LR model at 30mm for different diastolic inter-
vals. The upper graph shows that the correct LR action po-
tential shape and duration is preserved in the hybrid model
635
across a wide range of DI values (compare with Figure
3). The lower graphs compare upstrokes of the two mod-
els at the same LR element location for each of the three
DI values. The LR upstroke has been approximated well
throughout most of its height, although the peak is delayed
slightly (compare with Figure 3). Different peak voltages
(for different DIs), are also preserved accurately.
-70
-30
 10
 0  100  200  300
v 
(m
V)
time (ms)
36ms 56ms 225ms
-70
-30
 10
 0  2  4
v 
(m
V)
time (ms)
36ms
 0  2  4
time (ms)
56ms
 0  2  4
time (ms)
225ms
Figure 6. Top: Action potentials for a single LR element
embedded within a 1D FK model (red) shown with those
from the LR-only model (dashed), for three different DIs.
Bottom: Upstrokes of these three action potentials.
−70mV
−30mV
 +10mV
-70
-30
 10
 0  200  400  600  800
v 
(m
V)
time (ms)
Figure 7. A 2mm LR region embedded within a 8mm
square FK model. Top: Snapshots (with interval 1ms) of a
wavefront passing through the LR region. Bottom: Action
potentials generated by three stimuli at the left edge of the
2D hybrid model (red) at the point marked by •, and for a
8mm square LR-only model (dashed).
Figure 7 shows a 2D hybrid model which comprises an
8mm square isotropic tissue using the FK model, with a
2mm square LR region embedded at the centre. The tis-
sue was stimulated along the left-hand edge to give pla-
nar waves travelling through the FK and LR regions. The
snapshots show a single depolarisation front moving from
left to right through the tissue. The LR region can easily
be identified due to the differences in membrane potential
at the top of the upstroke and plateau phase (as illustrated
in Figure 2). The snapshots demonstrate that the LR re-
gion has little effect on the FK depolarisation front; the
near vertical wavefront shows only a slight deformation.
The lower graph compares the hybrid model with an 8mm
square LR-only model. Three stimuli were applied at the
left edge of each model to produce three propagating ac-
tion potentials with different DIs. The membrane potential
of the element identified by • in the snapshots is shown
against that of the corresponding element in the LR model.
The correct action potential shape of the LR model is faith-
fully preserved within the embedded LR region.
For 1s of this 2D simulation, the hybrid model, and an
FK-only model, take 16.8% and 13.6% of the runtime of
the LR-only model, respectively. The computational over-
head for coupling is less than 10% of the total runtime.
4. Discussion and conclusions
Hybrid models of cardiac tissue can provide a
computationally-realistic approach to the simulation of ac-
tion potential propagation. We have demonstrated the im-
portance of coupling at the interface for the hybrid ap-
proach to be useful, and developed a coupling method for
tissues employing the LR and FK action potential models.
This method is shown to successfully preserve the shape
and basic properties of the action potential within an em-
bedded LR region. Whilst being specific to these two com-
ponent models, the method appears sufficiently simple to
allow it to be adapted for use with other models.
References
[1] Poole MJ, Holden AV, Tucker JV. Hierarchical reconstruc-
tions of cardiac tissue. Chaos Solitons and Fractals 2002;
13:1581–1612.
[2] Defontaine A, Herna´ndez A, Carrault G. Multi-formalism
modelling and simulation: Application to cardiac modelling.
Acta Biotheoretica 2004;52:273–290.
[3] Luo CH, Rudy Y. A model of the ventricular cardiac ac-
tion potential. Depolarization, repolarization, and their inter-
action. Circulation Research 1991;68:1501–1526.
[4] Fenton F, Karma A. Vortex dynamics in three-dimensional
continuous myocardium with fiber rotation: Filament insta-
bility and fibrillation. Chaos 1998;8:20–47.
Address for correspondence:
Matthew Poole
School of Computing, University of Portsmouth
Lion Terrace, Portsmouth PO1 3HE, UK
matthew.poole@port.ac.uk
636


Simulation of cardiac action potential propagation using hybrid models

Portsmouth University Research Portal (Pure)

A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction.

Hierarchical reconstructions of cardiac tissue. Chaos Solitons and Fractals

Multi-formalism modelling and simulation: Application to cardiac modelling.

Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation. Chaos 1998;8:20–47. Address for correspondence:

https://researchportal.port.ac.uk/portal/services/downloadRegister/82182/0633.pdf

Simulation of cardiac action potential propagation using hybrid models

Abstract

Similar works

Full text

Available Versions

Name not available

Portsmouth University Research Portal (Pure)