 Research Article
 Open Access
 Published:
Period doubling cascades of limit cycles in cardiac action potential models as precursors to chaotic early Afterdepolarizations
BMC Systems Biologyvolume 11, Article number: 42 (2017)
Abstract
Background
Early afterdepolarizations (EADs) are pathological voltage oscillations during the repolarization phase of cardiac action potentials (APs). EADs are caused by drugs, oxidative stress or ion channel disease, and they are considered as potential precursors to cardiac arrhythmias in recent attempts to redefine the cardiac drug safety paradigm. The irregular behaviour of EADs observed in experiments has been previously attributed to chaotic EAD dynamics under periodic pacing, made possible by a homoclinic bifurcation in the fast subsystem of the deterministic AP system of differential equations.
Results
In this article we demonstrate that a homoclinic bifurcation in the fast subsystem of the action potential model is neither a necessary nor a sufficient condition for the genesis of chaotic EADs. We rather argue that a cascade of period doubling (PD) bifurcations of limit cycles in the full AP system paves the way to chaotic EAD dynamics across a variety of models including a) periodically paced and spontaneously active cardiomyocytes, b) periodically paced and nonactive cardiomyocytes as well as c) unpaced and spontaneously active cardiomyocytes. Furthermore, our bifurcation analysis reveals that chaotic EAD dynamics may coexist in a stable manner with fully regular AP dynamics, where only the initial conditions decide which type of dynamics is displayed.
Conclusions
EADs are a potential source of cardiac arrhythmias and hence are of relevance both from the viewpoint of drug cardiotoxicity testing and the treatment of cardiomyopathies. The modelindependent association of chaotic EADs with period doubling cascades of limit cycles introduced in this article opens novel opportunities to study chaotic EADs by means of bifurcation control theory and inverse bifurcation analysis. Furthermore, our results may shed new light on the synchronization and propagation of chaotic EADs in homogeneous and heterogeneous multicellular and cardiac tissue preparations.
Background
Chaos can be defined as an aperiodic longterm behaviour in a deterministic dynamical system (either a differential equation or an iterated map/difference equation) that shows sensitive dependence on the initial conditions [1]. Though biological systems are affected by intrinsic and external stochastic noise, experimentally recorded irregular dynamics in the action potential (the characteristic membrane voltage response to a superthreshold electric stimulus) of cardiomyocytes have still been shown to be of the chaotic nature [2, 3]. More precisely, it has been observed in [2, 3] that by increasing the frequency of the stimulating current (or correspondingly by decreasing the pacing cycle length (PCL)), the 1:1 entrainment of the action potential is lost and a sequence of different m:n rhythms with alterations in the action potential duration (APD) called AP alternans is obtained before the dynamics finally become irregular. Using an iterated map
derived from the experimentally obtained restitution curve (APD plotted against the diastolic interval (DI)), it has been shown that, as the PCL of the stimulus is decreased, the slope of F at the fixed point APD^{∗} of Eq. (1) becomes progressively steeper until ∂ F(APD^{∗},PCL_{crit})/∂APD=−1 at some critical value PCL_{crit}. At this flip bifurcation point PCL_{crit} the stability of the fixed point is lost and a period2 cycle of Eq. (1) is born. The latter marks the beginning of a cascade of perioddoubling bifurcations, also compare with [4], that is produced as PCL is further decreased and that finally leads to chaotic sequences of APD generated by Eq. (1). This is in accordance with the PD route to chaos in the general theory of iterated maps [1]. For later reference we emphasize that [2] studied chaotic APD variations in periodically forced cardiac pacemakers cells (i.e., they show spontaneous action potential oscillations in the absence of a stimulus), while [3] studied chaotic APD variations in periodically forced nonspontaneously active Purkinje fibre and ventricular muscle cells.
Chaotic early afterdepolarizations (EADs) are a different type of irregular cardiac action potential dynamics (different from chaotic AP alternans) that have been experimentally observed both in periodically stimulated ventricular cardiomyocytes [5] and more recently in spontaneously beating human induced pluripotent stem cell derived cardiomyocytes (hiPSCCMs) [6–8]. EADs are abnormal voltage oscillations during the repolarization phase of the action potential, characterized by one or more periods of positive voltage slope before the normal repolarization is completed. While the irregular appearance of EADs has also been associated with stochastic activities of the ion channels (that regulate the action potential formation) in [9], the first mathematical evidence for the chaotic nature of EADs in dependence of PCL has been given in [5]. Using simulations of a deterministic differential equation model of ventricular action potential dynamics with a voltage equation
a restitution curve APD=r(DI,PCL _{ crit }), similar to the one shown in Fig. 2PP, was constructed. This curve was then used to derive Eq. (1) with F(APD,PCL)=r(PCL−APD,PCL) due to
Finally, it has been argued in [5] that the steep positive slope of the restitution curve r before the peak (which translates into a steep negative slope of the map F at the fixed point APD ^{∗}) proves the chaotic nature of the EAD dynamics.
In this paper, we demonstrate that a steep slope of the restitution curve r cannot serve as a general explanation of chaotic EAD dynamics displayed by cardiac AP models. Indeed, we will show that chaotic EAD dynamics in Eq. (2) are possible even if extracted ADP and DI data points do not form a function r (hence, not even admitting to speak of a slope). However, the key contribution of our paper will be novel insight into chaotic EAD dynamics gained from mathematical bifurcation studies of differential equation models of the form Eq. (2). Using a separation into fast and slow time scale variables, bifurcation analysis has been previously applied in [10] for the illumination of EADs in (a variant of) the periodically driven LR91model [11] for ventricular cardiomyocytes. In particular, it was shown in [10] that the fast subsystem of Eq. (2) features a supercritical Hopf bifurcation from which stable limit cycles emerge until they terminate at a homoclinic bifurcation of a saddle equilibrium. Then, EAD behaviour is obtained if the model parameters are set such that the state trajectory of the full system Eq. (2) temporarily coils around the limit cycle surface spanned between the supercritical Hopf and the homoclinic bifurcations of the fast subsystem. Furthermore, the homoclinic bifurcation in the fast subsystem has been introduced in [10] as the reason for the chaotic EAD dynamics that could be obtained whenever the PCL was chosen appropriately in an EAD featuring parameter setting. Further affirmations of the statement that chaotic EAD dynamics in periodically triggered action potentials are due to a homoclinic bifurcation in the fast subsystem of Eq. (2) are given in [12–15].
Recently, we have shown in [16] that EADs may occur in action potential models Eq. (2) that do not feature a supercritical Hopf bifurcation in their fast subsystem. In this paper, we now demonstrate that a homoclinic bifurcation in the fast subsystem of Eq. (2) is neither a necessary nor a sufficient condition for obtaining chaotic EAD dynamics.
We rather argue that a PD cascade of limit cycles in the full action potential system Eq. (2) paves the way to chaotic EAD dynamics in a modelindependent manner and present examples with models of a) periodically paced and spontaneously active cardiomyocytes, b) periodically paced and nonactive cardiomyocytes as well as c) unpaced and spontaneously active cardiomyocytes. Furthermore, we reveal that chaotic EAD dynamics may coexist in a stable manner with regular action potential dynamics, where only the initial conditions decide about the dynamics displayed. The results of our article on chaotic EAD dynamics in single cardiomyocytes may shed new light on the synchronization of chaotic EADs [5, 17] and EADmediated fibrillation [18] in cardiac tissue, may open new paths for the control of cardiac chaos [19] and may be of relevance within the CiPAinitiative [20] for a new approach to preclinical drug cardiotoxicity testing with hiPSCCMs, which considers EADs as a potential mechanismbased metric for the assessment of proarrhythmic risk.
Methods
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 pacemaker 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 timedependence 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 pacemaker 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 nonactive 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 nonactive 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
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
with the calcium channel conductance G _{ Ca } and the dynamic inactivation variable f, as well as the outward potassium current
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 nonactive cells we chose
This model was previously used in [10] for the study of chaotic EADs and is a slightly modified version of the LR91model [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
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})
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 auto07p [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 nonautonomous models PP and PV, the latter were transformed into autonomous models of increased state dimension. To this end we introduced the twodimensional dynamical system
that admits an asymptotically stable periodic orbit [33], and appended it to Eq. (2) with the setting
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 crosschecked 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
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
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
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 S1S2protocol, 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 S1S2protocol, see the Additional file 2.
Results
Simulations of chaotic EADs via the numerical integration of Eq. (2) have so far only been reported [5, 10, 12] for the case of a periodic stimulation by an external current I _{ sti }, see Figs. 1PP and 1PV for examples. However, studying drug induced EADs in spontaneously beating human induced pluripotent stem cell derived cardiomyocytes, we found that chaotic EADs may also form in simulations of AP models without periodic stimulation, see Fig. 1UP. Note that a comparable high number of small oscillations during EADlike activity has been experimentally observed in [36]. As the occurence of chaotic EADs of the type shown in Fig. 1PP (i.e., the AP is only triggered by the external current after full repolarization) has been attributed in [5, 12] to the steep slope of the APD restitution curve, see Fig. 2PP, we first wondered if this argument may also apply to chaotic EADs of the types shown in Figs. 1PV (external stimulation also before full repolarization) and 1UP (no external stimulation at all). Having derived the corresponding APD restitution curves, see Figs. 2PV and 2UP, we realized that they strongly deviate from their previously published counterparts as exemplified in Fig. 2PP. In particular, due to the lack of continuity and differentiability properties the “APD restitution curves” of Figs. 2PV and 2UP do not allow to define maps Eq. (1) for the iteration of APD and hence do not contribute to the understanding of the chaotic EAD types shown in Figs. 1PV and 1UP.
In our attempt to find a common explanation for the chaotic EAD dynamics observed in Figs. 1PP, 1PV and 1UP, we next focused on the hypothesis featured in [10, 12–15] according to which chaotic EADs have their source in a saddlehomoclinic bifurcation in the fast AP subdynamics. While in [10] the homoclinic bifurcation occurs after a supercritical Hopf bifurcation, our analysis of the fast subystems of models PP and PV shows that in these examples the homoclinic bifurcation is rather accompanied by a subcritical Hopf bifurcation of which only unstable limit cycles emerge, see Figs. 3PP and 3PV. However, the bifurcation analysis of the fast subsystem Eq. (10) of the AP model UP reveals that in this case a homoclinic bifurcation is not involved at all, see Fig. 3UP. It has been previously shown in [16], that the reason for the occurence of EADs in model UP is a saddlefocus fixed point whose unstable manifold causes small scale oscillations of growing amplitudes. Since the AP model UP still features chaotic EADs, see Fig. 1UP, it is demonstrated that a homoclinic bifurcation in the fast AP subdynamics is in fact not a necessary condition for the occurence of chaotic EADs.
Neither the EADtheory based on the steepness of AP restitution curves nor the EADtheory based on homoclinic bifurcations in fast AP subsystems can attribute the chaotic EAD dynamics of Figs. 1PP, 1PV and 1UP to a common cause. In particular, none of these theories can shed light on the chaotic EAD dynamics of model UP, as it neither admits a steep AP restitution curve nor a homoclinic bifurcation in the fast subsystem. For gaining insight, we hence decided to perform a bifurcation analysis of the full AP system of the model UP with the potassium channel conductance G _{ K } as the continuation parameter. While in principle also other model parameters could be chosen for the continuation, the choice of G _{ K } is motivated by the strong focus of the established drug safety guidelines on potencies to block potassium currents. The corresponding bifurcation diagram of Fig. 4UP shows that, starting from an area of low G _{ K } values that do not admit spontaneous oscillatory activity but only attraction towards fixed points, stable limit cycles of comparatively small amplitudes emerge from a supercritical Hopf point as the value of G _{ K } is increased. At the period doubling point PD1, the singleperiod oscillation loses its stability and splits into stable doubleperiod oscillations. As G _{ K } is further increased, further period doublings PD2, PD3, PD4 and PD5 are numerically detected. Figure 5 illustrates corresponding periodic trajectories with single, double, fourfold and eightfold period before the motion becomes chaotic. Furthermore, the chaotic nature is also present after the transition from the small amplitude motion with a failed repolarization to the large amplitude motion of AP type, with one representative given by the chaotic EAD dynamics displayed in Fig. 1UP, before the AP type motion turns into periodic EAD dynamics (followed by, not shown, periodic AP dynamics and finally an attraction towards a steady state after another Hopf bifurcation). For projections of the trajectories onto the Vxplane of the state space, see the Additional file 3.
Having associated the chaotic EAD dynamics of the model UP with the perioddoubling route to chaos in ODE systems [30], we wondered if this phenomenon, numerically and experimentally observed in many other biological and physical systems including neuronal activity [37, 38], may also underlie the chaotic EAD dynamics in the nonautonomous models PP and PV. Transforming the latter into autonomous systems using Eqs. (7) and (8), we hence performed a numerical bifurcation analysis again with G _{ K } as the continuation parameter. The resulting bifurcation diagrams, displayed in Figs. 4PP and 4PV, reveal that also in the case of model PP and PV, the value of G _{ K } corresponding to the chaotic EAD dynamics shown in Fig. 1 in fact lies in close vicinity to a cascade of period doubling bifurcations. For illustrations of the corresponding periodic trajectories we refer to the Additional file 4.
Discussion
Chaotic EAD dynamics in AP models Eq. (2) have been previously observed in the case of periodic pacing [5, 10, 12] but have been attributed to either the steepness of APD restitution curves or the existence of homoclinic bifurcations in fast AP subdynamics. While none of these results is able to explain the chaotic EAD dynamics observed in the unforced AP model UP, our study suggests the existence of a cascade of period doubling bifurcations of limit cycles as a modelindependent explanation for chaotic EAD dynamics both in forced and unforced cardiac AP systems of the ODE type. We emphasize that period doubling bifurcations, though then observed in iterated APD maps Eq. (1) rather than in mechanistic cardiac AP models Eq. (2), have so far only been linked with chaotic APD alternans [2–4] (which constitutes a type of cardiac arrhythmia that is different from EADs).
The results of this study were obtained by means of bifurcation analysis applied to AP models Eq. (2) based on numerical continuation using the software packages matcont [31] and auto07p [32]. In the context of EADsresearch, numerical continuation was previously only applied to fast subsystems of Eq. (2), see [10, 16]. Though the analysis of fast subsystems can illuminate EAD generating mechanisms during a single AP [10, 16], its capability to study the occurence of EADs over a time span covering several APs may be limited as the conditions for a separation into fast and slow state variables [34] may not be met in the long term. In that regard, the proper consideration of the stimulating current in Eq. (9) or Eq. (10), which certainly constitutes a very fast component of the system, might be one hurdle to be taken. In contrast, the advantage of studying the dynamical behaviour by a bifurcation analysis of the full AP model Eq. (2) is that its long term behaviour defined by stable and chaotic attractors can be captured and tracked as model parameters are continuously varied. Besides of the detection of bifurcations such as Hopf, saddle node of cycles, homoclinic or period doubling bifurcations, another benefit of bifurcation analysis with numerical continuation is that unstable (i.e., nonattracting) dynamical structures of Eq. (2) can be revealed. While the interpretation of the unstable limit cycles illustrated in Fig. 4 goes beyond the scope of this study, the relevance of unstable structures in the context of AP modelling is highlighted in [39], in which the existence of an unstable chaotic invariant set suggests that the excitability of a membrane to fire an AP may be more complex than a smooth hypersurface that divides subthreshold and suprathreshold membrane potentials.
The transition between two chaotic states (chaotic low amplitude dynamics of failed repolarization and chaotic EAD dynamics) observed in the model UP, see Figs. 4 and 5, is reminiscent of the transition between chaotic spiking and chaotic bursting [37] in an unforced Hindmarsh–Rose model of neuronal activity. Still, a difference is the sharp change in amplitude of the neighboring stable limit cycles in Fig. 4UP as opposed to the comparatively same levels of amplitudes in neuronal spiking and bursting reported in [37]. Such sharp transitions in stable limit cycle amplitudes as observed in the unforced model UP are also evident in the bifurcation diagrams of the periodically forced cardiac AP models PP and PV, see Figs. 4PP and 4PV. However, the chaotic EAD traces displayed in Figs. 1PP and 1PV are obtained with parameter values of G _{ K } that lie far away from these transitions but rather close to different cascades of period doubling bifurcations. A common feature of the unforced and the periodically forced models of this study is that the corresponding PD cascades are all of the supercritical type and seem to obey Feigenbaum’s law [30]
where \(G_{K_{n}}\) is the parameter value of G _{ K } corresponding to the nperiod doubling PD _{ n }.
A further exploration of the bifurcation diagram of the model PV reveals that the PD cascade is accompanied by a stable branch of limit cycles, indicated by the solid line in Fig. 6 a at V≈−82.7 mV. This branch corresponds to limit cycles of periodically driven regular action potentials as displayed in Fig. 6 b and demonstrates the coexistence of regular AP dynamics with the EAD affected limit cycles of the PD cascade for a certain range of potassium channel conductances G _{ K }. Consequently, regular AP dynamics may also coexist with chaotic EAD dynamics if G _{ K } is chosen beyond the cascadelimit of periodic behaviour, see Figs 1PV and 6b which are obtained with the same value of G _{ K } and identical periodic forcing but with the two different initial conditions Eq. (6) and
This is similar to the coexistence of periodic spiking and chaotic bursting reported for a nonforced spontaneously active neuron model in [40], then attributed to a bifurcation of a saddlenode periodic orbit. Hence, even though the fast subsystem of model PV features a homoclinic bifurcation, the latter is not sufficient for the occurence of chaotic EAD dynamics as in addition to the model parameters also the initial conditions need to be properly set. This coexistence is of relevance both for numerical and experimental studies of cardiac action potentials which typically are based on the assumption that after a sufficiently long transient the (one and only) “steady state” periodic behaviour is observed independently of the pacing history.
Other than channel conductances such as G _{ K }, the dependency of cardiac AP dynamics on the pacing cycle length PCL of the stimulating current I _{ sti } is of high relevance for studies of cardioelectrophysiology. Typically, the action potential duration APD is derived from voltage traces that are simulated or experimentally recorded for a discrete sequence of different values of PCL, see [23]. Our transformation of the periodically forced and hence nonautonomous AP model Eq. (2) into an autonomous AP model of extended state based on Eqs. (7), (8) offers a complementary method for studying the impact of PCL on the dynamic repertoire of the cardiac AP model. In particular, the numerical bifurcation analysis based on a limit cycle continuation allows for a continuous PCL screening and ensures that critical PCL points or intervals are not missed as they possibly would be in case of only a discrete PCL sampling. As an example, Fig. 7 shows corresponding bifurcation diagrams for the model PV with PCL as the continuation parameter. As with the continuation of G _{ K }, both the period doubling route to chaos and the coexistence of two stable periodic orbits for one and the same parameter value can be detected. Furthermore, the spread of additional PD bifurcations over a wide range of PCL suggests that the PDroute to EAD chaos is a nonlocal phenomenon in periodically driven AP models, which is in accordance with the observation of APD chaos over a wide range of PCL in the parametric sweep studies of APD [5, 12].
We end this section with some limitations of our study. Our study focuses on the generation of chaotic EADs in ODE models Eq. (2) that represent the behaviour of single cardiomyocytes. Clearly, our work needs to be extended to PDE models of cardioelectrophysiology in order to take into account the coupling between cells and spatial AP wave propagation. Furthermore, we have not incorporated stochastic effects which also may have an impact on the bifurcation repertoire of dynamical systems. Finally, the discontinuities in the bifurcation diagrams due to a failure of the numerical continuation and the unraveled unstable limit cycle branches require further analysis to extract the full information about the dynamical repertoire of Eq. (2). Note, however, that our bifurcation approach led to the discovery of chaotic EADs in unforced cardiac AP models and furthermore offers an explanation via the PD route to chaos then also uniformly applicable to chaotic EADs in forced cardiac AP models.
Conclusions
EADs are pathological voltage oscillations during the repolarization phase of the AP of cardiomyocytes and are considered as potential triggers of cardiac arrhythmias in context of both ion channel diseases and drug cardiotoxicity testing. In this study, we have contributed to the theory of EAD dynamics by attributing their chaotic appearance to cascades of period doubling bifurcations of limit cycles in deterministic AP models. As demonstrated in this article, the detection of PD cascades via the numerical continuation of limit cycles is possible both for paced and unforced cardiac AP models and serves as a strong indicator of chaotic EAD dynamics that then take place in immediate vicinity in the parameter space. Hence, the automatically executable detection of PD cascades in the full AP dynamics defines a parametertooutput map F:p→PD(p) which might allow to formulate and solve associated inverse bifurcation problems [25] that address the avoidance or the control [19] of chaotic EAD dynamics in cardiomyocytes. Furthermore, PD cascades might serve as classifiers of proarrhythmicity that provide improved risk prediction in comparison with purely simulation based and unspecific markers such as APD90 or AP upstroke velocity. Also, the incorporation of drugion channel interaction models [41] into Eq. (2) such as, e.g., simple pore block
would allow to conduct the bifurcation analysis directly with respect to drug concentration D or IC50 of the ion channels of interest. In that regard, the findings of this study might be of relevance for the currently unfolding CIPA initiative to redefine the drug safety paradigm [20]. Though the latter considers mathematical modelling and simulation of cardiac APs as one of its three pillars (next to ion channel studies and experiments with hiPSCCMs), it seems to so far ignore the potential of bifurcation analysis for the illumination of arrhythmic and chaotic behaviour in dynamical systems.
Abbreviations
 AP:

Action potential
 APD:

Action potential duration
 CiPA:

Comprehensive Invitro proarrhythmic assay
 DI:

Diastolic interval
 EAD:

Eearly afterdepolarization
 hiPSCCM:

Human induced pluripotent stem cell derived cardiomyocyte
 ODE:

Ordinary differential equation
 PD:

Period doubling
 PCL:

Pacing cycle length
 PP:

Paced pacemaker
 PV:

Paced ventricular
 UP:

Unpaced pacemaker
References
 1
Strogatz SH. Nonlinear Dynamics and Chaos. Boulder: Westview Press; 2015.
 2
Guevara M, Glass L, Shrier A. Phase locking, perioddoubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells. Science. 1981; 214(4527):1350–1353. doi:10.1126/science.7313693. http://science.sciencemag.org/content/214/4527/1350.full.pdf.
 3
Chialvo DR, Gilmour Jr RF, Jalife J. Low dimensional chaos in cardiac tissue. Nature. 1990; 343:653–7.
 4
Guevara MR, Ward G, Shrier A, Glass L. Electrical alternans and perioddoubling bifurcations. IEEE Comput Cardiol. 1984; 562:167–70.
 5
Sato D, Xie LH, Sovari AA, Tran DX, Morita N, Xie F, Karagueuzian H, Garfinkel A, Weiss JN, Qu Z. Synchronization of chaotic early afterdepolarizations in the genesis of cardiac arrhythmias. Proc National Acad Sci. 2009; 106(9):2983–988. doi:10.1073/pnas.0809148106. http://www.pnas.org/content/106/9/2983.full.pdf+html.
 6
Navarrete EG, Liang P, Lan F, SanchezFreire V, Simmons C, Gong T, Sharma A, Burridge PW, Patlolla B, Lee AS, Wu H, Beygui RE, Wu SM, Robbins RC, Bers DM, Wu JC. Screening druginduced arrhythmia using human induced pluripotent stem cell–derived cardiomyocytes and lowimpedance microelectrode arrays. Circulation. 2013; 128(11 suppl 1):3–13. doi:10.1161/CIRCULATIONAHA.112.000570.
 7
Hoekstra M, Mummery C, Wilde A, Bezzina C, Verkerk A. Induced pluripotent stem cell derived cardiomyocytes as models for cardiac arrhythmias. Front Physiol. 2012; 3:346. doi:10.3389/fphys.2012.00346.
 8
Khan JM, Lyon AR, Harding SE. The case for induced pluripotent stem cellderived cardiomyocytes in pharmacological screening. British J Pharmacol. 2013; 169(2):304–17. doi:10.1111/j.14765381.2012.02118.x.
 9
Lerma C, KroghMadsen T, Guevara M, Glass L. Stochastic aspects of cardiac arrhythmias. J Stat Phys. 2007; 128(12):347–74.
 10
Tran DX, Sato D, Yochelis A, Weiss JN, Garfinkel A, Qu Z. Bifurcation and chaos in a model of cardiac early afterdepolarizations. Phys Rev Lett. 2009; 102(25):258103. doi:10.1103/PhysRevLett.102.258103.
 11
Luo CH, Rudy Y. A model of the ventricular cardiac action potential. Depolarization, repolarization, and their interaction. Circulation Res. 1991; 68(6):1501–26. doi:10.1161/01.RES.68.6.1501. http://circres.ahajournals.org/content/68/6/1501.full.pdf+html.
 12
Sato D, Xie LH, Nguyen TP, Weiss JN, Qu Z. Irregularly appearing early afterdepolarizations in cardiac myocytes: Random fluctuations or dynamical chaos?Biophysical J. 2010; 99(3):765–73. doi:10.1016/j.bpj.2010.05.019.
 13
Weiss JN, Garfinkel A, Karagueuzian HS, Chen PS, Qu Z. Early afterdepolarizations and cardiac arrhythmias. Heart Rhythm. 2010; 7(12):1891–899. doi:10.1016/j.hrthm.2010.09.017.
 14
Qu Z, Xie LH, Olcese R, Karagueuzian HS, Chen PS, Garfinkel A, Weiss JN. Early afterdepolarizations in cardiac myocytes: beyond reduced repolarization reserve. Cardiovascular Res. 2013; 99(1):6–15. doi:10.1093/cvr/cvt104. http://cardiovascres.oxfordjournals.org/content/99/1/6.full.pdf.
 15
KroghMadsen T, Christini DJ. Nonlinear dynamics in cardiology. Annu Rev Biomed Eng. 2012; 14(1):179–203. doi:10.1146/annurevbioeng071811150106. PMID: 22524390. http://dx.doi.org/10.1146/annurevbioeng071811150106
 16
Kügler P. Early afterdepolarizations with growing amplitudes via delayed subcritical Hopf bifurcations and unstable manifolds of saddle foci in cardiac action potential dynamics. PLoS ONE. 2016; 11(3):1–14. doi:10.1371/journal.pone.0151178.
 17
Qu Z. Chaos in the genesis and maintenance of cardiac arrhythmias. Prog Biophys Mol Biol. 2011; 105(3):247–57. doi:10.1016/j.pbiomolbio.2010.11.001. Muscle ExcitationContraction Coupling: Elements and Integration
 18
Zimik S, Vandersickel N, Nayak AR, Panfilov AV, Pandit R. A comparative study of early afterdepolarizationmediated fibrillation in two mathematical models for human ventricular cells. PLoS ONE. 2015; 10(6):1–20. doi:10.1371/journal.pone.0130632.
 19
Garfinkel A, Spano M, Ditto W, Weiss J. Controlling cardiac chaos. Science. 1992; 257(5074):1230–1235. doi:10.1126/science.1519060. http://science.sciencemag.org/content/257/5074/1230.full.pdf.
 20
Fermini B, Hancox JC, AbiGerges N, BridglandTaylor M, Chaudhary KW, Colatsky T, Correll K, Crumb W, Damiano B, Erdemli G, Gintant G, Imredy J, Koerner J, Kramer J, Levesque P, Li Z, Lindqvist A, ObejeroPaz CA, Rampe D, Sawada K, Strauss DG, Vandenberg JI. A new perspective in the field of cardiac safety testing through the comprehensive in vitro proarrhythmia assay paradigm. J Biomol Screen. 2016; 21(1):1–11. doi:10.1177/1087057115594589. http://jbx.sagepub.com/content/21/1/1.full.pdf+html.
 21
Noble D. A modification of the Hodgkin–Huxley equations applicable to Purkinje fibre action and pacemaker potentials. J Physiology. 1962; 160(2):317–52. doi:10.1113/jphysiol.1962.sp006849.
 22
Hund TJ, Rudy Y. Rate dependence and regulation of action potential and calcium transient in a canine cardiac ventricular cell model. Circulation. 2004; 110(20):3168–174. doi:10.1161/01.CIR.0000147231.69595.D3.
 23
O’Hara T, Virág L, Varró A, Rudy Y. Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation. PLoS Comput Biol. 2011; 7(5):1002061. doi:10.1371/journal.pcbi.1002061.
 24
Paci M, Hyttinen J, AaltoSetälä K, Severi S. Computational models of ventricular and atriallike human induced pluripotent stem cell derived cardiomyocytes. Ann Biomed Eng. 2013; 41(11):2334–348. doi:10.1007/s1043901308333.
 25
Engl HW, Flamm C, Kügler P, Lu J, Müller S, Schuster P. Inverse problems in systems biology. Inverse Problems. 2009; 25(12):123014.
 26
Kügler P. A sparse update method for solving underdetermined systems of nonlinear equations applied to the manipulation of biological signaling pathways. SIAM J Appl Math. 2012; 72(4):982–1001. doi:10.1137/110834780.http://dx.doi.org/10.1137/110834780
 27
Xie Y, Izu LT, Bers DM, Sato D. Arrhythmogenic transient dynamics in cardiac myocytes. Biophys J. 2014; 106(6):1391–397. doi:10.1016/j.bpj.2013.12.050.
 28
MATLAB. Version 8.3.0.532 (R2014a). Natick: The MathWorks Inc.; 2014.
 29
Kuznetsov YA, Vol. 112. Elements of Applied Bifurcation Theory. New York: SpringerVerlag; 2013.
 30
Seydel R. Practical Bifurcation and Stability Analysis. New York: Springer; 2010.
 31
Dhooge A, Govaerts W, Kuznetsov YA. matcont: A matlab package for numerical bifurcation analysis of odes. ACM TOMS. 2003; 29:141–64.
 32
Doedel EJ, Fairgrieve TF, Sandstede B, Champneys AR, Kuznetsov YA, Wang X. AUTO07P: Continuation and bifurcation software for ordinary differential equations. 2007.
 33
Ermentrout B. Simulating, Analyzing and Animating Dynamical Systems. Philadelphia: SIAM; 2002.
 34
Kuehn C, Vol. 191. Multiple Time Scale Dynamics. New York: Springer; 2015.
 35
Hegger R, Kantz H, Schreiber T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos: An Interdisciplinary Journal of Nonlinear Science. 1999; 9(2):413–35. doi:10.1063/1.166424. http://dx.doi.org/10.1063/1.166424
 36
Song Z, Ko C, Nivala M, Weiss J, Qu Z. Calciumvoltage coupling in the genesis of early and delayed afterdepolarizations in cardiac myocytes. Biophysical J. 2015; 108(8):1908–921. doi:10.1016/j.bpj.2015.03.011.
 37
Wang XJ. Genesis of bursting oscillations in the HindmarshRose model and homoclinicity to a chaotic saddle. Physica D: Nonlinear Phenomena. 1993; 62(1):263–74. doi:10.1016/01672789(93)90286A.
 38
Jia B, Gu H, Li L, Zhao X. Dynamics of perioddoubling bifurcation to chaos in the spontaneous neural firing patterns. Cognitive Neurodynamics. 2012; 6(1):89–106. doi:10.1007/s1157101191847.
 39
Guckenheimer J, Oliva RA. Chaos in the Hodgkin–Huxley model. SIAM J Appl Dyn Syst. 2002; 1(1):105–14. doi:10.1137/S1111111101394040. http://dx.doi.org/10.1137/S1111111101394040
 40
Shilnikov A, Calabrese RL, Cymbalyuk G. Mechanism of bistability: Tonic spiking and bursting in a neuron model. Phys Rev E. 2005; 71:056214. doi:10.1103/PhysRevE.71.056214.
 41
Brennan T, Fink M, Rodriguez B. Multiscale modelling of druginduced effects on cardiac electrophysiological activity. Eur J Pharm Sci. 2009; 36(1):62–77. doi:10.1016/j.ejps.2008.09.013.
Acknowledgements
Not applicable.
Funding
Not applicable.
Availability of data and materials
All information necessary for reproducing the computational results is given in the article and its additional files. The manuscript does not involve data or material of other type.
Authors contributions
PK designed the study and wrote the manuscript. PK, MB and AE conducted research and discussed the manuscript. All authors read and approved the manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Additional files
Additional file 1
Transformation of Periodically Driven AP Models into Autonomous Form. (PDF 315 kb)
Additional file 2
Restitution Curves Computed via S1S2Stimulation. (PDF 170 kb)
Additional file 3
Period Doubling Route to Chaos of Model UP. (PDF 391 kb)
Additional file 4
Illustration of Period Doubling Route to Chaotic EADs in Models PP and PV. (PDF 1421 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Nonlinear dynamics
 Chaos
 Bifurcation theory
 Cardiac action potential
 Early afterdepolarizations