Why Pacing Frequency Affects the Production of Early Afterdepolarizations in Cardiomyocytes: An Explanation Revealed by Slow/Fast Analysis of a Minimal Model

TL;DR

This paper uses a minimal mathematical model to explain how pacing frequency influences early afterdepolarizations in heart cells, revealing the underlying canard dynamics near folded node singularities.

Contribution

It demonstrates that EADs are canards near folded node singularities, providing a new analytical framework to predict EAD occurrence without extensive simulations.

Findings

Low-frequency pacing facilitates EAD production.

Intermediate-frequency pacing leads to complex alternans dynamics.

Canard analysis explains the number and timing of EADs.

Abstract

Early afterdepolarizations (EADs) are pathological voltage oscillations in cardiomyocytes that have been observed in response to a number of pharmacological agents and disease conditions. EADs are small voltage fluctuations that occur during the plateau of an action potential. Although a single-cell behavior, EADs can lead to tissue-level arrhythmias, including ventricular tachycardia. Much is currently known about the biophysical mechanisms (i.e., the roles of ion channels and intracellular calcium stores) for EADs, due partially to the development and analysis of mathematical models. This includes the application of slow/fast analysis, which takes advantage of timescale separation inherent in the system to simplify its analysis. We take this further, using a minimal 3D model to demonstrate that the phase-2 EADs are canards that are formed in the neighborhood of a folded node singularity. This knowledge allows us to determine the number of EADs that can be produced for a given parameter set without performing computer simulations, and provides guidance on parameter changes that can facilitate or inhibit EAD production. With this approach, we demonstrate why periodic stimulation, as would occur in an intact heart, preferentially facilitates EAD production when applied at low frequencies,. We also explain the origin of complex alternan dynamics that can occur with intermediate-frequency stimulation, in which varying numbers of EADs are produced with each stimulation. These revelations fall out naturally from an understanding of folded node singularities, but are hard or impossible to glean from a knowledge of the biophysical mechanism for EADs alone. Therefore, an understanding of the canard mechanism is a useful complement to an understanding of the biophysical mechanism that has been developed over years of experimental and computational investigations.

Equations

$$\begin{array}{rl}{C}_m\frac{dV}{dt}& =-\left(I_{\mathrm{K}}+I_{\mathrm{Ca}}\right)+I_{\mathrm{sti}},\\ \frac{dn}{dt}& =\frac{n_{\infty }(V)-n}{{\tau }_n},\\ \frac{dh}{dt}& =\frac{h_{\infty }(V)-h}{{\tau }_h},\end{array}$$

$$I_{\mathrm{K}}=g_{K}n\left(V-V_{\mathrm{K}}\right)\text{ and }{I}_{\mathrm{Ca}}=g_{Ca}{m}_{\infty }(V)h\left(V-V_{\mathrm{Ca}}\right),$$

$$I_{\mathrm{sti}}=40\sum _{k\in \mathbb{N}}\left[H\left(t-k\cdot \mathrm{PCL}\right)-H\left(t-(k\cdot \mathrm{PCL}+1)\right)\right]$$

$$x_{\infty }(V)=\frac{1}{1+\exp \left(\frac{V_x-V}{s_x}\right)},$$

$$h_{\infty }(V)=\frac{1}{1+\exp \left(\frac{V-V_h}{s_h}\right)}.$$

$$v=\frac{V}{k_V}\text{ and }{t}_s=\frac{t}{k_t},$$

$$\begin{array}{rl}\epsilon \frac{dv}{d{t}_s}& =-\left({\overline{g}}_{K}n\left(v-{\overline{V}}_{\mathrm{K}}\right)+{\overline{g}}_{Ca}{m}_{\infty }(v)h\left(v-{\overline{V}}_{\mathrm{Ca}}\right)\right)\equiv f(v,n,h),\\ \frac{dn}{d{t}_s}& =\frac{k_t}{{\tau }_n}\left(n_{\infty }(v)-n\right)\equiv g_1(v,n),\\ \frac{dh}{d{t}_s}& =\frac{k_t}{{\tau }_h}\left(h_{\infty }(v)-h\right)\equiv g_2(v,h),\end{array}$$

$$\begin{array}{rl}\frac{dv}{d{t}_f}& =f(v,n,h),\\ \frac{dn}{d{t}_f}& =\epsilon g_1(v,n),\\ \frac{dh}{d{t}_f}& =\epsilon g_2(v,h).\end{array}$$

$$\begin{array}{rl}\frac{dv}{d{t}_f}& =f(v,n,h),\\ \frac{dn}{d{t}_f}& =0,\\ \frac{dh}{d{t}_f}& =0,\end{array}$$

$$\begin{array}{rl}0& =f(v,n,h),\\ \frac{dn}{d{t}_s}& =g_1(v,n),\\ \frac{dh}{d{t}_s}& =g_2(v,h),\end{array}$$

$$S_0=\left\{(v,n,h):f(v,n,h)=0\right\}=\left\{(v,n,h):h=h_S(v,n)=-\frac{{\overline{g}}_{K}n(v-{\overline{V}}_K)}{{\overline{g}}_{Ca}{m}_{\infty }(v)(v-{\overline{V}}_{Ca})}\right\}.$$

$$L^{\pm }=\left\{(v,n,h)\in S_0:\frac{\partial f}{\partial v}=0\right\}.$$

$$\begin{array}{rl}\frac{dv}{d{t}_s}& =-{\left(\frac{\partial f}{\partial v}\right)}^{-1}\left(\frac{\partial f}{\partial n}{g}_1+\frac{\partial f}{\partial s}{g}_2\right),\\ \frac{dn}{d{t}_s}& =g_1,\end{array}$$

$$\begin{array}{rl}\frac{dv}{d{t}_d}& =\frac{\partial f}{\partial n}{g}_1+\frac{\partial f}{\partial s}{g}_2,\\ \frac{dn}{d{t}_d}& =-\left(\frac{\partial f}{\partial v}\right)g_1,\end{array}$$

$$M=\left\{(v,n,h)\in L^{\pm }:\frac{\partial f}{\partial n}{g}_1+\frac{\partial f}{\partial h}{g}_2=0\right\}.$$

$$s_{\max }=\lfloor \frac{\mu +1}{2\mu }\rfloor ,$$

Full text

Keywords  canard, action potential, early afterdepolarization, heart, myocyte

Introduction

In the normal heart, each heartbeat is associated with an action potential (AP). The cardiac AP consists of a depolarized phase in which the voltage is elevated; this is associated with transient increased permeability of the cell membrane to Na+ and Ca2+. The depolarized phase is followed by a repolarization to the resting membrane potential, associated with increased permeability to K+ ions. These changes in the membrane potential lead to a sequence of events that result in contraction of the heart muscle, thus allowing for the pumping of blood through the body.

Early afterdepolarizations (EADs) are pathological voltage oscillations that have been observed in heart muscle cells (cardiomyocytes) during the repolarizing phase of the cardiac AP under conditions in which the AP is elongated. EADs can be induced by hypokalemia [1, 2], as well as oxidative stress [3]. They have also been often observed following the administration of drugs that act on K+, Na+, or Ca2+ ion channels such as dofetilide [4], dl-sotalol [5], azimilide [5], bepridil [6, 7], isoproterenol [8, 9], quinidine [10], and BayK8644 [11, 2]. These drug-induced EADs can then lead to ventricular tachyarrhythmias [12, 13, 5]. Genetic defects in Na+ and K+ channels that prolong the action potential duration can also lead to an increased rate of EADs and risk of sudden death [14].

EADs have been associated with long QT syndrome [9], and have long been recognized as a mechanism for the generation of premature ventricular complexes (PVCs) in the electrocardiogram [15]. Different ventricular arrhythmias, including torsade de pointes, are thought to be initiated by PVCs stemming from EADs [16, 17, 18, 9]. That is, EADs at the myocyte level have been implicated as the primary mechanism promoting arrhythmias at the tissue level in acquired and congenital long-QT syndromes, including polymorphic ventricular tachycardia and ventricular fibrillation [19, 20, 5].

Numerous mathematical models have been constructed at the cellular level to study the genesis of EADs [21, 22, 2, 23, 24]. These have confirmed the importance of increased inward Ca2+ current and decreased outward K+ current in the production of EADs. They have also confirmed that reactivation of Ca2+ current is a key element of EAD production [24]. The importance of this “Ca2+ window current” in EAD production was later demonstrated through the use of the Dynamic Clamp technique [1], which is a hybrid between mathematical modeling and experimentation. Modeling at the tissue level has also been done, in this case to understand EAD propagation, synchronization and the genesis of arrthymia [25, 26, 27, 28]. These studies demonstrate that EADs at the cellular level can lead to arrhythmias at the tissue level, as has been suggested in experimental studies.

A useful analysis technique for understanding the behavior of models of excitable systems such as cardiomyocyte models separates system variables into those that change on a fast time scale and those that change on a slow time scale, and then analyzes the two subsystems and their interaction [29]. This slow/fast analysis has been used to understand the genesis of EADs, using a 3-variable model in which two variables were treated as “fast variables” and one variable treated as a “slow variable”. It was shown that EADs can arise via a delayed subcritical Hopf bifurcation of the fast subsystem of variables [23, 30]. This explanation, while providing insights, is limited in its descriptive capabilities. For example, it provides limited information on parameter sets for which EADs may occur, and it does not allow one to predict the number of EADs that are produced when they do occur.

Recently, it was demonstrated that EADs can be attributed to the existence of a folded-node singularity and the accompanying canard orbits [31]. This was done with the same three-dimensional model for cardiac action potentials, but now treating one variable as a fast variable and the other two as slow variables. Such a splitting provides the potential for insights that are not available with the 1-slow/2-fast splitting, as is demonstrated in [31] and in an earlier publication that focused on electrical bursting in pituitary cells [32]. For example, once it is established that the EADs are organized by a folded-node singularity, it is possible to determine regions of parameter space in which EADs can occur [31].

Ventricular cardiomyocytes are, in a physiological setting, subject to periodic stimulation from upstream cardiac cells, originating at the sinoatrial node. Prior experimental and modeling studies have demonstrated that EADs occur more readily at low pacing frequencies than at high frequencies [2, 24, 33]. At intermediate forcing frequencies the dynamics are very complex, consisting of alternans with varying numbers of EADs at each stimulus, a behavior described as “dynamical chaos” [2, 23]. The primary goal of this article is to provide an understanding for these phenomena. To achieve this, we use the same minimal cardiac action potential model that was developed in [2] and used recently in [31], and apply a 2-slow/1-fast splitting of the model. We demonstrate that the effects of periodic stimulation of the model cell can be understood precisely using the theory of folded-node singularities. In particular, we show that the number of EADs produced by a stimulus depends on where it injects the trajectory into the so-called “singular funnel”, and with this knowledge we demonstrate why low-frequency pacing is expected to yield more EADs than is high-frequency pacing. We also demonstrate the origin of the “dynamical chaos” that occurs at intermediate-frequency pacing. Finally, we demonstrate why drugs that inhibit the opening of K+ channels facilitate EADs, and why EADs can be induced by hypokalemia [1, 2, 5].

Action Potentials and EADs with the Minimal Model

We study a low-dimensional model for the electrical activity in a cardiomyocyte [2],

[TABLE]

where $I_{\rm K}$ is a repolarizing K+ current, $I_{\rm Ca}$ is a depolarizing Ca2+ current, and $I_{\rm sti}$ is an external pacmaking stimulus current. We note that (1) excludes the depolarizing Na+ current since prior studies have found that it has almost no effect on EADs (since it is inactivated during the plateau of the AP). Here, $V$ is the membrane potential across the cell, $n$ is the activation variable for the K+ channels, and $h$ is the inactivation variable for the L-type Ca2+ channels. The ionic currents are described by

[TABLE]

and we set

[TABLE]

where $H(\cdot)$ is the Heaviside function. That is, the stimulus current provides the system with square wave pulses of $1$ ms duration and $40$ $\mu$A/cm2 amplitude at a frequency set by the pacing cycle length (PCL). The steady state activation functions are

[TABLE]

where $x\in\{m,n\}$, and the steady state inactivation function is

[TABLE]

Standard parameter values are listed in Table Action Potentials and EADs with the Minimal Model; these have been tuned so that the model (1) periodically produces APs with EADs even in the absence of any stimulus current, as in [2].

The model cell (1) exhibits two distinct AP morphologies: regular APs and APs with EADs. For the remainder of the article, we use the Farey sequence notation, $1^{s}$, to denote a single large-amplitude AP with $s$ small-amplitude EADs during the repolarizing phase. Thus, a regular AP is denoted $1^{0}$ and an AP with 2 EADs is denoted $1^{2}$. More complicated rhythms are described using concatenations of these Farey sequences. For instance, a rhythm that periodically exhibits three regular APs followed by a single AP with 2 EADs is denoted $(1^{0})^{3}(1^{2})$.

Action Potential Duration and Number of EADs Increases with PCL

The model cell (1) is entrained to the periodic stimulus; for the parameter set in Table Action Potentials and EADs with the Minimal Model, the cell exhibits $1^{s}$ impulses with period set by the PCL. For small PCLs (i.e., high-frequency pulsing), the attractor is a $1^{2}$ rhythm (Fig. 1(a)). For intermediate PCLs ($1240$ ms $\lesssim{\rm PCL}\lesssim$ $1435~{}$ms), the cell exhibits complex EAD activity, including $1^{2}1^{3}$ alternans (Fig. 1(b)) and $1^{2}(1^{3})^{3}$ rhythms (Fig. 1(c)). For large PCLs (i.e., low-frequency pulsing), the cell is in a $1^{3}$ state (Fig. 1(d)).

We summarize the behaviour of the model cell and its response to periodic stimulation at various frequencies, by constructing a bifurcation diagram (Fig. 1(e)). To do this, we used a dynamic restitution protocol [34] in which the cell was paced at a fixed PCL until steady-state was reached, after which the action potential duration (APD) and PCL were recorded. We took the APD to be the amount of time the cell spends with $V>-70~{}$mV. With this choice of restitution protocol, the PCL is the sum of the APD and the diastolic interval, so our bifurcation diagram encodes the restitution curves (i.e., the plot of the APD as a function of the diastolic interval has the same qualitative features as shown in Fig. 1(e)).

The bifurcation diagram shows that the periodic stimulation elicits three types of behaviour. For high- and low-frequency stimulation, the model cell is in a purely $1^{2}$ or $1^{3}$ state, respectively. In the intermediate-frequency forcing range, the model cell has complex signature of the form $(1^{2})^{p}(1^{3})^{q}$, where $p$ and $q$ are integers. We observe that the AP signature becomes more complicated near the transition to the $1^{3}$ state. This increasing complexity of the AP signature near a transition is robust; it occurs for a wide range of $g_{K}$ and $g_{Ca}$ in (1) and has also been observed in other forced conductance-based cardiomyocyte models [2, 23].

Now that we have demonstrated the rich variety of dynamics present in the minimal model (1), we next investigate the dynamical mechanisms that underlie the observed rhythms. We use Geometric Singular Perturbation Theory [35, 36] as the basis of our analysis.

EADs Arise From Canard Dynamics

In this section, we show that the dynamical mechanisms responsible for the EADs are canards. To facilitate the analysis, we consider (1) with no stimulus input. A similar demonstration was provided by [31], but we elaborate on how the EADs emerge as the cell capacitance is increased from 0, (i.e., moving the system away from the singular limit), and we demonstrate how the underlying rotational sectors determine the number and duration of EADs. We first show that the model has a slow/fast structure. We use this slow/fast splitting to identify the geometric cast of characters involved in producing APs and EADs. We then demonstrate that folded node canards generate EADs, and these canards are robust in parameters. Finally, we demonstrate how drugs that inhibit K+ channels or a hypokalemic environment can facilitate EAD production.

The dynamics evolve over multiple timescales

A key observation is that the dynamics of the cell evolve over multiple timescales, with slow depolarized/hyperpolarized epochs interspersed with rapid transitions between them. We formally show this multi-timescale structure by introducing dimensionless variables, $v$ and $t_{s}$, via the rescalings

[TABLE]

where $k_{V}$ and $k_{t}$ are reference voltage and timescales, respectively. With these rescalings, the minimal model (1) becomes

[TABLE]

where $\overline{g}_{u}=\frac{g_{u}}{g_{\rm ref}}$ for $u\in\{\rm K,Ca\}$ denotes the dimensionless conductances with reference conductance $g_{\rm ref}$, $\overline{V}_{u}=\frac{V_{u}}{k_{V}}$ for $u\in\{\rm K,Ca\}$ denotes the dimensionless reversal potentials, and $0<\operatorname{\varepsilon}=\frac{C_{m}/g_{\rm ref}}{k_{t}}\ll 1$ is the ratio of the voltage timescale ($C_{m}/g_{\rm ref}$) to the reference timescale. The benefit of recasting the model in the dimensionless form (2) is that it reveals the typical timescales in the model. The voltage variable is fast with a timescale of $\frac{C_{m}}{g_{\rm ref}}\approx 5$ ms for $C_{m}=0.5~{}\mu$F/cm2 and $g_{\rm ref}=0.1$ mS/cm2. The activation variable, $n$, for the K+ channels is slow with timescale $\tau_{n}=80$ ms, and the inactivation variable, $h$, for the L-type Ca2+ channels is superslow with timescale $\tau_{s}=300~{}$ms. Thus, the system (2) is a three-timescale problem.

One effective approach to the analysis of multiple-timescale problems, as pioneered in the neuroscience context in [37], is Geometric Singular Perturbation Theory (GSPT). The idea of GSPT is to decompose a slow/fast system into lower dimensional slow and fast subsystems, analyze these simpler subsystems, and combine their information in order to understand the origin and properties of the dynamics of the original model. However, the GSPT approach is currently limited to two-timescale (i.e., slow/fast) problems. In three-timescale systems such as (2), a choice is usually made: to either group $v$ and $n$ together as ‘fast’, or to group $n$ and $h$ together as ‘slow’.

Prior studies of the minimal model chose to group $v$ and $n$ together as fast, whilst using $h$ as the sole slow variable [2]. In this 1-slow/2-fast approach, the EADs arise because the depolarized steady state of the $(v,n)$ subsystem loses stability via a Hopf bifurcation (with respect to $h$) leading to oscillations which are destroyed at a homoclinic bifurcation [27, 23, 38]. Whilst this mechanism is consistent with the in-vitro and in-silico observations that the EADs appear irregularly under periodic stimulation, it does not provide insight into how many EADs should be observed or why the number of EADs change with the PCL. Here, we take the alternative approach and treat $v$ as the only fast variable, whilst grouping $n$ and $h$ together as slow. We will show that this 2-slow/1-fast approach allows us to predict the maximal number of EADs that can be generated, and explain why the number of EADs changes with the PCL.

Underlying geometric structure

We now identify the geometric features that organize the EADs and APs. We begin by reformulating (2) in terms of the fast time, $t_{f}=\frac{1}{\operatorname{\varepsilon}}t_{s}$, which gives

[TABLE]

System (3) is equivalent to (2) in the sense that solutions of both systems trace out the same paths in the $(v,n,h)$ phase space, just at different speeds. We have seen that the dynamics of (1) alternate between slow and fast epochs. The rapid transitions between depolarized and repolarized phases are approximated by solutions of the 1D fast subsystem

[TABLE]

which is the approximation of (2) in which the slow variables move so slowly that they are fixed. (The fast subsystem is obtained by taking the singular limit $\operatorname{\varepsilon}\to 0$ in (3).) Similarly, the slow depolarized/repolarized segments of the dynamics are approximated by solutions of the 2D slow subsystem

[TABLE]

which is the approximation of (2) in which the fast voltage variable moves so rapidly that it (i) has already reached steady state and (ii) instantly adjusts to any changes in the slow gating dynamics. (The slow subsystem is obtained by taking the singular limit $\operatorname{\varepsilon}\to 0$ in (2).) Recall that the idea of GSPT is to analyze the 1D fast and 2D slow subsystems, and combine their information in order to understand the origin and properties of the dynamics in the full 3D system.

We begin with linear stability analysis of the 1D fast subsystem (4). The equilibria, $S_{0}$, of (4) form a cubic-shaped surface (in the $(v,n,h)$ space) called the critical manifold

[TABLE]

The outer sheets are stable and the middle sheet is unstable; these are separated by curves, $L^{\pm}$, of points corresponding to fold bifurcations of (4)

[TABLE]

For the cardiomyocyte model, the fold conditions (7) reduce to a set of lines on $S_{0}$ at constant voltage values (Fig. 2; red curves); $L^{+}$ denotes the fold curve at a depolarized voltage level, and $L^{-}$ denotes the fold curve at a hyperpolarized voltage that is the firing threshold. We note that the $V$-axis is also a fold curve (see ‘Folded Node and EAD Dynamics Are Robust’ section).

From the linear stability analysis, we conclude that most solutions of (4) end up on either the depolarized attracting sheet, $S_{0}^{a,+}$, or the hyperpolarized attracting sheet, $S_{0}^{a,-}$. Once trajectories reach one of these sheets, the slow dynamics dominate the evolution and the appropriate approximating system is the slow subsystem (5). The algebraic equation in (5) constrains the phase space to the critical manifold, whilst the differential equations describe the slow motions along $S_{0}$. Thus, the geometric singular perturbation analysis partitions the phase space into the fast dynamics away from the critical manifold together with the slow dynamics on the critical manifold. The critical manifold itself is the interface where the fast and slow subsystems interact.

For the slow evolution on $S_{0}$, we have differential equations to describe the motions of $n$ and $h$, whilst the algebraic equation implicitly describes the associated motions in $v$ (slaved to $S_{0}$; Fig. 2 black curves). To obtain an explicit description of the $v$-motions, we differentiate $f(v,n,h)=0$ with respect to the slow time, $t_{s}$, and use the graph representation of the critical manifold given in (6). This gives

[TABLE]

where $h$ has been replaced by $h_{S}(v,n)$. We stress that (8) is equivalent to (5); we have simply incorporated the restriction to $S_{0}$ explicitly by setting $h=h_{S}(v,n)$. In this formulation, it becomes clear that the slow flow is singular along the fold curves, $L^{\pm}$, where $\frac{\partial f}{\partial v}=0$. To deal with this finite-time blow-up of solutions, we perform the time rescaling $dt_{s}=-\frac{\partial f}{\partial v}dt_{d}$, which transforms the slow system (8) to the desingularized system,

[TABLE]

where again, $h=h_{S}(v,n)$. In this setting, the finite-time singularities of (8) along the fold curves have been transformed into nullclines of (9). Since the transformation that led to the desingularized system is phase space-dependent, some care must be taken when comparing trajectories of the desingularized system (9) with those of the true slow subsystem (8). On the attracting sheets, $S_{0}^{a,\pm}$, the flow of (9) is topologically equivalent to the flow of (8) since $\frac{\partial f}{\partial v}<0$ (and hence $t_{s}$ and $t_{d}$ have the same sign). On the repelling sheet, $S_{0}^{r}$, the flow of (9) is in the opposite direction to the flow of (8) since $\frac{\partial f}{\partial v}>0$ (and hence $t_{s}$ and $t_{d}$ have opposite signs).

With this relation between the slow and desingularized systems in mind, we now analyze the desingularized system in order to learn about the dynamics of the slow subsystem. The desingularized system possesses two types of equilibria or singularities. Ordinary singularities are isolated points such that $\{g_{1}=g_{2}=0\}$, and correspond to true equilibria of the desingularized system (9), of the slow subsystem (8), and of the original model (1). For the parameter set in Table Action Potentials and EADs with the Minimal Model, there is an ordinary singularity on $S_{0}^{r}$ (Fig. 2; cyan marker), corresponding to a saddle equilibrium.

Folded singularities, $M$, are isolated points on $L^{\pm}$ where the right-hand-side of the $v$-equation in (9) is zero, i.e.,

[TABLE]

Folded singularities correspond to equilibria of (9), however, they are not equilibria of the slow subsystem (8) or the original model (1). Instead, they are points where both the numerator and denominator of the right-hand-side of the $v$-equation in (8) vanish at the same time, so there may be a cancellation of a simple zero. This allows the possibility of solutions of the slow flow to cross the fold curves (via the folded singularity) with finite speed and move from an attracting sheet to the repelling sheet (or vice versa). Such solutions are called singular canards [39], and play important roles in applications. We refer to [40, 41] for extensive overviews of applications of folded singularities and canards in chemical, neural, and engineering contexts.

Folded singularities are classified as equilibria of the desingularized system. A folded singularity with real eigenvalues of the same sign is a folded node; with real eigenvalues of opposite signs it is a folded saddle; and with complex conjugate eigenvalues it is a folded focus. Folded nodes and folded saddles possess singular canards, whereas folded foci do not. The cardiomyocyte model possesses a folded node on $L^{+}$ for the standard parameter set (Fig. 2; green marker).

EADs originate from a folded node

We now demonstrate the origin of EADs in terms of the geometric structures identified in the prior section ‘Underlying geometric structure’. To motivate this, we first take a $1^{3}$ attractor of (1) (without periodic stimulation) and compare it to the critical manifold in the $(V,n,h)$ phase space (Fig. 3(a); magenta curve). The three EADs can be seen as small loops in the magenta trajectory about the upper fold curve, $L^{+}$. We observe that (i) the EADs are localized to the neighbourhood of the folded node; (ii) by decreasing $\operatorname{\varepsilon}$, or $C_{m}$, the EADs decrease in amplitude (compare curves of different colors in Fig. 3); (iii) by decreasing $\operatorname{\varepsilon}$, the location in phase space where the trajectory transitions from a depolarized state to a hyperpolarized state converges to the folded node. These observations lead us to conjecture that the EADs observed for $0<\operatorname{\varepsilon}\ll 1$ arise from the folded node itself.

How do the small oscillations emerge from the folded node? To answer this, we examine how the sheets, $S_{0}^{a,+}$ and $S_{0}^{r}$, of the critical manifold persist for small and nonzero $\operatorname{\varepsilon}$. As $\operatorname{\varepsilon}$ is increased away from zero, the attracting and repelling sheets, $S_{0}^{a,+}$ and $S_{0}^{r}$, perturb to attracting and repelling slow manifolds, $S_{\operatorname{\varepsilon}}^{a,+}$ and $S_{\operatorname{\varepsilon}}^{r}$, respectively [35, 36]. These slow manifolds are the surfaces to which the slow segments of trajectories of (1) are slaved. Both $S_{\operatorname{\varepsilon}}^{a,+}$ and $S_{\operatorname{\varepsilon}}^{r}$ are small and regular perturbations of $S_{0}^{a,+}$ and $S_{0}^{r}$, except in the neighbourhood of the folded node, where they instead twist around a common axis of rotation [39, 42]. The axis of rotation corresponds to the weak eigendirection of the folded node. The twisted slow manifolds are shown in Fig. 4 for various perturbations, corresponding to the $C_{m}$ values used in Fig. 3. (For the purposes of visualization, the slow manifolds have only been computed up to a plane, $\Sigma$, passing through the folded node. The method of computation is detailed in [43].) The twisting of the slow manifolds (and the slow flow on them) is confined to an $\mathcal{O}\left(\sqrt{\operatorname{\varepsilon}}\right)$ neighbourhood of the folded node [44]. Thus, the EADs arise from locally twisted slow manifolds around the folded node. This can be seen in the insets of Fig. 4, where the folded node is the intersection of the dashed blue curve (intersection of $S_{0}^{a,+}$ with $\Sigma$) and the dashed red curve (intersection of $S_{0}^{a,-}$ with $\Sigma$).

Canards organize the EADs

The local twisting of the slow manifolds results in a finite number of intersections between $S_{\operatorname{\varepsilon}}^{a,+}$ and $S_{\operatorname{\varepsilon}}^{r}$, called maximal canards. For the standard parameter set, we find that there are 5 maximal canards. The outermost, $\gamma_{0}$, is called the maximal strong canard and is the phase space boundary between those trajectories that exhibit EADs near the folded node and those that do not (Fig. 5). That is, any solution of the cardiomyocyte model (1) with initial condition to the left of $\gamma_{0}$ in Fig. 5 is a regular $1^{0}$ AP (Fig. 5(a) and (d); cyan curves). Any solution with initial condition between $\gamma_{0}$ and the secondary maximal canard $\gamma_{1}$ executes 1 EAD in the neighbourhood of the folded node (Fig. 5(b) and (d); beige curves). Any solution with initial condition enclosed by the secondary canards $\gamma_{1}$ and $\gamma_{2}$ exhibits 2 EADs around the folded node (Fig. 5(c) and (d); brown curves). In general, an orbit in the sector between the maximal secondary canards $\gamma_{k-1}$ and $\gamma_{k}$ will execute $k$ EADs. The innermost maximal canard, $\gamma_{w}$, is called the maximal weak canard and is the axis of rotation for both the slow manifolds and the other maximal canards. Thus, the maximal canards organize the EADs in phase space; the path taken by the trajectory relative to the maximal canards determines the number of EADs produced.

Folded Node and EAD Dynamics Are Robust

Given that the EADs arise from canard dynamics due to twisted slow manifolds around a folded node, is it possible to predict the number of maximal canards and associated EADs? The answer is ‘yes’, and it is encoded in the strong and weak eigenvalues, $\lambda_{s}<\lambda_{w}<0$, of the folded node (when considered as an equilibrium of the desingularized system). Let $\mu=\frac{\lambda_{w}}{\lambda_{s}}$ denote the eigenvalue ratio. Then, provided $\operatorname{\varepsilon}$ is sufficiently small and $\mu\gg\sqrt{\operatorname{\varepsilon}}$, the maximal number, $s_{\max}$, of EADs around the folded node is

[TABLE]

where $\lfloor\frac{\mu+1}{2\mu}\rfloor$ denotes the greatest integer less than or equal to $\frac{\mu+1}{2\mu}$ [44, 39]. The corresponding number of maximal canards is $s_{\max}+1$. For the folded node discussed in Figs. 2 – 5, we find $\mu\approx 0.13$ so that the maximal number of EADs that can be observed is $s_{\max}=4$, and there are 5 maximal canards (consistent with Fig. 4).

Not only does the formula (11) predict the number of EADs, it also predicts how the number of EADs changes with parameters. Bifurcations of maximal canards occur whenever $\mu^{-1}$ passes through an odd integer value [42]. That is, if the system parameters are varied so that $\mu^{-1}$ increases through $3$, then $s_{\max}$ increases from $1$ to $2$. If the system parameters are varied so that $\mu^{-1}$ increases through $5$, then $s_{\max}$ increases from $2$ to $3$, and so on.

There are two special cases, $\mu=0$ and $\mu=1$, where the folded node ceases to exist and hence the canard-induced EADs are eliminated. The resonance $\mu=1$ corresponds to the boundary where the folded node becomes a folded focus. Folded foci do not possess any canards. Hence, the $\mu=1$ resonance serves as the transition between regular $1^{0}$ APs and APs with EADs. This is illustrated in a two-parameter diagram, where the conductance of the K+ current ($g_{K}$) and the K+ Nernst potential ($V_{K}$) are varied and the asymptotic state of the system (1) is shown (Fig. 6). The red curve is the set of parameter values for which $\mu=1$. For parameter values within the region enclosed by the red curve the folded singularity is a folded focus, so only APs are produced (without EADs).

The dark green curves in Fig. 6 are parameter combinations such that $\mu=1/3$, so in the region delimited by these curves and the red $\mu=1$ curve there is a single maximal canard ($s_{\rm max}=1$) and APs with a single EAD are possible. On the olive curves $\mu=1/5$ and in the region delimited by these curves and the dark green curves APs with two EADs are possible. This process can be continued to higher odd integer values of $\mu^{-1}$; in the region between the olive curves and blue curves APs with three or more EADs are possible.

On the blue curves $\mu=0$. The $\mu=0$ resonance is known as a folded saddle-node (FSN) bifurcation and can occur in several distinct ways. The FSN bifurcation of type II (FSN II) is a bifurcation of the desingularized system in which a folded singularity and an ordinary singularity coalesce and swap stability in a hybrid transcritical bifurcation [45, 46]. That is, for $\mu>0$, the folded singularity on $L^{+}$ is a folded node and the ordinary singularity on $S_{0}^{r}$ is a saddle equilibrium. For $\mu<0$, the folded singularity on $L^{+}$ is a folded saddle and the ordinary singularity has moved to $S_{0}^{a,+}$ where it is a stable node. Hence, the FSN II bifurcation corresponds to the transition between EADs and stable depolarized steady states (Fig. 6; left blue curve).

The other way in which the FSN bifurcation can occur in the desingularized system is via a true transcritical bifurcation of folded singularities. That is, for $\mu>0$, there is a folded node on $L^{+}$ and there is a folded saddle on the $V$-axis. At $\mu=0$, the folded node and folded saddle coalesce, and for $\mu<0$, the folded singularity on $L^{+}$ is a folded saddle whereas the folded singularity on the $V$-axis is a folded node. The slow flow around the folded node on the $V$-axis is directed away from the $V$-axis, and so EADs will not be observed. Thus, for $\mu<0$, orbits of the slow flow encounter regular fold points on $L^{+}$, and the corresponding rhythm exhibits regular APs (without EADs). Hence, this FSN bifurcation corresponds to the transition between EADs and regular APs (Fig. 6; right blue curve). We note that, to the best of our knowledge, this type of FSN bifurcation (in which a pair of folded singularities undergo a true transcritical bifurcation) has not yet been reported or studied.

The two-parameter diagram (Fig. 6) illustrates that, in this model, there is a large set of $g_{K}$, $V_{K}$ parameters in which EADs can be produced. Thus, the behavior is generic, not limited to small regions of parameter space. It also illustrates the precision that GSPT provides in the determination of when EADs are possible, and the maximum number of EADs that are possible. Finally, the diagram shows that decreasing the K+ conductance, as is done with drugs like azimilide that act as K+ channel antagonists, can induce EADs (thick leftward arrow). Also, increasing the magnitude of the K+ Nernst potential, as in hypokalemia, can induce EADs (thick downward arrow). These observations are consistent with experimental studies [1, 2, 5].

Periodic Stimulation & Mixed-Mode Oscillations

We have established that EADs originate from canard dynamics around a folded node, and that the canards organize the EADs in both phase and parameter space. In this section, we restore the periodic stimulation and study the stimulus-driven EAD attractors. Our aim is to explain the bifurcation diagram in Fig. 1. We will show that the variety of AP morphologies exhibited under various PCLs can be explained by the canards.

High- and low-frequency pacing: canard-induced mixed-mode oscillations

Recall that when there is periodic stimulation, $I_{\rm sti}$, the system entrains to the driving oscillator. For low PCLs (i.e., high-frequency pacing), the attractor is a $1^{2}$ AP with EADs (see Figs. 1(a) and (e)). Using the results of our geometric analysis from the ‘EADs Arise From Canard Dynamics’ section above, we now deconstruct the $1^{2}$ rhythm (Fig. 7) and find that it consists of

1. (i)

canard-induced EADs around the folded node due to twisted slow manifolds, 2. (ii)

a fast transition from the depolarized folded node region to the hyperpolarized slow manifold, and 3. (iii)

a stimulus-driven transition from the hyperpolarized slow manifold to the depolarized slow manifold.

A representative example is shown in Fig. 7, where we compare the $1^{2}$ attractor to the slow and fast subsystems (panel (a)) and to the twisted slow manifolds (panel (b)). Note in Fig. 7(a) that the weak canard is approximately given by the stable manifold of the (cyan) saddle, and that the EADs are centered on this weak canard (i.e., the weak canard is the axis of rotation).

The periodic stimulus provides the mechanism for returning orbits to the neighbourhood of the folded node. More specifically, the stimulus switches ‘on’ during the slow hyperpolarized segment of the trajectory. This drives the orbit away from the hyperpolarized sheet before it can reach the lower firing threshold $L^{-}$. Moreover, the amplitude of the stimulus pulse is large enough that it pushes the orbit past the repelling sheet of the critical manifold and into the basin of attraction of the depolarized sheet, $S_{\operatorname{\varepsilon}}^{a,+}$. The timing of the stimulus is also such that the orbit is injected into the rotational sector enclosed by the maximal canards $\gamma_{1}$ and $\gamma_{2}$, and hence exhibits 2 EADs. This combination of a local canard mechanism (for the EADs) and a global (stimulus-induced) return mechanism is known as a canard-induced mixed-mode oscillation (MMO) [44].

Similarly, we find that for large PCLs (i.e., low-frequency pacing), the stimulus-driven $1^{3}$ attractor is a canard-induced MMO with period set by the PCL (see Fig. 1(d) and (e)). The $1^{3}$ MMO attractor consists of (local) canard-induced EADs around the folded node combined with a global stimulus-driven return that projects orbits from the hyperpolarized sheet into the rotational sector enclosed by the canards $\gamma_{2}$ and $\gamma_{3}$ (Fig. 8).

Intermediate-frequency pacing: EAD alternans due to reinjection into adjacent rotational sectors

In Fig. 1(e), we found that there is a band of intermediate pacing frequencies for which the stimulus-driven attractor is a $1^{2}1^{3}$ alternator (see Fig. 1(b) for a representative time series). We now compare the $1^{2}1^{3}$ attractor with the underlying geometric structures of the model cell (Fig. 9). As in the low- and high-frequency forcing cases, we find that the $1^{2}$ and $1^{3}$ segments are each canard-induced MMOs. The difference here is that the timing of the stimulus is such that the orbit visits different (contiguous) rotational sectors on each stimulus pulse.

The $1^{2}1^{3}$ alternans attractor, $\Gamma$, decomposes as follows. Starting on the hyperpolarized sheet, the first stimulus pulse (Fig. 9(a); leftmost cyan segment) projects the orbit $\Gamma$ into the rotational sector enclosed by $\gamma_{2}$ and $\gamma_{3}$ (Fig. 9(b); inset – black marker above $\gamma_{2}$). Thus, $\Gamma$ exhibits 3 EADs. After these 3 EADs are completed, the orbit transitions to the hyperpolarized sheet where it slowly drifts towards the firing threshold $L^{-}$. Before it can reach $L^{-}$, the next stimulus pulse (Fig. 9(a); rightmost cyan segment) projects the orbit into the rotational sector enclosed by $\gamma_{1}$ and $\gamma_{2}$ (Fig. 9(b); inset – black marker below $\gamma_{2}$), and thus $\Gamma$ exhibits only 2 EADs. The orbit then returns to the hyperpolarized sheet where it again slowly drifts towards $L^{-}$. Since $\Gamma$ only underwent 2 EADs, the APD is shorter (compared to the previous one) and the corresponding diastolic interval (DI) is longer. As such, the orbit is able to drift further along the hyperpolarized sheet before the next stimulus occurs. Once the stimulus ‘switches on’, the process repeats periodically, thus producing the $1^{2}1^{3}$ attractor.

Intermediate-frequency pacing: dynamical chaos and intermittency due to sensitivity near maximal canards

In Fig. 1(e), we found a band of intermediate pacing frequencies for which the model cell exhibited seemingly chaotic and intermittent behaviour. Here, we show that the complex EAD signatures arise from the crossing of a maximal canard. We do this for a representative $1^{2}(1^{3})^{3}$ attractor (Fig. 10), which we denote by $\Gamma$. As before, the individual APs with EADs are canard-induced MMOs. The variability in the number and magnitude of the EADs is due to the stimulus, which perturbs the orbit away from the hyperpolarized sheet at different locations on each pulse.

Let $\Gamma_{2}$ denote the $1^{2}$ segment of $\Gamma$, and $\Gamma_{3,j},j=1,2,3$ denote the $1^{3}$ segments of $\Gamma$, i.e., $\Gamma=\Gamma_{2}\cup\Gamma_{31}\cup\Gamma_{32}\cup\Gamma_{33}$. Starting with $\Gamma_{2}$, the stimulus (Fig. 10(c); cyan) induces a fast transition to the depolarized sheet close to $\gamma_{2}$ and in the sector between $\gamma_{1}$ and $\gamma_{2}$ (Fig. 10(b)), and hence there are 2 EADs. The intrinsic dynamics of the model cell return $\Gamma_{2}$ to $S_{\operatorname{\varepsilon}}^{a,-}$ where it slowly drifts to smaller $n$. The next stimulus initiates $\Gamma_{31}$ and causes the orbit to enter the depolarized phase in the rotational sector bound by $\gamma_{2}$ and $\gamma_{3}$. The additional EAD produced in $\Gamma_{31}$ extends the APD compared to that of $\Gamma_{2}$ (Fig. 10(e)). As such, the DI of $\Gamma_{31}$ is shorter than that of $\Gamma_{2}$. This means $\Gamma_{32}$ is initiated on $S_{\operatorname{\varepsilon}}^{a,-}$ at a larger $n$ value (Fig. 10(c)), and enters $S_{\operatorname{\varepsilon}}^{a,+}$ closer to $\gamma_{2}$. Since $\Gamma_{32}$ follows $\gamma_{2}$ more closely than $\Gamma_{31}$, (i) the resulting EADs are larger amplitude (Fig. 10(d)), (ii) the APD is longer, and (iii) the DI is shorter. Consequently, $\Gamma_{33}$ is initiated on $S_{\operatorname{\varepsilon}}^{a,-}$ at a larger $n$ value, enters $S_{\operatorname{\varepsilon}}^{a,+}$ closer to $\gamma_{2}$, and hence exhibits the (i) largest EADs, (ii) longest APD, and (iii) shortest DI.

The other complex MMO signatures reported in Fig. 1(e) emerge by the same mechanism. That is, the $(1^{2})^{p}(1^{3})^{q}$ attractors for $p,q\in\mathbb{N}$ arise because the PCL is such that the orbit enters the depolarized sheet close to the maximal canard $\gamma_{2}$. Since the behaviour of trajectories near a maximal canard is exponentially sensitive [42], small changes in the PCL manifest as significant changes in the number, amplitude, and duration of EADs on each pulse. Likewise, at such PCL values, small changes in initial conditions have large effects on the $V$ time course, the hallmark of chaos.

Discussion

It has been demonstrated previously that early afterdepolarizations produced by a simple cardiomyocyte model [2], a reduction of the Luo-Rudy 1 model [47], are the consequence of carnard dynamics in the vicinity of a folded node singularity [31], a result further illuminated through the geometric analysis shown in Figs. 2–5. We showed that these dynamics are robust in the $(g_{K},V_{K})$ parameter plane (Fig. 6). These parameters were chosen since they can be modulated by drugs or environment; $g_{K}$ is reduced by K+ channel antagonists such as azimilide, while $V_{K}$ is increased in magnitude in hypokalemia. Figure 6 predicts that both manipulations can induce EADs, and indeed both manipulations have been shown to do this in experiments [1, 2, 5].

The second set of results from our study involves the paced system, which receives periodic depolarizing stimuli (Fig. 1). Each stimulus pushes an orbit into the basin of attraction of the depolarized attracting sheet, triggering an action potential that can be a mixed-mode oscillation if EADs are produced. For high- and low-frequency pacing, orbits land in the rotational sectors delimited by the maximal canards and stay far from any of the maximal canards, so that the voltage time course exhibits regular, periodic behavior. At high stimulus frequencies, the orbits land in the rotational sector with 2 EADs, so each AP is a mixed-mode oscillation with 2 EADs. At low stimulus frequencies, the orbits land in the rotational sector with 3 EADs, so each AP is a mixed-mode oscillation with 3 EADs. The number of EADs depends upon the rotational sector in which the orbit lands in response to the stimulus (Figs. 7–10). The EAD alternans observed for intermediate-frequency pacing emerge because the stimulus current alternately projects the orbit into different rotational sectors on each pulse. In some cases, the outcome can be quite complex, with a sequence of mixed-mode oscillations of different durations and numbers of EADs. This behavior is what was referred to as “dynamical chaos” in earlier publications [2, 23].

The advantage of the minimal model for the analysis presented here is its low dimensionality. More realistic cardiomyocyte models can have 40 or more dimensions, reflecting many types of ionic currents and in many cases equations for Ca2+ handling in the cytosol, the sarcoplasmic reticulum (SR), and the subspace between the SR and the cell membrane [22, 48, 49, 50, 51]. One major advantage of these larger models is that they have more biological detail that allows for simulation of, for example, the application of pharmacological agents that act as antagonists for specific types of ion channels, such as inward-rectifying K+ channels, while the minimal model incorporates only a single type of K+ current and a single type of Ca2+ current. With the correct parameterizations, these more complete models are capable of reproducing the various forms of EADs that have been characterized, each with different, but partially overlapping, biophysical mechanisms [52], while the minimal model was developed to produce EADs of a particular type. EADs are divided broadly into types according to the timing of the events: “phase-2 EADs” occur during the plateau of an elongated AP, and “phase-3 EADs” occur during the falling phase of the AP. There are also “depolarizing afterdepolarizations” that occur after the completion of the action potential. The analysis that we performed herein on a minimal model suggests that the dynamics underlying some phase-2 EADs are canard induced, and we speculate that this will be the case in more complete biophysical models. While the full geometric singular perturbation analysis done with the minimal model is not possible with the high-dimensional models, it is possible to perform a less complete analysis, such as determining the existence of folded node singularities. Indeed, such an analysis is important for establishing that canard dynamics are the basis of phase-2 EADs in more complete models, and is currently being undertaken by our group.

Why does it matter whether EADs are due to canard dynamics near a folded node singularity? Although it sounds very abstract, the ramifications of knowing this can be very important and useful. As we have demonstrated, if the EADs are associated with a folded node singularity, then one can simply analyze the eigenvalues of the reduced desingularized system at the folded node to determine how many EADs are possible. Also, through analysis of the eigenvalues, one can determine parameter changes that will enhance EAD production or eliminate the EADs. In particular, one can determine regions of parameter space where canard-induced EADs are not possible, without the need to perform any numerical integrations (as in Fig. 6 and [31]). So once EADs are linked to folded node singularities, one gains a great deal of predictive capability. In addition to this, knowing the dynamical mechanism for the EADs helps in the understanding of complex behavior, such as dynamical chaos, that would be hard or impossible to understand from the viewpoint of interacting ionic currents (i.e., a biophysical interpretation). Knowing which ion channels are key players in EADs is of course important, and can provide targets for pharmacological or genetic manipulation, but the complexity of the multiscale nonlinear dynamical system provides limitations to interpreting behavior without mathematical tools such as GSPT.

The theory of folded singularities has been applied to numerous biological systems. This includes intracellular Ca2+ dynamics [53], the electrical activity of neurons [54, 55, 56] and pituitary cells [57], and mixed-mode oscillations that are likely canard-induced have been observed in the oxidation of platinum [58], dusty plasmas [59], and chemical oscillations [60, 61]. The demonstration that some forms of EADs are canard-induced, at least in a minimal cardiomyocyte model, adds cardiac cells to the growing list of the biological and chemical systems whose dynamics are organized by folded singularities. Our system is novel, however, in that it is periodically forced under normal (i.e., physiological) conditions, where the forcing is initiated at the sinoatrial node. As we demonstrated here, this forcing can lead to complicated dynamics due to the injection of the orbit into different rotational sectors, so that the number of EADs produced following each stimulus can vary. The result can appear to be unpredictable, and chaotic, and sensitive to small changes in the forcing frequency and initial conditions. Whether this complex behavior is exhibited in a physiological setting, within an intact heart, is unclear. It is generally accepted that EADs can lead to arrythmias [16, 17, 18, 9], including ventricular tachycardia, but it has not been establshed that complex, chaotic behavior at the single myocyte level contributes to this.

Conclusions

In this report, we showed the benefits of a 2-slow/1-fast analysis of a model for cardiac early afterdepolarizations. Knowing that the small EADs are due to canards organized around a folded node singularity not only explains the origin of the EADs, but provides a viewpoint through which one can comprehend important behaviors. For example, an analysis of the eigenvalues of the folded-node singularity provides information on the number of EADs that are possible for different parameter sets. It also explains why inhibition of K+ channels or a hypokalemic environment facilitates EAD production. Finally, it provides a solid basis for understanding the effects of periodic stimulation of cardiomyocytes. We used this technique to show why more EADs are generated at low-frequency pacing than at a higher pacing frequency. The technique was also used to explain the origin of complex alternan behavior that occurs with intermediate-frequency pacing. Overall, the use of slow-fast analysis provides information on the dynamics of this multi-timescale system that are hard or impossible to comprehend from a purely biophysical analysis (i.e., in terms of the effects of different ionic currents) or from computer simulations alone.