Early afterdepolarizations in cardiac action potentials as mixed mode oscillations due to a folded node singularity.

Early afterdepolarizations (EADs) are pathological voltage oscillations during the repolarization phase of cardiac action potentials. They are considered as potential precursors to cardiac arrhythmias and have recently gained much attention in the context of preclinical drug safety testing under the Comprehensive in vitro Proarrhythmia Assay (CiPA) paradigm. From the viewpoint of multiple time scales theory, the onset of EADs has previously been studied by means of mathematical action potential models with one slow ion channel gating variable. In this article, we for the first time associate EADs with mixed mode oscillations in dynamical systems with two slow gating variables and present a folded node singularity of the slow flow as a novel mechanism for EADs genesis. We derive regions of the pharmacology parameter space in which EADs occur using both the folded node analysis and a full system bifurcation analysis, and we suggest the normal distance to the boundary of the EADs region as a mechanism-based risk metric to computationally estimate a drug's proarrhythmic liability.

Introduction

Early afterdepolarizations (EADs) are abnormal depolarizations during the repolarization phase of the cardiac action potential and may be caused by drugs, oxidative stress or ion channelopathies. EADs are an important cause of cardiac arrhythmias such as polymorphic ventricular tachycardia (PVT) or torsades de pointes (TdP) [1], [2]. Recently, EADs propensity, expressed in terms of the net charge carried by major ionic currents during an action potential, has been chosen as an in-silico biomarker [3] for TdP risk evaluation of drugs within the CiPA (Comprehensive in vitro Proarrhythmia Assay) initiative [4], [5], [6] to overhaul the current cardiac drug safety regulations.

Computational models of cardiac cells [7] form the basis for the mathematical analysis of EADs. The most common approach is to numerically simulate single cell action potential models after on purpose or random model parameter variations such that EADs occur [8], [9], [3], [10], [11], [12]. Numerical simulation studies have also been performed at the cardiac tissue and cardiac organ level in order to analyze the synchronization of EADs and the triggering of arrhythmic patterns due to EADs [13], [14], [15]. Furthermore, numerical continuation has been used in [16] to explore the bifurcation structure of non-paced human ventricular myocyte models and to relate it to parameter regions of EADs occurence in the periodically paced model counterparts. EADs have also been investigated by bifurcation analysis applied to parametrized families of fast subsystems that are obtained from separating a single slow variable, related to ion channel gating, from the full cell model: while in [17], [18], [19], the onset of EADs has been linked to the existence of a supercritical Hopf bifurcation in the fast subsystem, it was shown in [20] that EADs may also arise in the presence of a subcritical Hopf bifurcation as well as in the complete absence of Hopf bifurcations in the fast subsystem.

In this study, we consider action potentials with EADs as mixed mode oscillations (MMOs), see [21] for a review, of signature 1 that are comprised of 1 large amplitude oscillation and s small amplitude oscillations. As an example, see Fig 1A for a 12 pattern corresponding to an action potential (AP) distorted by EADs and Fig 1B for a 10 pattern representing an undisturbed AP. We again take advantage of the difference in time scales between the dynamic variables of a cardiac cell model but apply—as opposed to previous multiple time scale investigations of EADs—a fast/slow analysis technique that features two slow variables. We demonstrate that the corresponding slow flow, which is restricted to a surface called critical manifold, admits a folded node singularity characterized by two real eigenvalues of the Jacobian that have the same sign. Then, the theory of MMOs [22], [23] implies that the trajectories of the full system are organized, locally near the folded node, by twisted slow manifolds [24]. This twisting organizes the small amplitude oscillations which in combination with a global return mechanism determined by the critical manifold gives rise to the MMO patterns that underlie cardiac action potentials with EADs.

Fig 1

Failure of standard fast-slow analysis to clarify disappearance of EADs.

A: Action potential with EADs in the form of a 12 mixed mode oscillation, obtained with the default parameter setting from Table 1. B: Action potential with normal repolarization in the form of a 10 mixed mode oscillation, obtained after changing G from 0.04 to 0.06. C: Standard fast-slow analysis for the default parameter setting. The bifurcation diagram for (2) with G as continuation parameter contains a z-curve with three branches that are connected via two saddle node bifurcations SN2 and SN1: a bottom branch of stable equilibria (solid line) of (2), a middle branch of unstable equilibria (dashed line), and a top branch of both stable and unstable equilibria of (2). At the top branch, the stability changes at a subcritical Hopf bifurcation subHB from which unstable limit cycles emerge that terminate at a saddle-homoclinic bifurcation HC. The blue dashed lines denote the maximum and minimum voltage values of the unstable limit cycles. The red line shows the projection of the 12-trajectories from A onto the bifurcation diagrams, the yellow line is the nullcline of the variable x. D: Standard fast-slow analysis for G = 0.06. The arrangement of the z-curve, the subcritical Hopf bifurcation and the x-nullcline is similar to C and does not explain why the EADs disappear in the transition from G = 0.04 to G = 0.06.

Table 1

Default parameter values.

Cm	ECa	GCa	EK	GK	VTf	kf	τf	VTx	kx	τx	VTd	kd
1	100	0.025	-80	0.04	-20	8.6	80	-40	-5	300	-35	-6.24
μF/cm2	mV	mS/cm2	mV	mS/cm2	mV	mV	ms	mV	mV	ms	mV	-mV

MMOs due to a folded node have been previously discussed, e.g., in the context of electrical bursting in pituitary cells [25], [26], [27] or climate change models [28], [29]. However, to the best of our knowledge, this study is the first that connects this particular MMO mechanism with EADs in cardiac action potentials and hence makes a novel contribution to EADs analysis. Furthermore, the conditions for the occurrence of such MMOs due to a folded node, namely existence of a folded node singularity and return of the large amplitude oscillation to the basin of attraction of the small amplitude oscillatory state, define borders of the EADs region in the model parameter space, which—further away from the singular limit—may be adapted by a complementary bifurcation analysis of the full system. This allows us to introduce the normal distance of a pharmaceutical compound to these EADs boundaries as a mechanism-based proarrhythmic risk metric, which may enrich the current discussion of computational biomarkers for the classification of a drug’s proarrhythmic liability within the CiPA initiative.

Materials and methods

Cardiac action potential model

Computational models of cardiac action potentials describe the dynamics of the transmembrane voltage V in dependence of the ionic currents and, for more than 50 years [30], [31] are an important tool in studying cardiac electrophysiology. Modern cardiac AP models [32], [9], [3], [12] consist of dozens of state variables and hundreds of model parameters, a complexity that results from the tremendous insight gained under the holistic paradigm of systems physiology but often hampers model analysis and validation. An alternative approach is to use parsimonious cardiac AP models [33], [34], [35], [36] that are low dimensional and only permit to address one or two phenomena but are amenable to mathematical analysis that goes beyond pure numerical simulation. In the context of EADs analysis, the model has been introduded in [37] and subsequently also used in [38], [20]. In particular, this ODE model of state dimension n = 3 features one inward calcium current with the calcium channel conductance G and the dynamic inactivation variable f, as well as one outward potassium current with the potassium channel conductance G and the dynamic activation variable x. The corresponding steady state variables are given by

The default parameter values from [37] are given in Table 1 and correspond to an action potential with EADs in form of a 12 mixed mode oscillation, see Fig 1.

Numerical methods

The model Eq (1) (and all of its subsystems, see the following subsections) were written in Matlab R2017b [39]. Initial value problems were integrated using the ode15s solver with relative and absolute error tolerances of 1e-10 and analytically derived Jacobian matrix, boundary value problems were solved using the bvp4c solver. Numerical curve continuation and bifurcation analysis was performed using both MATCONT 6.10 [40] and AUTO-07P [41].

Preliminaries for standard fast-slow analysis

The state variables of cardiac action potential models typically change at different time scales. With respect to the model (1) and the default parameter setting from Table 1, the time constants of the variables f and x are given by τ = 80 and τ = 300, respectively. The time constant of the variable V can be estimated as . Hence, τ < τ < τ, such that x is the slowest and V is the fastest variable.

The standard approach of studying EADs by multiple time scale analysis [17], [18], [20], [19] is to separate a single variable—namely the slowest one—from the faster ones. In case of (1), this yields a (2, 1)-fast-slow system with the fast subsystem where τ denotes the fast time scale. Hence, the slow variable x acts as a constant model parameter in the equations for the fast variables V and f. Formally, the fast subsystem (2) is obtained by taking the singular limit τ → ∞ in (1).

Preliminaries for (1, 2)-fast-slow analysis

An alternative to the standard fast-slow approach is to treat (1) as a (1, 2)-fast-slow system, in which not only x but also f is considered as a slow variable. Taking the singular limit C → 0, the trajectories of (1) converge during fast episodes to solutions of the fast subsystem or the layer equations where τ = t/C is the fast time scale. The equilibrium set of the system (3) is called the critical manifold

During slow episodes, the trajectories of (1) rather converge to solutions of the slow flow or the reduced system

Hence, the slow flow is described by a differential algebraic system whose phase space is given by the critical manifold (4).

Fig 2 gives an example of the critical manifold (4) for the AP model (1), and shows that, for the default parameter setting of Table 1, it is a folded surface with respect to the fast variable V. It consists of three sheets , S, separated by two disjoint fold curves F+, F−, i.e.,

Fig 2

Critical manifold C0 and singular periodic orbit s.

The figure shows the critical manifold (4) associated with the cardiac AP model (1) for two different values of the model parameter G. C0 is a folded surface in with attracting , and repelling sheets S separated by two fold lines F+ and F−. The projections of F+ and F− onto and are denoted by P(F+) and P(F−). Also shown is the singular orbit s that consists of slow segments on , and fast jumps between them. Inside the funnel region, which is bounded by the singular strong canard γ and F+, all trajectories approach the folded node FN allong the eigendirection e2 associated with the singular weak canard γ. SF denotes a saddle-focus equilibrium of the full system (1). A: G = 0.04. The singular orbit is projected into the funnel, which is a key reason for the appearance of EADs in the solution of the full system (1). B: G = 0.06. the singular orbit s lands outside of the funnel. Consequently, EADs do not appear in the solution of the full system (1).

Here, and are the attracting upper and lower sheets of C0 that are formed by stable hyperpolarized and stable depolarized steady states of (3). While for all points on or , the repelling sheet S is characterized by which corresponds to the unstable steady states of (3). Furthermore, the fold curves F+ and F− consist of all points on C0 with that satisfy a nondegeneracy condition , i.e.,

Taking the total time derivative of the constraint h(V, f, x) = 0, i.e, the slow flow (5) can also be described by restricted to C0. The system (8) is singular at the fold curves F+, F−, which means that the velocity of a trajectory goes to infinity as it approaches F+ or F−. The removal of the singularity by rescaling time according to leads to the desingularized system restricted to C0. The flow of (9) is equivalent to that of (8) on and but is reversed on .

The equilibrium points of the desingularized flow (9) are classified into ordinary singularities defined by (then also equilibrium points of the full system (1)) and into folded singularities that lie on a fold curve and satisfy

Given a folded singularity p* = (V*, x*), the eigenvalues λ1, λ2 of the Jacobian of (9) evaluated at p* decide whether p* is a folded node ( with λ1 ⋅ λ2 > 0), a folded saddle ( with λ1 ⋅ λ2 < 0) or a folded focus ( with ).

A key observation of this study, see the Results section, is that the desingularized slow flow associated with the cardiac action potential model (1) in its default parameter setting features a folded singularity of the folded node type.

Results

Failure of standard fast-slow analysis

A simulation of the cardiac model (1) with the default parameter values from Table 1 yields an action potential with EADs in form of a 12 mixed mode oscillation, see Fig 1. When the potassium channel conductance G is increased from 0.04 to 0.06, the EADs disappear and a normal action potential corresponding to a 10 oscillatory pattern is obtained. In a first attempt to understand the disappearance of the EADs we performed a standard fast-slow analysis [17], [18], [20], [19] and studied the bifurcations of the fast subsystem (2) with the slow variable x as numerical continuation parameter. The bifurcation diagram, see Fig 1, features a z-curve with three branches that are connected via two saddle node bifurcations SN2 and SN1: a bottom branch of stable equilibria (solid line) of (2), a middle branch of unstable equilibria (dashed line), and a top branch of both stable and unstable equilibria of (2). At the top branch, the stability changes at a subcritical Hopf bifurcation subHB from which unstable limit cycles emerge that terminate at a saddle-homoclinic bifurcation HC. The blue dashed lines denote the maximum and minimum voltage values of the unstable limit cycles. Next, we projected the 12- and the 10-trajectories onto the corresponding bifurcation diagrams and complemented the latter by the nullcline of the variable x (solid yellow line), which is defined by g2(V, f, x) = 0 or x = x∞(V), respectively. Above the x-nullcline, the trajectories move to the right due to dx/dt > 0, below, they move to the left due to dx/dt < 0 until they pass SN2 and are rejected, but in total, they only roughly follow the z-curve. Both for G = 0.04 and G = 0.06, the stable equilibria at the top branch are of the focus type, i.e., the eigenvalues are conjugate-complex with negative real part, but only for G = 0.04 the trajectory feels the attraction of the foci and is forced into a transient spiraliform movement that gives rise to the EADs. While both for G = 0.04 and G = 0.06 the x-nullcline crosses the z-curve in vicinity of SN1 to generate an unstable equilibrium of the full system (1), only for G = 0.04, the x-nullcline also intersects the branch of unstable limit cycles, which produces an unstable limit cycle of the the full system (1). Still, the EADs start at x-values much lower than that of this intersection such that based on Fig 1 the latter is not involved in the EADs genesis. Overall, the arrangement of the z-curve, the subcritical Hopf bifurcation and the x-nullcline is similar for G = 0.04 and G = 0.06, such that the standard fast-slow analysis fails to explain the destruction of the EADs in the transition from G = 0.04 to G = 0.06.

(1, 2)-fast-slow analysis reveals EADs genesis due to a folded node singularity

As an alternative to the standard fast-slow analysis, we next related the intermediate variable f to the slow variable x and determined the key objects of the resulting (1, 2)-fast-slow analysis.

Critical manifold, layer problem and desingularized system

The equation h = 0 that defines the critical manifold (4) is linear in f and can be solved for f in dependence on V and x, hence

For the default parameter setting from Table 1, (11) corresponds to a folded surface in that consists of two attracting sheets , and one repelling sheet S, see Fig 2. Note that the shape of C0 does not qualitatively change if G is set to 0.06, as G only acts as a scaling factor in (11). With respect to the fold curves F+ and F− from (7), we considered

Consequently, the constraint is satisfied if x = 0 or if

As for x = 0 the nondegeneracy condition is violated, we concentrated on (12), which is independent of G, and obtained the two solutions for the default parameter setting. Hence, the two fold curves are given by

Fig 2 shows how F+ and F− organize the partitioning (6) of C0 into attracting and repelling sheets.

The flow on the attracting sheets , is described by the desingularized system while transitions from one attracting sheet to another are described by the layer problem

In particular, trajectories of (15) that start on the fold curves F+ and F− land on the corresponding projections P(F+) and P(F−), see Fig 2.

Folded node singularity, funnel and singular periodic orbit

Next, we focused on the the folded singularities of the desingularized system (14) and constrained the defining Eq (10) to F±. This led to with V+ and V− given by (13) and f(V, x) defined in (11). Solving for x, we obtained and consequently two folded singularities

An evaluation of the Jacobian matrix of (14) at (V+, x+) and (V−, x−) showed that FS+ is a folded node singularity FN both for the default parameter values and the setting G = 0.06 due to with λ1 ⋅ λ2 > 0, see Table 2, while FS− is a folded focus that lies outside of the physiological range due to f(V−, x−) > 1.

Table 2

Folded singularities of the desingularized system.

GK	FS+	J(V+, x+)	λ1	λ2	FS−
0.04	(−24.79, 0.57, 0.68)	(-9.43·10^-4-2.02·10^-24.47·10^-6-2.7·10^-20)	−8.34 ⋅ 10−4	−1.08 ⋅ 10−4	(−73.51, 1.35, 0.05)
0.06	(−24.79, 0.43, 0.34)	(-9.75·10^-4-3.04·10^-27.44·10^-62.79·10^-20)	−5.96 ⋅ 10−4	−3.79 ⋅ 10−4	(−73.51, 1.34, 0.03)

Both for G = 0.04 and G = 0.06, the singularity FS+ is a folded node due to with λ1 ⋅ λ2 > 0. FS− lies outside of the pyhsiological range due to f(V−, x−) > 1, an evaluation of the Jacobian of (14) yields conjugate-complex eigenvalues λ1, λ2 and hence a folded focus singularity.

The so-called singular strong canard γ is a particular trajectory of the desingularized system (14) that reaches FN on along the eigendirection e1 related to the strong eigenvalue λ1, where strongness follows from |λ1| > |λ2|. The singular strong canard γ can be calculated by a shooting method and is plotted in Fig 2. Likewise, the singular weak canard γ is associated with the eigendirection e2 corresponding to the weak eigenvalue λ2. the singular strong canard γ and the fold line F+ bound a region of trajectories that is referred to as the singular funnel. Inside the funnel, all trajectories approach the folded node FN along the weak eigendirection e2, see the dashed line in Fig 2.

The final key object of the (1, 2)-fast-slow analysis is the singular periodic orbit s, that acts as a global return mechanism and is generated via a continuous concatenation of trajectories of the layer problem and the desingularized system. Starting at the folded node FN, the layer problem (15) is solved until the trajectory hits P(F+). From there, the desingularized system (14) is solved to continue the trajectory along until the fold line F− is reached. The jump from F− to P(F−) is obtained by again solving (15), after which the orbit is closed by returning the trajectory back to FN via the slow flow on , see Fig 2. As will be discussed next, the decisive difference between Fig 2A and 2B is that the singular periodic orbit s lands within the funnel for G = 0.04 while it lands and stays outside of it for G = 0.06. That way the value of G will decide whether EADs occur or whether they do not.

EADs due to the folded node

The critical manifold C0 and the singular periodic orbit s are objects that are defined in the singular limit C → 0. The theory of MMOs with multiple time scales [21], [42] characterizes perturbations of these objects away from the singular limit, i.e., for C > 0, and allows us to draw conclusions about the dynamics of the full system (1) with C > 0. In particular, for C > 0 and away from the folded node FN, the critical manifold smoothly perturbs to a slow manifold with attracting manifolds , and a repelling manifold according to geometric singular perturbation theory [43]. However, the theory of MMOs [44], [23], [24] implies that near the folded node FN, the critical manifold rather perturbs to twisted sheets. Fig 3B illustrates these twisted sheets in vicinity of the folded node of (14) for G = 0.04, which can be computed by continuation of solutions to boundary value problems associated with (1), see [24], [42]. Furthermore, the singular periodic orbit s, if injected into the funnel, perturbs to a trajectory of the full system (1) that flows from to in a rotating manner, see Fig 3. That way, a folded node singularity, in combination with a suitable global return mechanism, may give rise to cardiac action potentials with EADs. If s passes by the funnel, as in case of G = 0.06, then s perturbs to a relaxation oscillation of type 10, see Fig 1B, and a normal AP without EADs is obtained.

Fig 3

Twisted slow manifold near the folded node.

A: Periodic orbit of (1) with G = 0.04 in phase space. The default parameter setting leads to an action potential with EADs in form of a 12-MMO. B: The folded node causes a twisting of the slow manifolds and of the trajectory of (1) such that EADs occur. The attracting manifold and the repelling manifold were computed up to a surface Σ containing the folded node and transverse to the fold line F+. The red line is the projection of the trajectory from A.

Parameter region of EADs occurence

If one denotes the ratio of the weak and strong eigenvalues by the condition for the existence of a folded node singularity can be formulated as 0 < μ < 1. In this notation, the folded node is either destroyed at μ = 0, where λ2 passes through 0 and the folded node turns into a folded saddle, or at μ = 1, where λ1 and λ2 coalesce and the folded node turns into a folded focus. The second condition for the occurence of EADs, i.e., the reinjection of s into the funnel, can be quantified in terms of the distance δ of s from the singular strong canard γ. More precisely, δ measures the distance of the landing point of s on the projection P(F−) from γ along P(F−), see Fig 4 for an illustration of δ, and is positive if the landing point is inside of the funnel. Given these definitions of μ and δ, we constructed a bifurcation diagram that separates the (G, G)-parameter space into areas of different qualitative behaviour, see Fig 4. The region for which the folded node theory predicts the appearance of EADs is bounded by the curves defined by μ = 0 and δ = 0. Above the curve μ = 0, the theory predicts convergence to a depolarized steady state, while below the curve δ = 0, a relaxation oscillation of type 10 is expected.

Fig 4

EADs boundaries according to folded node theory.

A: The distance of the singular periodic orbit s from the singular strong canard γ, measured along the projection P(F−) of the fold curve F−, is denoted by δ. For δ > 0 (as depicted), s enters the funnel and perturbs to an action potential with EADs away from the singular limit. B: Two-parameter bifurcation diagram predicting the region of parameter space for which EADs appear. The folded node singularity exists for 0 < μ < 1, according to folded node theory the occurence of EADs in addition requires δ > 0. For parameter combinations of G and G that lie above the curve μ = 0, the theory predicts convergence to a depolarized steady state. For parameter combinations of G and G that lie below the curve δ = 0, a relaxation oscillation of type 10 is predicted.

Discussion

EADs in cardiac action potential models with different time scales may be considered as mixed mode oscillations of the form 1. MMOs arise in a variety of scientific areas [25], [26], [27], [28], [29], and MMO theory [21], [42] discusses at least four different local mechanisms that give rise to such a behaviour: i) the tourbillon mechanism of a dynamic Hopf bifurcation, ii) a saddle focus equilibrium that goes through a singular Hopf bifurcation, iii) the passage through a folded node, and iv) three-time-scale problems with a singular Hopf bifurcation. Previous results of EADs analysis can either be linked with the tourbillon mechanism, see [17], [18], [20] for the case of a dynamic supercritical Hopf bifurcation and [20] for the case of a dynamic subcritical Hopf bifurcation, or with the saddle focus mechanism, see [20] and also the discussion below. Here, we for the first time have linked EADs with the folded node mechanism by showing that cardiac AP models may feature a folded node singularity in the desingularized slow subflow that is combined with a rejection of the singular periodic orbit s into the funnel area. Away from the singular limit, the critical manifold near the folded node perturbs to a twisted slow manifold, see Fig 3, and its twisting finally organizes the small amplitude oscillations observed as EADs in the trajectory of the full system. One conceptual difference between i) and iii) is that i) is based on the separation of a single slow variable from the full system while in iii) two slow variables are separated. In the example discussed in this paper, the (1, 2)-fast-slow-analysis explained the occurrence of EADs in the parameter transition from G = 0.06 to G = 0.04, while the alternative (2, 1)-fast-slow-analysis could not give a proper illumination.

Overestimation of EADs region by folded node theory

One shortcoming of the (1, 2)-fast-slow-analysis is that the predictive power of singular perturbation theory is only guaranteed for C sufficiently small. While the sufficient smallness condition could be made slightly more precise in terms of , see [21], with ε = C/(G ⋅ τ) obtained from a nondimensionalisation of (1), it turns out that the region of EADs occurence is overestimated in case of the default parameter value C = 1, see Fig 5 for an illustration based on simulations of the full AP model (1) with different parameter values. For example, both (G, G) = (0.05, 0.025) and (G, G) = (0.034, 0.025) lie inside of the shaded EADs region of Fig 4, which is based on the folded node theory, but a simulation of (1) with C = 1 leads to a 10-oscillation corresponding to a normal action potential in the first case, and convergence to a hyperpolarized steady state in the second case. However, the simulations are in accordance with the prediction of Fig 4, if one reduces the capacitance to, e.g., C = 0.4, and action potentials with EADs of the form 11 and 19 are obtained. Note that in the singular limit C → 0, the trajectories converge to the singular periodic orbit s, see Fig 2. Hence, as C is further reduced, the amplitudes of the small scale oscillations become smaller and smaller until they no longer can be numerically detected.

Fig 5

Predictive power of folded node theory depends on smallness of C.

The figure shows simulations of the full AP model (1) with different combinations of G, G and C. A: G = 0.05, G = 0.025. The point (0.05, 0.025) lies in the EADs area of Fig 4, still, for C = 1 one obtains an oscillation of the l0-pattern. However, for C = 0.04 one obtains EADs in accordance with the folded node theory. B: G = 0.034, G = 0.025. The point (0.034, 0.025) lies in the EADs area of Fig 4, still, for C = 1 the trajectory reaches a hyperpolarized steady state. However, for C = 0.04 one obtains EADs in accordance with the folded node theory. C: Zoom into EADs area of B, which shows a 19-MMO pattern for C = 0.4.

In order to determine the actual EADs boundaries in the (G, G)-parameter space, we first performed a bifurcation analysis of the full system (1) with G as continuation parameter. Fig 6A shows the bifurcation diagram for the default parameter setting from Table 1 with C = 1. The curve corresponds to the equilibrium solutions of (1), which are unstable (dashed line) between the subcritical Hopf subHB and the supercritical Hopf bifurcation supHB. The coloured curves below and above represent the minimum and maximum voltages of stable (solid lines) and unstable (dashed lines) limit cycles of (1), which coalesce at saddle node of limit cycle bifurcations. Of particular interest is the sequence SNLC, s = 0, 1, 2, …, as its members mark the transition from 1 to 1-MMOs, that is from APs with s EADs to s + 1 EADs, as G is reduced. For instance, the default setting G = 0.04 is located in the red area between SNLC1 and SNLC2 and hence features a 12-limit cycle as shown in Fig 1. The unstable limit cycles that emerge from SNLC terminate at period doubling bifurcations of limit cycles PD. While cascades of PDs often are associated with chaos, see [45] for an illustration of the PD-route to chaotic EADs dynamics, we did not observe any chaotic behaviour in the particular system at hand. Due to numerical difficulties, only the first two unstable limit cycle branches could be detected. Still, we conjecture that the branching pattern is repeated with ever smaller distances as G is decreased until the limit cycle behaviour finally terminates at subHB. Given that EADs appear in the G-interval between SNLC0 and subHB, we next studied how that interval changes if the second channel conductance G is varied. Performing a sequence of curve continuations and tracking the locations of SNLC0 and subHB in the (G, G)-space gave rise to the region of EADs occurence shown in Fig 6C, which also visualizes the previously addressed overestimation obtained by the folded node theory.

Fig 6

Actual EADs boundaries derived from full system analysis.

A: Bifurcation diagram of the full system (1) with default parameter setting Table 1 and G as continuation parameter. Stable limit cycle oscillations of the type 10 are born at the supercritical Hopf bifurcation supHB, turn into 1-MMOs at the saddle node of limit cycle bifurcation SNLC and terminate at the subcritical Hopf bifurcation subHB. In particular, EADs occur in the G-range between SNLC0 and subHB. B: Zoom into A. C: Two-parameter bifurcation diagram showing the region of parameter space for which EADs actually appear. For comparison, dashed lines depict the EADs boundaries predicted by the folded node theory, see Fig 4B.

Distance to bifurcation as proarrhythmic risk metric

One possible benefit of calculating the region of EADs behaviour from computational cardiac AP models is that distances from its boundaries could be used in quantifying the proarrhythmic risk of pharmaceutical compounds. Currently, the CiPA initiative (Comprehensive in vitro Proarrhythmic Assay) [4], [5] features the use of mathematical modelling of cardiac APs as part of future preclinical drug safety testing. Recently, a computational proarrhythmic risk metric called qNet was introduced in [3] and defined as the net charge carried by the key ionic currents integrated from the beginning to the end of the simulated AP beat. Using statistical classification algorithms, qNet was shown to correctly separate reference compounds into different risk categories, where drugs of high risk are linked to low values of qNet and drugs of low risk are linked to high values of qNet, see Fig 5 in [3]. Furthermore, qNet was found to correlate with the system’s robustness against EADs in the sense that EADs propensity is higher the lower the value of qNet. In our opinion, a manifest and mechanism-based measure of EADs propensity is the normal distance in parameter space from the SNLC0-bifurcation at which EADs behaviour starts, see Fig 7A. The closer the system is to the boundary, the higher is the proarrhythmic risk. Surprisingly, CiPA’s risk metric qNet, if computed with the AP model (1), is a monotonically decreasing function of the normal distance, see Fig 7B. This would associate lower proarrhythmic risk with lower values of qNet, which is in complete contrast to the association found in [3]. While this observation only constitutes a preliminary result that also might be influenced by the simplicity of (1), it still suggests a model dependency of the link between qNet and EADs propensity and calls for further investigations.

Fig 7

Normal distance to bifurcation as proarrhythmic risk metric.

The figure illustrates the normal distance from the parameter region of EADs occurrence as a mechanism-based candidate measure of proarrhythmic risk. A: For instance, the AP system (1) with (G, G) = (0.042, 0.00835) is in normal distance 0.01 ms/cm2 from EADs occurrence. B: qNet computed with the AP model (1) is a monotonically decreasing function of the normal distance from SNLC0. The result is independent of the location chosen on the SNLC0-curve and in controversy to the qNet property found in [3].

Concurrence of folded node with saddle focus equilibrium

Finally, we return to Fig 5C and address the occurence of EADs in form of small scale oscillations with increasing amplitude. This behaviour is not an issue of the smallness of C but is also present for the default setting with C = 1 and, e.g., G = 0.035, see Fig 8A for a corresponding simulation of (1). While the first few oscillations are still organized by the twisted slow manifold associated with the folded node singularity FN, the growing oscillations are due to the saddle-focus equilibrium SF of the full system (1). More precisely, the trajectory approaches SF along its one-dimensional stable manifold, which is associated with the eigenvector e of the real eigenvalue λ < 0, and then spirals away along the unstable manifold with its tangential space M spanned by the pair of complex conjugate eigenvectors with positive real part, see Fig 8C. This mechanism has previously been related to the appearance of EADs in [20], but then independent of the folded node mechanism. As illustrated in Fig 8, the two mechanims may also coexist and consecutively shape the small scale oscillations, a phenomenon that was discussed in [21] in the context of the Koper model. Note that FN and SF also coexist for the parameter settings represented in Figs 2 and 3. However, the SF mechanism is only activated if the trajectory is properly conveyed to the SF after passage through the slow manifold, as exemplified in Fig 8B.

Fig 8

Concurrence of folded node with saddle focus equilibrium.

A: An action potential with EADs in form of a 115-MMO, obtained for the default parameter setting with C = 1 but G = 0.035. B: Visualization of the perodic orbit of A in the phase space. First, the EADs are organized by the twisted slow manifold due to the folded node FN, then EADs are caused by the spiraling along the unstable manifold of the saddle-focus equilibrium SF with tangential space M. C: Zoom into B.

Conclusion

In this paper, we have studied EADs in cardiac action potentials from the mathematical view point of mixed mode oscillations with multiple time scales. While the standard fast-slow analysis based on a single slow variable failed to explain the appearance of EADs in a parsimonious action potential model, a (1, 2)-fast-slow analysis based on two slow variables revealed that EADs may be caused by a folded node singularity if combined with a suitable global return mechanism. Furthermore, we have shown that the folded node coexists with a saddle focus equilibrium of the full system. Hence, EADs may result from the concurrence of two dynamical mechanisms, where the small scale oscillations first are induced by the twisted slow manifold associated with the folded node and then caused by spiraling along the unstable manifold of a saddle focus.

Our results form a novel contribution to EADs theory and may find application, e.g., in the context of preclinical cardiac drug safety testing. Both the MMO theory of folded nodes and a bifurcation analysis of the full system allow to compute regions in the ion channel conductance parameter space where EADs occur. Pharmacology data of a test compound such as IC50 values then defines its location in that space, and the normal distance from the boundary of the EADs region could be used as a mechanism based proarrhythmic risk metric. Interestingly since in contrast to [3], the normal distance associates EADs propensity with high values of the risk metric qNet, which is currently proposed by CiPA and defined based on statistical classification rather than on a EADs generating mechanism. In order to further investigate the suitability of the normal distance for preclinical cardiac safety testing, the concept needs to be tested by means of the the optimized IKr-dyn ORd model [3], which was used to define and validate the qNet risk metric. However, this in return requires to develop novel theoretical and computational tools for the MMO analysis of medium-to-large scale AP models.

State-of-the-art models [32], [9], [3], [12] comprise dozens of state variables and hundreds of model parameters. Hence, one limitation of our study is that we only dealt with a parsimonious cardiac action potential model, which in particular does not account for the sodium current, a known contributor to EADs [19]. However, it is the model’s simplicity that made it amenable to multiple time scale analysis [42], and the dynamical EADs mechanisms characterized by its help are likely to also be present in the more complex AP models. Future work will include the performance of tailored experiments with human induced pluripotent stem cell derived cardiomyocytes to challenge the EADs hypothesis generated in this paper. Furthermore, we aim to investigate whether a three-time-scale analysis of cardiac AP models can attribute EADs also with the remaining MMO-mechanism iv) of the listing given at the beginning of the discussion section.