Cardiac action potential models
Modelling of cardiac action potentials with systems of nonlinear ordinary differential equations (ODEs) (in which the voltage Eq. (2) is coupled to differential equations that describe the dynamics of the ion channel currents I
ion
) dates back to the work of Denis Noble [21] on both Purkinje fibre and pace-maker cells. Modern cardiac AP models for animal [22], human adult [23] and human induced pluripotent stem cell derived [24] cardiomyocytes comprise dozens of state variables and hundreds of model parameters. Typically, the only explicit time-dependence of cardiac AP models is due to the stimulating current I
sti
such that without stimulation, i.e., I
sti
≡0, the AP model forms an autonomous ODE system. Depending on the actual equations and the particular choice of the model parameters, the autonomous ODE system may show limit cycle behaviour (then modelling spontaneous beating activity such as in pace-maker cells) or steady state behaviour (then modelling cells that show no spontaneous AP activity such as ventricular cardiomyocytes). Using inverse bifurcation analysis with sparsity promoting penalization [25, 26], models for spontaneously beating cells can be transferred into models for non-active cells and vice versa. Anyhow, the application of a periodic stimulus I
sti
with pacing cycle length PCL leads to AP models of either periodically paced and spontaneously active or periodically paced and spontaneously non-active cells. In our study we used the following representatives for each of the three available classes of AP models. If the values of the model parameters are not explicitly mentioned in the text, they were chosen to be identical to those used in the original publications.
AP Model PP –Periodically paced pacemaker cell
As an example of an AP model of periodically paced and spontaneaously active cells we chose
$$\begin{array}{@{}rcl@{}} C \frac{dV}{dt} &=& - G_{Ca} d_{\infty}(V) f (V-E_{Ca}) - G_{K} x (V-E_{K})\\ &&+\, I_{sti}(t,PCL), \\ \frac{df}{dt} &=& \frac{f_{\infty}(V)-f}{\tau_{f}}, \\ \frac{dx}{dt} &=& \frac{x_{\infty}(V)-x}{\tau_{x}}, \end{array} $$
(4)
as introduced in [12] and subsequently also used in [16, 27]. In particular, this ODE model of state dimension n=3 includes the inward calcium current
$$I_{Ca} = G_{Ca} d_{\infty}(V) f (V-E_{Ca}) $$
with the calcium channel conductance G
Ca
and the dynamic inactivation variable f, as well as the outward potassium current
$$I_{K} = G_{K} x (V-E_{K}) $$
with the potassium channel conductance G
K
and the dynamic activation variable x. In our discussion we used this model with the original setting τ
f
=80 and τ
x
=300 for the relaxation variables [12] and the initial conditions y
0=(V
0,f
0,x
0)=(−79.46,0.9989,0.1153).
AP Model PV – Periodically paced ventricular cell
As an example of an AP model of periodically paced and spontaneously non-active cells we chose
$$\begin{array}{@{}rcl@{}} C \frac{dV}{dt} &=& - I_{Ca}(V) - I_{K}(V) - I_{Na}(V) - I_{0}(V)\\&& +\, I_{sti}(t,PCL), \\ \frac{dd}{dt} &=& \frac{d_{\infty}(V)-d}{\alpha \tau_{d}(V)}, \; \frac{df}{dt} = \frac{f_{\infty}(V)-f}{\beta \tau_{f}(V)}, \\ \frac{dx}{dt} &=& \frac{x_{\infty}(V)-x}{\gamma \tau_{x}(V)},\\ \frac{dh}{dt} &=& \frac{h_{\infty}(V)-h}{ \tau_{h}(V)}, \; \frac{dm}{dt} = \frac{m_{\infty}(V)-m}{ \tau_{m}(V)}, \\ \frac{dj}{dt} &=& \frac{j_{\infty}(V)-j}{ \tau_{j}(V)}. \end{array} $$
(5)
This model was previously used in [10] for the study of chaotic EADs and is a slightly modified version of the LR91-model [11] for ventricular cardiomyocytes in which the intracellular calcium now is set constant. Here, the currents depending on the dynamic gating variables d, f, x, h, m and j are
$$\begin{array}{@{}rcl@{}} I_{Ca} &=& G_{Ca} d f (V-E_{Ca}) \\ I_{K} &=& G_{K} x \bar{x}(V) (V-E_{K}) \\ I_{Na} &=& G_{Na} m^{3} h j (V-E_{Na}), \end{array} $$
while I
0(V) summarizes those currents that are not directly driven by the dynamic gating variables. In our study we chose the relaxation parameters α=1, β=1 and γ=2.5 (as in Figure 4 of [10]) and the initial conditions y
0=(V
0,d
0,f
0,x
0,h
0,m
0,j
0)
$${} {\small{\begin{aligned} y_{0} = (-84.5286,3\cdot 10^{-6},1,0.3158,0.9832,0.0017,0.995484) \end{aligned}}} $$
(6)
as also considered in [10].
AP Model UP – Unpaced pacemaker cell
As an example of an AP model of unforced and spontaneously active cells we chose Eq. (4) with the setting τ
f
=18, τ
x
=100 (from [27], see also [16]) and I
sti
≡0, and the initial conditions y
0=(V
0,f
0,x
0)=(−75.16,0.9984,0.03976).
Numerical simulation of action potential models
For the numerical simulation of the AP models PP (paced pacemaker), PV (paced ventricular) and UP (unpaced pacemaker) we used the MATLAB [28] solver ode15s for stiff ODE systems with the relative tolerance option set to 10−12. All models were simulated for a time span of 2000 s. Then, the first 500 s were discarded in order to eliminate possible transient effects, such that only the remaining 1500 s of the simulations were used for the analysis.
Numerical bifurcation analysis of AP models
Due to the nonlinearities involved, cardiac AP models Eq. (2) may show a variety of complex dynamical behaviour. Bifurcation analysis [29, 30] is the study of the dynamical repertoire and its dependence on the model parameters. Bifurcations are qualitative changes in the dynamics as model parameters are varied, and the parameter values at which they occur are referred to as bifurcation points. In the bifurcation analysis of the AP models PP, PV and UP we used the potassium channel conductance G
K
as the primary bifurcation parameter and, for models PP and PV, the pacing cycle length PCL as the secondary one. The bifurcation diagrams of the AP models PP, PV and UP were obtained by means of the numerical continuation packages matcont [31] and auto-07p [32]. Continuation is the technique of following a particular solution (such as a fixed point or a limit cycle) of an autonomous system as the continuation parameter changes. One advantage of continuation is that both stable and unstable solution branches can be calculated while bifurcations of fixed points and limit cycles can be simultaneously detected. For the application of this technique to the non-autonomous models PP and PV, the latter were transformed into autonomous models of increased state dimension. To this end we introduced the two-dimensional dynamical system
$$\begin{array}{@{}rcl@{}} \frac{du_{1}}{dt} & = & u_{1}\left(1-u_{1}^{2}-u_{2}^{2}\right) - \frac{2 \pi}{PCL} u_{2}, \\ \frac{du_{2}}{dt} & = & u_{2}\left(1 - u_{1}^{2} - u_{2}^{2}\right) + \frac{2 \pi}{PCL} u_{1}, \end{array} $$
(7)
that admits an asymptotically stable periodic orbit [33], and appended it to Eq. (2) with the setting
$${} {\small{\begin{aligned} I_{sti} = \frac{A}{1+\text{exp} \left[ 5 \cdot 10^{6} \left\{ (1-u_{1}) \cos \left(\frac{d \pi}{PCL} \right) - u_{2} \sin \left(\frac{d \pi}{PCL} \right) \right\} \right] }. \end{aligned}}} $$
(8)
Here, A and d denote the amplitude and the duration of the step pulse, see the Additional file 1 for further details.
The stable branches of the bifurcation diagrams were cross-checked and partially complemented (in case of a numerical continuation failure) by a parametric sweep, i.e., by numerical simulations of Eq. (2) that were repeated for a variety of different parameter values.
Multiple time scale analysis of AP models
Often, cardiac AP models Eq. (2) have state variables with time derivatives of much smaller magnitude than those of other variables. Then, multiple time scale analysis [34] can be used to study the flow of the full system Eq. (2) by splitting the variables into slow and fast ones and by analyzing the corresponding slow and fast subsystems of Eq. (2). One underlying rationale is that the trajectory of the full system under certain conditions evolves along the bifurcation scaffold built by the fast subsystem in which the slow variables serve as bifurcation parameters. In particular, the multiple time scales approach was followed in [10] to discuss the genesis of EADs. More precisely, the authors first eliminated the sodium current I
Na
from Eq. (5) (as of minor relevance for the AP repolarization phase) in order to arrive at an AP model with state dimension n=4. Then, considering the potassium gating variable x as much slower (and also eliminating the stimulating current I
sti
), the fast subsystem
$${} {{\begin{aligned} C \frac{dV}{dt} &= - G_{Ca} d f (V-E_{Ca}) - G_{K} x \bar{x}(V) (V-E_{K}) - I_{0}(V), \\ \frac{dd}{dt} &= \frac{d_{\infty}(V)-d}{\alpha \tau_{d}(V)}, \\ \frac{df}{dt} &= \frac{f_{\infty}(V)-f}{\beta \tau_{f}(V)} \end{aligned}}} $$
(9)
was derived in which x acts as a model parameter. Finally, a bifurcation analysis of Eq. (9) for the setting α=0.1 and β=1.1 and with x as the continuation parameter revealed that Eq. (9) features a supercritical Hopf bifurcation followed by a homoclinic bifurcation to a saddle fixed point. While the Hopf bifurcation was considered in [10] to be necessary for the genesis of EADs in the full system Eq. (5), the homoclinic bifurcation was introduced as the reason for the occurence of chaotic EAD dynamics in Eq. (5) under an appropriate periodic stimulation.
While, with respect to chaotic EADs, our paper suggests to rather explore the full AP dynamics of Eq. (2) than fast subsystems, we also studied for the sake of comparison the fast subsystem
$$\begin{array}{@{}rcl@{}} C \frac{dV}{dt} &=& - G_{Ca} d_{\infty}(V) f (V-E_{Ca}) - G_{K} x (V-E_{K}),\\ \frac{df}{dt} &=& \frac{f_{\infty}(V)-f}{\tau_{f}} \end{array} $$
(10)
of the cardiac AP models PP and UP, then again with the slow variable x as the bifurcation parameter. The submodel Eq. (10) was previously used in [16] in order to demonstrate that a supercritical Hopf bifurcation in a fast subsystem is actually not necessary for EAD genesis.
Calculation of Lyapunov exponents
Cardiac AP models Eq. (2) may feature both regular and chaotic EAD dynamics. In order to test whether a given trajectory is actually chaotic, we calculated the largest Lyapunov exponent λ from the simulated voltage time series data using the TISEAN package [35]. The central idea is to consider a trajectory on the presumably chaotic attractor and to pick a point \(\mathbf {s}_{n_{1}}\) in the state space. Then, a second point \(\phantom {\dot {i}\!}\mathbf {s}_{n_{2}}\) in close distance \(\delta _{0} = \|\mathbf {s}_{n_{1}}-\mathbf {s}_{n_{2}}\|\phantom {\dot {i}\!}\) is chosen, and the separation of the two trajectories over time is measured in terms of
$$\| \mathbf{s}_{n_{1} + \Delta n}-\mathbf{s}_{n_{2} + \Delta n} \| \approx \delta_{0} e^{\lambda \Delta n}. $$
A positive value of λ corresponds with an exponentially fast growth of the initial perturbation δ
0 and means that the respective trajectory is chaotic.
Calculation of restitution curves
As considered in [5] and [12] as an alternative to the S1-S2-protocol, we derived the restitution curves directly from the simulated voltage traces. To this end, all time points at which the voltage crosses the threshold value V
th
were recorded during the simulations using the MATLAB event option of the ODE solver. Action potential durations APDs were then calculated as the time spans during which the voltage lies above V
th
, while diastolic intervals DIs were constructed as the time spans with voltage values below V
th
. Finally, pairs of DI and subsequent APD were built and plotted as APD
n+1 vs. DI
n
. Our motivation for the use of this direct method was its straightforward applicability to both paced and spontaneaously active models without the need of introducing perturbations to the latter. For the sake of completeness, we also applied the S1-S2-protocol, see the Additional file 2.
