Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 10771.
Published online 2017 September 7. doi:  10.1038/s41598-017-11355-1
PMCID: PMC5589958

Hysteretic Dynamics of Multi-Stable Early Afterdepolarisations with Repolarisation Reserve Attenuation: A Potential Dynamical Mechanism for Cardiac Arrhythmias


Some cardiovascular and non-cardiovascular drugs frequently cause excessive prolongation of the cardiac action potential (AP) and lead to the development of early afterdepolarisations (EADs), which trigger lethal ventricular arrhythmias. Combining computer simulations in APs with numerical calculations based on dynamical system theory, we investigated stability changes of APs observed in a paced human ventricular myocyte model by decreasing and/or increasing the rapid (I Kr) and slow (I Ks) components of delayed rectifying K+ current. Upon reducing I Kr, the APs without EADs (no-EAD response) showed gradual prolongation of AP duration (APD), and were annihilated without AP configuration changes due to the occurrence of saddle-node bifurcations. This annihilation caused a transition to an AP with EADs as a new stable steady state. Furthermore, reducing repolarisation currents (repolarisation reserve attenuation) evoked multi-stable states consisting of APs with different APDs, and caused multiple hysteretic dynamics. Depending on initial ion circumstances within ventricular myocytes, these multi-stable AP states might increase the local/global heterogeneity of AP repolarisations in the ventricle. Thus, the EAD-induced arrhythmias with repolarisation reserve attenuation might be attributed to the APD variability caused by multi-stability in cardiac AP dynamics.


Early afterdepolarisation (EAD)1, which is believed to trigger lethal arrhythmias such as torsade de pointes (TdP) in patients with long QT syndromes (LQTS)2, 3 and heart failure4, is a repolarisation abnormality in the cardiac action potential (AP). In cardiomyocytes, outward ionic currents that play a critical role in AP repolarisation are comprised of multiple components, including the slow (I Ks) and rapid components (I Kr) of the delayed-rectifier K+ channel currents, the inward rectifying K+ channel current (I K1), and Na+-K+ pump (NaK) current (I NaK). Although these multiple components for repolarisation appear somewhat redundant when viewing cardiomyocytes as a system, this redundancy confers robustness to AP repolarization; should diseases or drugs diminish any one of the outward currents (e.g., I Ks or I Kr), the others form an available reserve for AP repolarisation. This concept is referred to as the repolarisation reserve5, 6. Thus, a fundamental understanding of the relationship between the repolarisation reserve and EAD is essential if we are to gain new insights into the prevention of lethal arrhythmias, including drug-induced arrhythmias7.

Numerous experimental812 and theoretical studies1315 have shown that attenuating the repolarisation reserve caused the prolongation of AP duration (APD) and subsequent EAD development. Although it is believed that excessive APD prolongation produces EAD via destabilisation of the membrane potential (V m)10, 16, the effects of EAD generation on configuration changes in AP remain unclear. Moreover, a recently published study17 demonstrated that pharmacological inhibition of I Kr in rabbit ventricular myocytes (VMs) caused spontaneous and intermittent switching behaviour in the APs with and without EADs. Thus, the dynamical mechanism of EAD formation resulting from disruptions in the repolarisation reserve might be more complex than currently believed18, 19.

In the present study, we investigated dynamical stability changes of APs observed in a paced human VM model with I Kr and I Ks manipulation by combining computer simulations of APs with numerical calculations based on dynamical system theory, i.e. bifurcation theory20, 21. We show that EAD formation is caused not by V m destabilization resulting from the excessive prolongation of APs without EADs (no-EAD response), but by abrupt transition to a different steady state via a bifurcation due to the reduction of I Kr and/or I Ks. Our findings would provide a theoretical background for the prevention and treatment of EAD-induced arrhythmias in patients with LQTS and heart failure, as well as drug-induced arrhythmias.


Hysteretic dynamics of AP responses in the VM model

Figure 1 shows steady-state AP trains observed in the paced VM model when the I Kr conductance (%G Kr), which was expressed as a percentage of the maximal conductance value of I Kr (see Methods), was reduced or increased. First, %G Kr was set to 100% as the control condition, and the steady-state AP response was acquired. To investigate the effects of repolarisation reserve attenuation on AP response, %G Kr was gradually reduced in 1% intervals. For instance, when 1% was reduced from the control condition of %G Kr (i.e., 99%G Kr), the steady-state AP response was acquired by starting with the steady-state variable values of the AP response obtained just before the reduction of %G Kr (i.e., 100%G Kr; Supplementary Table S1) as an initial condition. The reduction of %G Kr until 11% caused AP prolongation in the no-EAD response, elevation of transient peaks and diastolic intracellular Ca2+ concentration ([Ca2+]i), and decline in intracellular Na+ concentration ([Na+]i) (Fig. 1, top). Subsequently, a slight reduction of %G Kr to 10% caused the no-EAD response to transform into AP with EAD, which in turn caused a single small depolarised potential during AP phase 2 (EAD1), where the transient [Ca2+]i peaks markedly increased and [Na+]i further decreased. Such a transition between different AP responses was also elicited between 2%G Kr and 1%G Kr; the EAD1 response transformed to an EAD2 response, characterised by two small depolarisations during AP phase 2.

Figure 1
Simulated hysteretic switching dynamics of early afterdepolarisations (EADs) in the ventricular myocyte model during I Kr manipulation. Steady-state action potential (AP) trains (top row), intracellular Ca2+ concentration, [Ca2+]i, (middle row), and intracellular ...

In contrast, when %G Kr was increased from 0%G Kr at which EAD2 response was observed, a transition from EAD2 to EAD1 occurred between 5%G Kr and 6%G Kr (Fig. 1, bottom). With increasing %G Kr from 5% to 6%, the diastolic [Ca2+]i and transient [Ca2+]i peaks declined while [Na+]i elevated. Further increasing %G Kr to 16% caused the disappearance of the EAD1 response, resulting in a transition to the no-EAD response. Of note, two distinct AP responses appeared at the same %G Kr value (refer to 11%, 5%, and 2%G Kr in Fig. 1). Such responses are known as a bi-stable phenomenon in which both of two distinct APs appears with the same parameter set, depending on an initial condition. Furthermore, the transitions between distinct AP responses observed in the decrease and increase of %G Kr exhibited hysteretic characteristics. This suggests that the AP state was highly dependent on ion concentrations (especially in [Na+]i) within the VM just prior to the manipulation of %G Kr.

Effects of changes in repolarisation reserve on dynamical stabilities of AP responses

Next, we investigated the effects of attenuating and enhancing the repolarisation reserve on the dynamical stabilities of AP. Figure 2 shows one-parameter bifurcation diagrams for each AP with 100%G Ks as a function of %G Kr. When I Kr was at the control condition, i.e. 100%G Kr (asterisk in Fig. 2Aa), a normal no-EAD response existed uniquely as a steady-state AP response. As %G Kr decreased (see horizontal blue arrows in each panel of Fig. 2), the no-EAD response prolonged the APD measured at 90% repolarisation (APD90) (Fig. 2B), while diastolic [Na+]i (Fig. 2C) and [Ca2+]i (Fig. 2D) declined and increased, respectively. The occurrence of saddle-node (SN) bifurcation at 10.8%G Kr (SN1) led to the annihilation of the no-EAD response. Immediately after the disappearance of the no-EAD response, EAD1 (as in Fig. 1) manifested as another stable periodic oscillation (vertical blue arrows annotated with a dagger in each panel of Fig. 2). Then, the number of transiently depolarised membrane potentials during AP phase 2–3 (#TDMP) was increased by one (Fig. 2Ab). At the same time, the transition to EAD1 caused the abrupt prolongation of APD, a marked decline in diastolic [Na+]i, and an increase in diastolic [Ca2+]i (Fig. 2B–D). As %G Kr further reduced, EAD1 disappeared together with an unstable solution due to the occurrence of another SN bifurcation at 1.43%G Kr (SN2), causing the transition to EAD2 (see vertical blue arrows annotated with double daggers in each panel of Fig. 2). This transition also led to further AP prolongation, decreased diastolic [Na+]i and increased diastolic [Ca2+]i. Notably, bifurcation analysis found that EAD3, a third type of AP response characterised by three small depolarisations during AP phase 2, coexisted with the EAD2 response. Whether EAD2 or EAD3 was elicited depended on the initial conditions, such as [Na+]i and [Ca2+]i.

Figure 2
Effects of decreasing I Kr on early afterdepolarisation (EAD) formation. One-parameter bifurcation diagrams of the number of transiently depolarised membrane potentials (#TDMP) during (A) action potential (AP) phase 2–3, (B) AP duration (APD) ...

In contrast, when %G Kr was increased (horizontal magenta arrows in each panel of Fig. 2) from a condition of I Kr loss (0%G Kr), the EAD2 response was destabilised via a period-doubling (PD) bifurcation at 5.2%G Kr (PD2). This caused a transition to the EAD1 (see vertical magenta arrows annotated with section signs in each panel of Fig. 2). Further increasing %G Kr led to the EAD1 destabilisation due to the occurrence of a PD bifurcation at 15.08%G Kr (PD1), causing a transition to a no-EAD response (pilcrows in each panel of Fig. 2). Thus, transitions between distinct AP responses were accompanied by jumps in APD90 and in [Na+]i and [Ca2+]i (vertical blue and magenta arrows). Furthermore, changes in dynamical stabilities of AP responses when %G Kr was reduced (horizontal blue arrows) or increased (horizontal magenta arrows) created hysteretic loops.

Ionic mechanism of bi-stable dynamics and EAD formation

Each corresponding AP response could be observed, depending on initial conditions, in the %G Kr ranges at which each solid black line overlapped in Fig. 2. We referred to those observations as bi-stable and/or multi-stable phenomena. Next, we investigated how the same %G Kr elicited bi-stable AP dynamics without and with EAD. Figure 3 shows each AP profile, [Na+]i, [Ca2+]i and ion currents associated with their ion concentration variabilities. Initial values, except for [Na+]i, were set to the same values at which the VM model exhibited the steady-state no-EAD response at 12%G Kr (Supplementary Table S2). When [Na+]i was set to a value (8.03 mM) near the state transition threshold, evoked APs converged to no-EAD responses (cyan trace in Fig. 3A), and [Na+]i increased due to a gradual accumulation of Na+ (cyan trace in Fig. 3B). This elevated [Na+]i enhanced NaK activity, which resulted in relatively large I NaK (Fig. 3C) via an increase in Na+ extrusion (Fig. 3D) during AP phase 2. Therefore, the outward current during the AP phase 2 was sufficient for the completion of AP repolarisation. In contrast, when [Na+]i was set to 8.02 mM (black trace in Fig. 3A), the evoked APs were gradually prolonged with each stimulation, and eventually converged to an EAD1 response, which caused a marked decline in [Na+]i (black trace in Fig. 3B). A prolongation of the active duration of I NaK due to AP prolongation increases Na+ extrusion. In particular, emerging EAD (asterisk in Fig. 3A) markedly increased Na+ extrusion (asterisk in Fig. 3D), consequently leading to a reduction in [Na+]i (Fig. 3B). On one hand, the Na+/Ca2+ exchanger (NCX) current, I NCX (Fig. 3E), contributed to increases in [Na+]i via Ca2+ extrusion and Na+ loading, particularly during the diastolic phase (daggers in Fig. 3D and F). However, increased [Na+]i was restricted because the diastolic interval (DI) was shortened by AP prolongation due to emerging EAD (see double daggers in Fig. 3A,D and F). Eventually, [Na+]i in the EAD1 response converged to a lower value than that in the no-EAD response (black traces in Fig. 3B). This lowered [Na+]i diminished NaK activity, thereby causing I NaK to decrease during AP phase 2 (Fig. 3C), and accelerating AP prolongation.

Figure 3
Ionic mechanism of the bi-stable dynamics between action potentials (APs) without and with early afterdepolarisation (EAD). Simulated APs and the changes in (A) membrane potential (V m), (B) [Na+]i, (C) I NaK, (D) net Na+ flux, J Na,net, (E) I NCX, ( ...

The AP prolongation mediated by the [Na+]i decline slowed AP repolarisation (Fig. 3A and Supplementary Fig. S1Aa), which delayed the voltage-dependent reductions and increases of the activation (d L) and inactivation gating variable (f L), respectively, in the L-type Ca2+ channel (LTCC) current (I CaL) (Supplementary Fig. S1Ab and Ac). The I CaL was defined by the product of d L, f L and the Ca2+-dependent inactivation (f CaL) state (Supplementary Fig. S1Ad); hence the deactivation delay of the d L-gate caused an augmentation of the I CaL window current during late AP phase 2 with each stimulation (section signs in Fig. 3G and Supplementary Fig. S1Ae). Furthermore, the AP prolongation increased Ca2+ influx via the LTCC (Fig. 3F), and led to the elevation of [Ca2+]i (Fig. 3H). The [Ca2+]i elevation enhanced NCX activity, consequently leading to the large inward shift of I NCX during AP phase 2 (Fig. 3E). The augmented I CaL window current and large inward shift of I NCX enhanced an inward current component in the net ionic current (I net) during late AP phase 2 (Supplementary Fig. S1Af). As the inward current components in I net increased with each stimulation, the I net temporarily caused a balance between the inward and outward currents (sharps in Supplementary Fig. S1A and S1B) that interrupted AP repolarisation. Then, the predominance of inward currents caused a transition from AP repolarisation to depolarisation during the AP phase 2. The initiation of transient depolarisation during the AP phase 2 in the EAD1 response coincided with the timing of the inward-outward balance of I net and the reactivation of the d L in I CaL (red dashed line in Supplementary Fig. S1A and S1B). Once EAD emerged, Ca2+ influx via LTCC markedly increased (Fig. 3G and H), causing elevation of Ca2+ concentration in the network ([Ca2+]NSR) and junctional sarcoplasmic reticulum ([Ca2+]JSR) (Supplementary Fig. S1Ag and Ah). This in turn led to the elevation of [Ca2+]i via increases in Ca2+ release from junctional sarcoplasmic reticulum to the myoplasm, i.e. Ca2+-induced Ca2+ release (CICR) (Supplementary Fig. S1Ai), resulting in the elevation of transient [Ca2+]i peaks (Fig. 3H).

Modification of EAD formations by attenuation of the repolarisation reserve

Next, we investigated the combined effects of decreases in both I Kr and I Ks on dynamical behaviours in the VM model. As I Ks decreased, all bifurcation points gradually shifted toward higher %G Kr values (Fig. 4 and Supplementary Fig. S2). The I Ks reduction also resulted in a right shift of the %G Kr range at which black solid lines overlapped in Fig. 2, which caused increasingly complicated multi-stability and hysteresis. The reason for the rightward shift of %G Kr, which causes EADs, can be intuitively understood as follows: the potential I Ks decline leads to the AP prolongation, increasing the I CaL window current and I NCX inward shift during AP phase 2. These inward current increases tend to cause inward-outward current balance in the I net during AP phase 2. Thus, the repolarization reserve attenuation accompanied by the decrease in I Ks facilitates the EAD formation with the I Kr reduction. Figure 5 shows representative examples of tri-stable dynamics that depended on initial conditions (Supplementary Table S3). The transitions among AP responses were caused by the dissipative or accumulative perturbations in [Na+]i. Reducing I Kr along with the decrease in I Ks formed a new low-amplitude voltage oscillation (LAVO), as shown in Fig. 5B. The LAVO was a V m oscillation at near plateau potential, and caused repolarisation failure, i.e. arrest at the depolarised potential. Four distinct stable APs (no-EAD, EAD1, EAD2, and LAVO) coexisted at a range of %G Kr (orange region in Fig. 4), indicating tetra-stable dynamics (Supplementary Fig. S3A). The VM model also exhibited tetra-stable dynamics (Supplementary Fig. S3B) at 70%G Ks (orange region in Supplementary Fig. S2), which converged to either no-EAD, EAD1, EAD2 or EAD3 responses, depending on the initial condition (Supplementary Table S4).

Figure 4
Effects of decreased I Ks on the formation of early afterdepolarisations. One-parameter bifurcation diagrams of the number of transiently depolarised membrane potentials (#TDMP) during (A) action potential (AP) phase 2–3, (B) AP duration measured ...
Figure 5
Representative examples of tri-stable action potential (AP) dynamics observed in the ventricular myocyte model. Simulated AP trains (top) and changes in intracellular Na+ concentration ([Na+]i) (bottom) at (A) 62%G Kr, and (B) 64%G Kr with 50%G Ks. The ...

The reduction of %G Ks caused complicated transitions between attractors. In Fig. 4 and Supplementary Fig. S2, the vertical blue and magenta arrows show the relationship of AP response transitions that occur as a result of a bifurcation when the %G Kr was reduced (horizontal blue arrows) or increased (horizontal magenta arrows). The annihilation of the no-EAD response caused by the SN1 bifurcation resulted in the transition to a LAVO response (see solid vertical blue arrow in Fig. 4), while the SN1 bifurcation of the no-EAD response in Supplementary Fig. S2 could transition to either EAD1, EAD2, or EAD3 responses (dashed vertical blue arrows), depending on the state variables just after the SN1 bifurcation. Furthermore, when the EAD1, EAD2, and EAD3 responses were annihilated by the occurrence of SN2, SN3, and SN4 bifurcations, respectively, the emerging AP response was not uniquely determined due to multi-stability (dashed vertical blue arrows in Fig. 4 and Supplementary Fig. S2). For example, the annihilation of the EAD1 response in Fig. 4 could potentially cause a transition to the no-EAD or LAVO response. In contrast, upon increasing %G Kr, the EAD1, EAD2, and EAD3 responses destabilised due to the occurrence of PD1, PD2, and Neimark-Sacker (NS) bifurcations, respectively. Furthermore, the PD4 or SN5 bifurcation destabilised or annihilated the LAVO response, respectively; for example, destabilisation of the EAD1 response causes a transition to either no-EAD, EAD2, or LAVO (dashed vertical magenta arrows in Fig. 4), and EAD2 destabilisation potentially transitions to a no-EAD or a LAVO response.

Contribution of [Na+]i and [Ca2+]i to multi-stable dynamics

To evaluate the impact of [Na+]i variability on multi-stable AP behaviour in the VM model, we performed additional bifurcation analyses in the non-autonomous [Na+]i-fixed system. Figure 6A shows one-parameter bifurcation diagrams of AP responses observed in the model with [Na+]i fixed at several values. When [Na+]i is fixed at 10 mM, neither multi-stable AP responses nor EAD formations occurred (Fig. 6Aa). Reducing %G Kr at relatively low [Na+]i (fixed at 8 mM and 6 mM) caused AP with EADs and bi-stable behaviour, but multi-stable AP responses did not occur. This suggests that the emergence of higher order multi-stable states depended on [Na+]i variation.

Figure 6
Contribution of [Na+]i- and [Ca2+]i-variability to the multi-stable dynamics. (A) One-parameter bifurcation diagrams of the number of transiently depolarised membrane potentials (#TDMP) during action potential (AP) phase 2–3 (top) and diastolic ...

Next we conducted AP simulations to examine differences between the mono-stable ([Na+]i = 10 mM) and bi-stable AP behaviours ([Na+]i = 8 mM) observed in the [Na+]i-fixed system at the same I Kr value (Fig. 6B). Except for [Ca2+]JSR and [Ca2+]NSR, initial values were adopted from a set of steady-state values of the state variables, at which both 10 and 8 mM [Na+]i-fixed systems exhibited no-EAD response at 15%G Kr (Supplementary Table S5). In both systems with low Ca2+ conditions ([Ca2+]JSR = [Ca2+]NSR = 0.01 mM), evoked APs were consistently no-EAD responses, including the transient response (black dashed V m traces in Fig. 6Ba and Bb), until steady state was obtained (black solid V m traces in Fig. 6Ba and Bb). However, the APD in the 8 mM [Na+]i-fixed system was longer than that of the 10 mM [Na+]i-fixed system (651 vs. 528 ms, respectively), even at the same I Kr (black V m traces in Fig. 6Ba and Bb). This was because the low [Na+]i (8 mM) attenuated NaK activity (black I NaK traces in Fig. 6Ba and Bb), decreased the outward current component of I net and augmented forward-mode NCX activity during AP phase 2 (black I NCX traces in Fig. 6Ba and Bb); the inward shift of I NCX increased the inward current component of I net. The low Ca2+ condition simulates a Ca2+-depleted myocyte, where a large beat-to-beat Ca2+ influx (black J Ca,net traces in Fig. 6Ba and Bb) leads to [Ca2+]i accumulation. Likewise, the [Ca2+]i peak increased with each stimulation and converged to a single value (Fig. 6Ca and Cb).

In contrast, the [Na+]i-fixed systems exhibited different responses in high Ca2+ conditions ([Ca2+]JSR = [Ca2+]NSR = 6 mM) depending on the [Na+]i-fixed value. The AP response in the 10 mM [Na+]i -fixed system transiently exhibited an EAD1 response (dashed cyan V m trace in Fig. 6Ba) that converged to a no-EAD response (solid cyan V m trace in Fig. 6Ba). In contrast, the AP response in the 8 mM [Na+]i-fixed system converged to an EAD1 response (cyan V m trace in Fig. 6Bb). In high Ca2+ conditions, Ca2+ efflux became dominant at the earlier phase of pacing stimuli (cyan dashed J Ca,net traces in Fig. 6Ba and Bb) because the high Ca2+ concentration in the myocyte augmented NCX activity (cyan dashed I NCX traces in Fig. 6Ba and Bb). In addition, high I NaK activity in the 10 mM [Na+]i-fixed system (I NaK in Fig. 6Ba) contributed to the increase in the outward current component of I net, leading to AP shortening and DI prolongation. This DI prolongation further increased Ca2+ efflux, causing a rapid decline in [Ca2+]i peaks (Fig. 6Ca). The low [Na+]i (8 mM) diminished I NaK (I NaK in Fig. 6Bb) and the high Ca2+ condition augmented forward-mode I NCX activity (I NCX in Fig. 6Bb); consequently, AP was prolonged and the I CaL window current was increased in late AP phase 2. Increases in these inward current components of I net eventually led to transient depolarisation during AP phase 2. AP prolongation with emerging transient depolarisation increased Ca2+ influx via reactivated I CaL (black I CaL and J Ca,net traces in Fig. 6Bb); hence the [Ca2+]i peaks converged to a relatively high value compared to the low [Ca2+]i condition (Fig. 6Cb). These results suggest that the [Na+]i variability was involved in the emergence of higher order multi-stable AP responses, and that which steady-state response appeared in bi-stable AP dynamical behaviour was mainly determined by the NCX activity strongly influencing [Na+]i and [Ca2+]i states.


The major findings of the present study are as follows: (1) attenuating the repolarisation reserve (I Kr and/or I Ks) caused prolongation of no-EAD responses, and the AP was annihilated without EAD formation via SN bifurcation. The EADs observed after the annihilation of AP emerged as a transition from the no-EAD response to a new stable state. (2) The destabilisation of AP responses that was mediated by the occurrence of PD and/or NS bifurcations caused the disappearance of APs with EADs during I Kr recovery. (3) These transitions among attractors with hysteresis were accompanied by the discontinuous change in various states, including [Na+]i, [Ca2+]i and APD, consequently creating the bi- and multi-stable states. (4) The development of multi-stable AP dynamics was strongly related to I NaK and I NCX changes, which were mediated by [Na+]i variability, and [Ca2+]i dynamics under the multi-stable AP behaviour contributed to EAD formation and APD variability. Such APD variability resulting from multi-stable AP dynamics in cardiomyocytes might be responsible for arrhythmias arising from repolarisation reserve attenuation.

The repolarisation reserve attenuation resulted in EAD responses with different numbers of transiently depolarised membrane potentials during AP phase 2–3 (Figs 1, ,55 and Supplementary Fig. S3). The evoked EAD patterns resembled those of erythromycin induced EADs recorded from the mid-myocardial cells in the canine left ventricle10. The EAD formations in our VM model were not due to destabilisation of V m during AP phase 2–310, 16 following excessive APD prolongation in no-EAD response with repolarisation reserve attenuation. The EAD emergence was due to the transition to the potentially-existing AP with EADs as a result of the annihilation of no-EAD response; therefore, EADs appeared to occur suddenly during %G Kr reduction. This shows that the AP with EADs is not directly linked to V m destabilization in the no-EAD response via AP prolongation. Furthermore, we found that the presented VM model exhibited multi-stable AP dynamics with hysteretic loops (Fig. 2, ,44 and Supplementary Fig. S2). In our previous study22, we confirmed the bi-stable AP behaviour was also observed in the more sophisticated human VM model proposed by O’Hara et al.23 (referred to ORd model). In rabbit VMs and the in silico model, Xie et al.17 have also observed intermittent spontaneous switching between AP responses with and without EADs that are similar to bi-stable AP behaviour with hysteresis, and they demonstrated that the [Na+]i variation played important roles in the occurrence of bi-stable dynamics. Such a hysteretic response in myocytes should be validated experimentally, but not yet at the moment. If a fine control of the I Kr/I Ks inhibition using drug concentration changes is technically feasible, it might demonstrate the presence of hysteretic dynamics in cardiomyocyte.

We showed that whether an AP with or without EAD were evoked under the bi-stable condition depended on the initial [Na+]i state (Fig. 3). Once steady-state AP responses in higher-order multi-stable behaviour was achieved, it was necessary to apply a relative large [Na+]i perturbation to cause state transitions among the different steady-state AP responses (Fig. 5 and Supplementary Fig. S3). Furthermore, the [Na+]i clamp in the [Na+]i-fixed system eliminated multi-stable behaviour in AP dynamics (Fig. 6A); thus, a wide [Na+]i dynamic range was required for APs with different APDs to coexist as multi-stable states. These results suggest that [Na+]i variability, mediated via NaK and NCX dynamics, is the underlying mechanism responsible for multi-stable AP dynamics. The NaK and NCX activity are relevant for Na+ extrusion amount upon AP development, and for the amount of Na+ loading during DI, respectively. Furthermore, AP prolongation under periodic stimuli results in DI shortening. Thus, AP prolongation increases the amount of Na+ extrusion via I NaK, whereas DI shortening decreases the amount of Na+ loading via I NCX. This extrusion and loading of Na+ generates the variation in [Na+]i. These imply that multi-stable AP dynamics depends significantly on the I NaK and I NCX models.

The initiation mechanism for the transiently depolarised membrane potential during AP phase 2–3 in the AP with EAD in the present study is essentially the same I CaL-dependent process that has been suggested in many previous studies2428. The low [Na+]i led to AP prolongation that was mediated by I NaK decrease, as shown in previous experimental29 and theoretical studies30. This AP prolongation caused [Ca2+]i elevation due to increase in Ca2+ influx via LTCCs resulting from the I CaL window current increase (Fig. 3G). Consequently, the [Ca2+]i elevation augmented forward-mode NCX activity (Fig. 3E), which increases Ca2+ extrusion as well as Na+ loading. This led to further prolongation of APD in a positive feedback manner17, resulting in an increase in AP repolarization delay. As the I CaL window current and inward I NCX shift increased with the AP repolarization delay, the inward-outward balance in I net was created, forming a transient depolarisation during AP phase 2 (Supplementary Fig. S1B). Based on this result, the suppression of inward I NCX shift might prevent the occurrence of transient depolarisation during AP phase 2 (i.e., EAD formation). Bourgonje et al.31 showed that combined NCX and LTCC inhibition was effective against TdP in dogs with chronic atrioventricular blocks. Thus, our results provide evidence that suppressing NCX activity might reduce the risk of EAD-induced arrhythmias.

The multi-stable AP behaviours may have significant implications for EAD-related arrhythmias. In addition to the I Kr inhibition in patients with LQTS type 1 (I Ks decrease) and type 2 (I Kr decrease)32, 33, the unintended I Kr inhibition by some drugs7, 10 may augment the risk of EAD development (Figs 2 and and4).4). Moreover, the hysteresis dynamics implies that, once EADs are elicited by I Kr inhibition, greater I Kr recovery (i.e. more than when the EADs were initiated) is required for suppression of the arrhythmogenic response. In addition, the emergence of multi-stable dynamics (Fig. 5 and Supplementary Fig. S3) may augment the APD heterogeneity of each myocyte in the ventricle. These predictions are based on the dynamical mechanism of EAD development, and are limited to the single cell model. As demonstrated by previous studies19, 3436, the relationship between EAD development in single cell and arrhythmogenicity in multi-cellular (tissue) models becomes more complicated due to the existence of electrotonic interaction (i.e. gap-junction coupling). Although the EAD development in mid-myocardial cell may augment the transmural (global) APD heterogeneity as they are much more vulnerable to EAD formation than the endocardial and epicardial myocytes13, 37, 38, the local APD heterogeneity in well-coupled myocytes in the intact heart may decrease due to the averaging effect resulting from the electrotonic interaction34, 35. In contrast, diminishing electrotonic interaction, for example with ageing-related fibrosis, is known to facilitate the EAD formations36, 39 and thereby enhance the local and global heterogeneity in AP repolarisation. This may create suitable conditions for reentry10 and lead to the generation of lethal arrhythmias such as TdP7, 40. Therefore, the effects of multi-stable AP behaviours on EAD-induced arrhythmias need to be studied in more realistic two- and three-dimensional ventricular models.

Using the relatively simple VM model was advantageous as we were able to utilize bifurcation analysis as a means to elucidate the dynamical mechanism of EAD emergences, but the incompleteness of the human VM model was a major limitation. For instance, our model lacked some ionic currents, such as late I Na, and the I Ks amplitude was relatively large compared to the experimental data obtained for basal human VMs12, 23. Our previous study22, however, showed that many bifurcation phenomena (i.e., mathematical structures in dynamical AP response) observed in the ORd model23 without pacing were essentially same as those of the presented VM model. The similarity of bifurcation structures between those models implies that the major dynamical properties of AP responses in the physiologically relevant human VM model may be represented by the relatively simple AP model. Nevertheless, more complex human VM models23, 4143 with refined formulas for I Ks and other components may be useful for understanding the detailed physiological and pathophysiological processes that are responsible for multi-stable AP phenomena. However, many of these complex human VM models are high dimensional systems, and to our knowledge, bifurcation analyses in such large-scale non-autonomous systems have not yet been successful. Such bifurcation analysis is theoretically possible, but is also practically challenging due to the extreme difficulty in obtaining an analytical description of the variational equations required for the evaluation of dynamical stability changes of cardiac AP responses (Supplementary Methods). Thus, dynamical system approaches in complex VM models would be highly challenging, and more research will be necessary to elucidate the precise roles of multi-stable AP behaviours in EAD-induced arrhythmias.


Ventricular myocyte model

As in our previous study22, we used the mid-myocardial cell version of a human VM model proposed by Kurata et al.44 that could reproduce phase 2 EADs during the inhibition of either I Kr or I Ks (i.e. attenuating the repolarisation reserve). The membrane currents include, I K1, I Kr, I Ks, I NaK, I CaL, I NCX, 4-aminopyridine-sensitive transient outward current (I to), Na+ channel current (I Na), Ca2+ pump current (I pCa) and background Na+ (I Nab) and Ca2+ (I Cab) currents. Time-dependent changes in the membrane potential (V m) are given by the formula:

dVm/dt = Istim − (ICaLIKrIKsItoINaIK1INabICabINaKINCXIpCa), 

where I stim represents the stimulus current (in pA/pF). The dynamics in the internal concentrations of Ca2+ ([Ca2+]i), Na+ ([Na+]i) and K+ ([K+]i) are modelled as first-order differential equations including material balance expressions. Details on the expressions of the VM model have been described previously22, 44. The model was implemented in an XML-based Physiological Hierarchy Markup Language (PHML), which is available at as an open-access resource. The external concentrations of Ca2+ ([Ca2+]o), Na+ ([Na+]o) and K+ ([K+]o) were fixed at 2, 140 and 5.4 mM, respectively. Furthermore, the [K+]i was fixed at 140 mM as per previous studies22, 45.

Bifurcation and analyses

The AP responses in the VM model were elicited with 80 pA/pF, 1 ms current pulses delivered at 0.5 Hz. Such a system can be defined as a periodic non-autonomous system, and dynamical responses in this system become periodic oscillations. Based on dynamical system theory20, 21, the study of qualitative properties of periodic oscillations in non-autonomous systems can be reduced to a diffeomorphism known as a Poincare map. Bifurcation occurs when a qualitative property (i.e. dynamical stability) of periodic oscillations changes due to altered system parameter values, e.g. the maximum conductance (g Kr) of I kr. We performed numerical calculations to detect bifurcation points using the homemade C language program as previously described4648. Details are provided in the Supplementary Methods. In briefly, bifurcation analysis proceeded as follows: first, we identified a stable AP response using numerical simulation. Numerical integration continued until the state variable values at each application of the stimulus converged to a steady-state value. The steady-state value was defined as the value when the difference between each sampled value of state variables obtained from k-th and (k + 1)-th APs in a simulated AP train was below 1 × 10−8. Next, we computed accurate state variable values of the steady-state AP using a shooting method (e.g. Newton’s method) by using the sampled state variable values obtained previously as the initial condition. The accurate state variable values were used to simultaneously compute the characteristic multipliers determining the dynamical stability of the AP response (see Supplementary Methods). Then, we changed the parameter value (e.g. g Kr) slightly, and executed Newton’s method again using the preceding result as the new initial condition. If Newton’s method converges, we can obtain new state variable values and characteristic multipliers for the changed parameter value. The shooting method was performed for every changes in system parameter values until bifurcations were detected with the assessment of the characteristic multipliers; this method is known as continuation21. The bifurcations that may arise in the paced VM model are saddle-node (SN), period-doubling (PD), and the Neimark-Sacker (NS) bifurcations21. At the SN bifurcation point, two periodic APs (e.g. stable and unstable AP responses) coalesce and annihilate. The PD and NS bifurcations cause the stability change in a periodic AP response. As a side effect, another AP response with a doubled period is generated around the AP that caused the PD bifurcation. In NS bifurcation, a quasi-periodic AP response may occur in the VM model. Detailed descriptions for each bifurcation are provided in the Supplementary Methods.

Bifurcation diagrams

In the present study, we constructed one-parameter bifurcation diagrams of periodic AP responses for changes in the g Kr of I kr. The maximum conductances of I Ks and I Kr per unit area were set to 0.0257 nS/pF and 0.00738 nS/pF, respectively22, which we defined as the control values (g Ks and g Kr). Throughout the article, the maximum conductance of I Ks and I Kr was expressed as a percentage of the control values, i.e. %G Ks and %G Kr, respectively. Each of the state variables could be plotted on the vertical axis in a one-parameter bifurcation diagram. In the present study, diastolic [Na+]i and/or [Ca2+]i values that were sampled during each application of the stimulus were plotted as a function of %G Kr. Furthermore, we depicted one-parameter bifurcation diagrams of the %G Kr in relation to APD90 and the number of transiently depolarised membrane potentials during AP phase 2.

Data availability

All data generated or analysed during this study are included in this published article (and its Supplementary Information files).

Electronic supplementary material


This work was supported by JSPS KAKENHI Grant Numbers 24790214, 16KT0194 (K.T.) and 17H04018 (Yo.K.), the Ichiro Kanehara Foundation (K.T.), the Takeda Science Foundation (K.T.), the Hiroshi and Aya Irisawa Memorial Promotion Award for Young Physiologists from the Physiological Society of Japan (K.T.), MEXT KAKENHI Grant Numbers 25136720 (Ya.K) and 22136002 (Yo.K.), Research on Regulatory Harmonization and Evaluation of Pharmaceuticals, Medical Devices, Regenerative and Cellular Therapy Products, Gene Therapy Products, and Cosmetics from Japan Agency for Medical Research and Development (AMED), Grant Number 16mk0104007h9903 (K.F.).

Author Contributions

Author Contributions

K.T. and Ya.K. conceived and designed the experiments; K.T. conducted the simulations, performed numerical calculations and prepared figures; K.T., Ya.K., and K.F. analysed the results; K.T., Ya.K., K.F. and Yo.K. drafted and edited the manuscript. All authors reviewed the final version of manuscript.


Competing Interests

The authors declare that they have no competing interests.


Electronic supplementary material

Supplementary information accompanies this paper at doi:10.1038/s41598-017-11355-1

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Contributor Information

Kunichika Tsumoto,

Yoshihisa Kurachi,


1. Cranefield PF. Action potentials, afterpotentials, and arrhythmias. Circ Res. 1977;41(4):415–423. doi: 10.1161/01.RES.41.4.415. [PubMed] [Cross Ref]
2. Chiang CE, Roden DM. The long QT syndromes: genetic basis and clinical implications. J Am Coll Cardiol. 2000;36(1):1–12. doi: 10.1016/S0735-1097(00)00716-6. [PubMed] [Cross Ref]
3. Jons C, et al. Use of mutant-specific ion channel characteristics for risk stratification of long QT syndrome patients. Sci Transl Med. 2011;3(76):76ra28. doi: 10.1126/scitranslmed.3001551. [PubMed] [Cross Ref]
4. Undrovinas AI, Maltsev VA, Kyle JW, Silverman N, Sabbah HN. Gating of the late Na+ channel in normal and failing human myocardium. J Mol Cell Cardiol. 2002;34(11):1477–1489. doi: 10.1006/jmcc.2002.2100. [PubMed] [Cross Ref]
5. Roden DM. Taking the “idio” out of “idiosyncratic”: predicting torsades de pointes. Pacing Clin Electrophysiol. 1998;21(5):1029–1034. doi: 10.1111/j.1540-8159.1998.tb00148.x. [PubMed] [Cross Ref]
6. Nattel S, Maguy A, Le Bouter S, Yeh YH. Arrhythmogenic ion-channel remodeling in the heart: heart failure, myocardial infarction, and atrial fibrillation. Physiol Rev. 2007;87(2):425–456. doi: 10.1152/physrev.00014.2006. [PubMed] [Cross Ref]
7. Belardinelli L, Antzelevitch C, Vos MA. Assessing predictors of drug-induced torsade de pointes. Trends Pharmacol Sci. 2003;24(12):619–625. doi: 10.1016/ [PubMed] [Cross Ref]
8. Roden DM, Hoffman BF. Action potential prolongation and induction of abnormal automaticity by low quinidine concentrations in canine Purkinje fibers. Relationship to potassium and cycle length. Circ Res. 1985;56(6):857–67. doi: 10.1161/01.RES.56.6.857. [PubMed] [Cross Ref]
9. Takanaka C, Singh BN. Barium-induced nondriven action potentials as a model of triggered potentials from early afterdepolarizations: significance of slow channel activity and differing effects of quinidine and amiodarone. J Am Coll Cardiol. 1990;15(1):213–221. doi: 10.1016/0735-1097(90)90205-4. [PubMed] [Cross Ref]
10. Antzelevitch C, Sun ZQ, Zhang ZQ, Yan GX. Cellular and ionic mechanisms underlying erythromycin-induced long QT intervals and torsade de pointes. J Am Coll Cardiol. 1996;28(7):1836–1848. doi: 10.1016/S0735-1097(96)00377-4. [PubMed] [Cross Ref]
11. Burashnikov A, Antzelevitch C. Acceleration-induced action potential prolongation and early afterdepolarizations. J Cardiovasc Electrophysiol. 1998;9(9):934–948. doi: 10.1111/j.1540-8167.1998.tb00134.x. [PubMed] [Cross Ref]
12. Jost N, et al. Restricting excessive cardiac action potential and QT prolongation: a vital role for IKs in human ventricular muscle. Circulation. 2005;112(10):1392–1399. doi: 10.1161/CIRCULATIONAHA.105.550111. [PubMed] [Cross Ref]
13. Viswanathan PC, Rudy Y. Cellular arrhythmogenic effects of congenital and acquired long-QT syndrome in the heterogeneous myocardium. Circulation. 2000;101(10):1192–1198. doi: 10.1161/01.CIR.101.10.1192. [PubMed] [Cross Ref]
14. Silva J, Rudy Y. Subunit interaction determines IKs participation in cardiac repolarization and repolarization reserve. Circulation. 2005;112(10):1384–1391. doi: 10.1161/CIRCULATIONAHA.105.543306. [PMC free article] [PubMed] [Cross Ref]
15. Vandersickel N, et al. A study of early afterdepolarizations in a model for human ventricular tissue. PLoS One. 2014;9(1):e84595. doi: 10.1371/journal.pone.0084595. [PMC free article] [PubMed] [Cross Ref]
16. Liu GX, et al. Differential conditions for early after-depolarizations and triggered activity in cardiomyocytes derived from transgenic LQT1 and LQT2 rabbits. J Physiol. 2012;590(5):1171–1180. doi: 10.1113/jphysiol.2011.218164. [PubMed] [Cross Ref]
17. Xie Y, Liao Z, Grandi E, Shiferaw Y, Bers DM. Slow [Na]i changes and positive feedback between membrane potential and [Ca]i underlie intermittent early afterdepolarizations and arrhythmias. Circ Arrhythm Electrophysiol. 2015;8(6):1472–1480. [PMC free article] [PubMed]
18. Qu Z, et al. Early afterdepolarizations in cardiac myocytes: beyond reduced repolarization reserve. Cardiovasc Res. 2013;99(1):6–15. doi: 10.1093/cvr/cvt104. [PMC free article] [PubMed] [Cross Ref]
19. Weiss JN, Garfinkel A, Karagueuzian HS, Chen PS, Qu Z. Early afterdepolarizations and cardiac arrhythmias. Heart Rhythm. 2010;7(12):1891–1899. doi: 10.1016/j.hrthm.2010.09.017. [PMC free article] [PubMed] [Cross Ref]
20. Guckenheimer, J., & Holmes, P. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. 462 (Springer-Verlag, New York, 1983).
21. Kuznetsov, Y. A. Elements of applied bifurcation theory (3rd Ed.) 632 (Springer-Verlag, New York, 2004).
22. Kurata Y, et al. Dynamical mechanisms of phase-2 early afterdepolarizations in human ventricular myocytes: insights from bifurcation analyses of two mathematical models. Am J Physiol Heart Circ Physiol. 2017;312(1):H106–H127. doi: 10.1152/ajpheart.00115.2016. [PubMed] [Cross Ref]
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):e1002061. doi: 10.1371/journal.pcbi.1002061. [PMC free article] [PubMed] [Cross Ref]
24. January CT, Riddle JM, Salata JJ. A model for early afterdepolarizations: induction with the Ca2+ channel agonist Bay K 8644. Circ Res. 1988;62(3):563–571. doi: 10.1161/01.RES.62.3.563. [PubMed] [Cross Ref]
25. January CT, Riddle JM. Early afterdepolarizations: mechanism of induction and block. A role for L-type Ca2+ current. Circ Res. 1989;64(5):977–990. doi: 10.1161/01.RES.64.5.977. [PubMed] [Cross Ref]
26. Zeng J, Rudy Y. Early afterdepolarizations in cardiac myocytes: mechanism and rate dependence. Biophys J. 1995;68(3):949–964. doi: 10.1016/S0006-3495(95)80271-7. [PubMed] [Cross Ref]
27. Guo D, et al. L-type calcium current reactivation contributes to arrhythmogenesis associated with action potential triangulation. J Cardiovasc Electrophysiol. 2007;18(2):196–203. doi: 10.1111/j.1540-8167.2006.00698.x. [PubMed] [Cross Ref]
28. Madhvani RV, et al. Targeting the late component of the cardiac L-type Ca2+ current to suppress early afterdepolarizations. J Gen Physiol. 2015;145(5):395–404. doi: 10.1085/jgp.201411288. [PMC free article] [PubMed] [Cross Ref]
29. Wang DY, Chae SW, Gong QY, Lee CO. Role of aiNa in positive force-frequency staircase in guinea pig papillary muscle. Am J Physiol. 1988;255(6Pt1):C798–C807. [PubMed]
30. Faber GM, Rudy Y. Action potential and contractility changes in [Na+]i overloaded cardiac myocytes: a simulation study. Biophys J. 2000;78(5):2392–2404. doi: 10.1016/S0006-3495(00)76783-X. [PubMed] [Cross Ref]
31. Bourgonje VJ, et al. Combined Na+/Ca2+ exchanger and L-type calcium channel block as a potential strategy to suppress arrhythmias and maintain ventricular function. Circ Arrhythm Electrophysiol. 2013;6(2):371–379. doi: 10.1161/CIRCEP.113.000322. [PubMed] [Cross Ref]
32. Shimizu W. Update of diagnosis and management of inherited cardiac arrhythmias. Circ J. 2013;77(12):2867–2872. doi: 10.1253/circj.CJ-13-1217. [PubMed] [Cross Ref]
33. Shimizu W, Horie M. Phenotypic manifestations of mutations in genes encoding subunits of cardiac potassium channels. Circ Res. 2011;109(1):97–109. doi: 10.1161/CIRCRESAHA.110.224600. [PubMed] [Cross Ref]
34. Sato D, et al. Synchronization of chaotic early afterdepolarizations in the genesis of cardiac arrhythmias. Proc Natl Acad Sci USA. 2009;106(9):2983–2988. doi: 10.1073/pnas.0809148106. [PubMed] [Cross Ref]
35. Huelsing DJ, Spitzer KW, Pollard AE. Electrotonic suppression of early afterdepolarizations in isolated rabbit Purkinje myocytes. Am J Physiol Heart Circ Physiol. 2000;279(1):H250–H259. [PubMed]
36. Xie Y, Sato D, Garfinkel A, Qu Z, Weiss JN. So little source, so much sink: requirements for afterdepolarizations to propagate in tissue. Biophys J. 2010;99(5):1408–1415. doi: 10.1016/j.bpj.2010.06.042. [PubMed] [Cross Ref]
37. Antzelevitch C, et al. The M cell: its contribution to the ECG and to normal and abnormal electrical function of the heart. J Cardiovasc Electrophysiol. 1999;10(8):1124–1152. doi: 10.1111/j.1540-8167.1999.tb00287.x. [PubMed] [Cross Ref]
38. Clancy CE, Rudy Y. Cellular consequences of HERG mutations in the long QT syndrome: precursors to sudden cardiac death. Cardiovasc Res. 2001;50(2):301–313. doi: 10.1016/S0008-6363(00)00293-5. [PubMed] [Cross Ref]
39. Nguyen TP, Xie Y, Garfinkel A, Qu Z, Weiss JN. Arrhythmogenic consequences of myofibroblast-myocyte coupling. Cardiovasc Res. 2012;93(2):242–251. doi: 10.1093/cvr/cvr292. [PMC free article] [PubMed] [Cross Ref]
40. Surawicz B. Electrophysiologic substrate of torsade de pointes: dispersion of repolarization or early afterdepolarizations? J Am Coll Cardiol. 1989;14(1):172–184. doi: 10.1016/0735-1097(89)90069-7. [PubMed] [Cross Ref]
41. Iyer V, Mazhari R, Winslow RL. A computational model of the human left-ventricular epicardial myocyte. Biophys J. 2004;87(3):1507–1525. doi: 10.1529/biophysj.104.043299. [PubMed] [Cross Ref]
42. Asakura K, et al. EAD and DAD mechanisms analyzed by developing a new human ventricular cell model. Prog Biophys Mol Biol. 2014;116(1):11–24. doi: 10.1016/j.pbiomolbio.2014.08.008. [PubMed] [Cross Ref]
43. Himeno Y, et al. A human ventricular myocyte model with a refined representation of excitation-contraction coupling. Biophys J. 2015;109(2):415–427. doi: 10.1016/j.bpj.2015.06.017. [PubMed] [Cross Ref]
44. Kurata, Y., Hisatome, I., Matsuda, H. & Shibamoto, T. Dynamical mechanisms of pacemaker generation in IK1-downregulated human ventricular myocytes: insights from bifurcation analyses of a mathematical model. Biophys J. 89(4), 2865–2887 (2005). [PubMed]
45. Kurata, Y., Matsuda, H., Hisatome, I. & Shibamoto, T. Regional difference in dynamical property of sinoatrial node pacemaking: role of Na+ channel current. Biophys J95(2), 951–977 (2008). [PubMed]
46. Kawakami H. Bifurcation of periodic responses in forced dynamic nonlinear circuits: Computation of bifurcation values of the system parameters. IEEE Trans Circuits and Systems. 1984;31(3):248–260. doi: 10.1109/TCS.1984.1085495. [Cross Ref]
47. Tsumoto K, Yoshinaga T, Iida H, Kawakami H, Aihara K. Bifurcations in a mathematical model for circadian oscillations of clock genes. J Theor Biol. 2006;239(1):101–122. doi: 10.1016/j.jtbi.2005.07.017. [PubMed] [Cross Ref]
48. Tsumoto K, Ueta T, Yoshinaga T, Kawakami H. Bifurcation analyses of nonlinear dynamical systems: From theory to numerical computations. Nonlinear Theory and Its Applications, IEICE. 2012;3(4):458–476. doi: 10.1587/nolta.3.458. [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group