Multiple Dynamical Mechanisms of Phase-2 Early Afterdepolarizations in a Human Ventricular Myocyte Model: Involvement of Spontaneous SR Ca2+ Release

Early afterdepolarization (EAD) is known to cause lethal ventricular arrhythmias in long QT syndrome (LQTS). In this study, dynamical mechanisms of EAD formation in human ventricular myocytes (HVMs) were investigated using the mathematical model developed by ten Tusscher and Panfilov (Am J Physiol Heart Circ Physiol 291, 2006). We explored how the rapid (IKr) and slow (IKs) components of delayed-rectifier K+ channel currents, L-type Ca2+ channel current (ICaL), Na+/Ca2+ exchanger current (INCX), and intracellular Ca2+ handling via the sarcoplasmic reticulum (SR) contribute to initiation, termination and modulation of phase-2 EADs during pacing in relation to bifurcation phenomena in non-paced model cells. Parameter-dependent dynamical behaviors of the non-paced model cell were determined by calculating stabilities of equilibrium points (EPs) and limit cycles, and bifurcation points to construct bifurcation diagrams. Action potentials (APs) and EADs during pacing were reproduced by numerical simulations for constructing phase diagrams of the paced model cell dynamics. Results are summarized as follows: (1) A modified version of the ten Tusscher-Panfilov model with accelerated ICaL inactivation could reproduce bradycardia-related EADs in LQTS type 2 and β-adrenergic stimulation-induced EADs in LQTS type 1. (2) Two types of EADs with different initiation mechanisms, ICaL reactivation–dependent and spontaneous SR Ca2+ release–mediated EADs, were detected. (3) Termination of EADs (AP repolarization) during pacing depended on the slow activation of IKs. (4) Spontaneous SR Ca2+ releases occurred at higher Ca2+ uptake rates, attributable to the instability of steady-state intracellular Ca2+ concentrations. Dynamical mechanisms of EAD formation and termination in the paced model cell are closely related to stability changes (bifurcations) in dynamical behaviors of the non-paced model cell, but they are model-dependent. Nevertheless, the modified ten Tusscher-Panfilov model would be useful for systematically investigating possible dynamical mechanisms of EAD-related arrhythmias in LQTS.


INTRODUCTION
Early afterdepolarization (EAD) is well known to trigger lethal ventricular arrhythmias, called Torsades de Pointes (TdP), in patients with long QT syndrome (LQTS) Shimizu and Horie, 2011;Shimizu, 2013). For prevention and treatment of ventricular arrhythmias in LQTS patients, therefore, elucidating the mechanisms of initiation and termination of EADs and how to suppress EADs is of crucial importance. There are many experimental studies regarding the mechanisms of EAD formation in cardiomyocytes, suggesting major contribution of reactivation of the L-type Ca 2+ channel current (I CaL ) to the initiation of EADs during the action potential (AP) phase 2 (e.g., January et al., 1988;January and Riddle, 1989;Weiss et al., 2010;Xie et al., 2010;Milberg et al., 2012a;Shimizu, 2013). However, recent experimental studies suggested the major role in EAD formation of the spontaneous Ca 2+ release from the sarcoplasmic reticulum (SR) (Volders et al., 2000;Choi et al., 2002;Zhao et al., 2012). In our recent theoretical study  using two human ventricular myocyte (HVM) models developed by Kurata et al. (2005) and O'Hara et al. (2011), referred to as K05 and O11 models, respectively, we could find EAD formations resulting from the I CaL reactivation, but not the spontaneous SR Ca 2+ release-mediated EADs. With respect to the termination of EADs (AP repolarization), theoretical studies Qu et al., 2013) using a guinea-pig ventricular myocyte model (Luo and Rudy, 1991) suggested the slowly activating delayed-rectifier K + channel current (I Ks ) as a key current to cause termination of EADs. However, our preceding study  suggested that the mechanisms of EAD termination were model-dependent, not necessarily requiring I Ks . Thus, despite many experimental and theoretical studies, how individual membrane and intracellular components contribute to the initiation, termination and modulation of EADs remains controversial.
The aims of this study were (1) to determine whether the ten Tusscher and Panfilov model (ten Tusscher and Panfilov, 2006; referred to as the TP06 model) for HVMs, which has often been used for simulations of reentrant arrhythmias in the human ventricle (ten Tusscher et al., 2007;Adeniran et al., 2012;Zimik et al., 2015;Kazbanov et al., 2016), could reproduce EAD formation in LQTS (validation of the model cell for EAD reproducibility), and (2) to define the contributions of individual sarcolemmal and intracellular components to the initiation, termination, and modulation of phase-2 EADs in the TP06 model in comparison with those in other HVM models (evaluation of model dependence for EAD mechanisms). As in our preceding study Tsumoto et al., 2017), we examined parameter-dependent changes in stabilities of steady states and AP dynamics in the HVM model from the aspect of bifurcation phenomena, which are parameterdependent qualitative changes in dynamical behaviors, in nonlinear dynamical systems (Guckenheimer and Holmes, 1983;Parker and Chua, 1989;Kuznetsov, 2003). Conditions and dynamical mechanisms of EAD formation in the paced model cell were determined in relation to bifurcations of the nonpaced model cell.
With respect to the dynamical mechanisms of EAD formation, we particularly focused on (1) whether and how contributions of each cellular component to occurrences of EADs and bifurcations in the TP06 model are different from those in the K05 and O11 models; (2) whether spontaneous SR Ca 2+ release-mediated EAD initiation, which did not occur in the K05 or O11 model, can be reproduced by the TP06 model in connection with a bifurcation (destabilization) of intracellular Ca 2+ dynamics; and (3) how slow I Ks activation, as well as I CaL inactivation and other slow factors, contributes to EAD termination. This study would further provide a theoretical background for experimental and simulation studies on mechanisms of EAD formation and EADtriggered reentrant arrhythmias in the LQTS human ventricle, as well as for prevention and treatments of life-threatening arrhythmias, like TdPs, in LQTS.

Base Mathematical Model
In this study, we tested the mid-myocardial (M) cell version of the TP06 model for HVMs (ten Tusscher and Panfilov, 2006), which could reproduce phase-2 EADs during inhibition of I Ks and/or the rapidly activating delayed rectifier K + channel current (I Kr ) or enhancement of I CaL . The M cell version was chosen because it has smaller I Kr and I Ks and thus more vulnerable to EAD formation than the epicardial or endocardial version, as suggested experimentally as well (Antzelevitch et al., 1999), and a few modifications were made for the M cell version of the TP06 model. Figure 1 shows simulated behaviors of APs, sarcolemmal ionic currents and intracellular Ca 2+ concentrations in the original and modified M cell versions of the TP06 model with various g Ks and g Kr values. Inconsistent with experiments for HVMs that observed only small prolongation of AP duration (APD) by I Ks inhibition (Jost et al., 2005;O'Hara and Rudy, 2012), the original version of the TP06 model, which has relatively large I Ks , exhibited marked APD prolongation during I Ks inhibition, and failed to reproduce greater APD prolongation and phase-2 EADs during I Kr inhibition (see Figure 1A). In addition, the Ca 2+ concentration in the SR (Ca SR ) (3-4 mM during 1-Hz pacing) was higher than the experimentally observed values of 1-2 mM for rabbit ventricular myocytes (Shannon et al., 2003(Shannon et al., , 2004. Therefore, the modified version, referred to as the "mTP06a" model, underwent the following modifications: (1) 60% reduction of the maximum I Ks conductance (g Ks ) with 50% increment of the maximum I Kr conductance (g Kr ) to reproduce the I Kr /I Ks inhibition experiments, and (2) 40% reduction in the SR Ca 2+ uptake rate (P up ) to reduce the Ca 2+ concentration in the SR during pacing under control conditions. These modifications yielded the experimentally observed small APD prolongation by I Ks inhibition and smaller Ca SR of 1.3-2.6 mM during pacing at 0.5-1 Hz, but not EAD formation ( Figure 1B). Therefore, we have developed another version of the modified TP06 model referred to as the "mTP06b" model with halved time FIGURE 1 | Simulated behaviors of APs (EADs), sarcolemmal ionic currents (I Ks , I Kr , I CaL , I NCX ) and intracellular Ca 2+ concentrations (Ca ss , Ca i , Ca SR ) in the M cell versions of the original TP06 (A), mTP06a (B), and mTP06b (C) models. To mimic the pathological conditions of LQT1 and LQT2, g Ks and g Kr values, respectively, were decreased by 30-100%; individual APs are labeled by the numbers representing the residual g Ks or g Kr (%Control), with ionic currents and Ca 2+ concentrations for each g Ks and g Kr value shown by the same colors. The horizontal dashed lines denote the 0 mV, zero current and zero concentration levels. Current amplitudes for the y-axis scale bars are given in pA/pF. The model cells were paced at 0.5 Hz, i.e., with the cycle length (CL) of 2 s for 60 min; AP waveforms, ion currents and Ca 2+ concentrations after the last stimulus are shown as steady-state behaviors under each condition. constant of I CaL inactivation (τ fL ) and doubled maximum I CaL conductance (g CaL ) on the basis of a previous theoretical study by Vandersickel et al. (2014) that required acceleration of the voltage-dependent inactivation of I CaL for reproducing EADs in the TP06 model. As shown in Figure 1C, the mTP06b model could reproduce the experimentally observed responses of HVMs to reductions of I Kr or I Ks , with EADs generated during I Kr reductions. Maximum conductance of the ionic channels, densities of transporters, and SR Ca 2+ uptake/release rates for the modified versions, as well as for the original version, are given in Supplementary Table S1.
The TP06 model for the normal activity of single HVMs is described as a non-linear dynamical system of 19 firstorder ordinary differential equations. The membrane current system includes the Na + channel current (I Na ), I CaL , I Kr , I Ks , 4aminopyridine-sensitive transient outward current (I to ), inwardrectifier K + channel current (I K1 ), background K + (I pK ), Na + (I bNa ) and Ca 2+ (I bCa ) currents, Na + -K + pump current (I NaK ), Na + /Ca 2+ exchanger current (I NCX ), and Ca 2+ pump current (I pCa ). Time-dependent changes in the membrane potential (V m ) are described by the equation, dV m /dt = I stim − (I Na + I CaL + I Kr + I Ks + I to + I K1 + I pK + I bNa + I bCa + I NaK + I NCX + I pCa ) (1) where I stim represents the stimulus current (in pA/pF).
The basic model systems include material balance expressions to define the temporal variations in concentrations of myoplasmic K + (K i ), Na + (Na i ) and Ca 2+ (Ca i ), and subspace Ca 2+ (Ca ss ), while external concentrations of K + , Na + and Ca 2+ were fixed at 5.4, 140, and 2.0 mM, respectively. For bifurcation analyses, K i was fixed at 140 mM for the removal of degeneracy (Krogh-Madsen et al., 2005;Kurata et al., 2008); effects of parameter-dependent changes in K i (∼5 mM) on EAD formation and bifurcation phenomena in the model cell were much smaller than those of the same amount of changes in Na i . Na i was unfixed unless otherwise stated, but fixed at 6 mM in some cases (e.g., for the slow-fast decomposition analysis and for voltage-clamped cells, as described later); changes in Na i during AP phase 2 and EAD formation in paced model cells were slow and relatively small.
Details on expressions, standard parameter values, and dynamics of the TP06 model are provided in the original article (ten Tusscher and Panfilov, 2006), and the original TP06 model has been implemented in a cellML-based open resource for public access 1 . In addition, the original TP06, mTP06a, and mTP06b models have been implemented in PhysioDesigner as XML-based Physiological Hierarchy Markup Language (PHML) 2 models. These models can be referred from PHML database (ID938 to 940) 3 , and simulations of their temporal behaviors can be performed using the software, Flint 4 .
Mutations of I Ks and I Kr channels (Roden et al., 1996;Chouabe et al., 1997;Anderson et al., 2006;Wiener et al., 2008;Kondo et al., 2016), as well as their pharmacological inhibitions (e.g., Carmeliet, 1992;Volders et al., 2003;Jost et al., 2005), are known to cause a wide range of channel conductance changes, leading to congenital or acquired LQTS type 1 (LQT1) and type 2 (LQT2), respectively. We developed LQT1-and LQT2-type model cells by continuously reducing g Ks and g Kr , respectively, from unity to zero. As illustrated in Figure 1, the TP06b model, but not the original version (with larger I Ks ) or the mTP06a model, reproduced phase-2 EADs (and AP repolarization failure) when g Kr became smaller as in LQT2 cardiomyocytes. In contrast, g Ks -reduced LQT1 model cells did not exhibit EADs but showed only slight prolongation of APDs under the basal condition, consistent with the recent experimental results from HVMs (Jost et al., 2005;O'Hara and Rudy, 2012).

Simulating Conditions of β-Adrenergic Stimulation
To simulate the condition of β-adrenergic stimulation (β-AS) as a major trigger of EADs and TdP in LQT1 patients, we modified the maximum conductance of ion channels and density of transporters based on previous reports (Zeng and Rudy, 1995;Volders et al., 2003;Kuzumoto et al., 2008), as described in our preceding article . g CaL and g Ks were increased up to 250% and 200%, respectively, according to previous reports for their changes during β-AS (Veldkamp et al., 2001;Saucerman et al., 2003;Himeno et al., 2008;Maltsev and Lakatta, 2010;Briston et al., 2014). Modifications of parameters for simulating β-AS are listed in Supplementary Table S2.

Basic Methods
Dynamic behaviors of the model cells were determined by numerically solving a set of non-linear ordinary differential equations including Eqn. 1. AP responses were elicited by 1-ms current stimuli of 60 pA/pF. When phase-2 EADs occurred at higher frequencies of the pacing, a complete AP repolarization was preceded by the next stimulus. Thus, the pacing cycle length (CL) was usually set to longer values of 2-5 s, except for analyses of the rate dependence. Numerical integration was performed by using MATLAB (The MathWorks, Inc., Natick, MA, United States) ODE solvers, ode15s and ode45, with the maximum relative error tolerance for the integration methods of 1 × 10 −8 .
Initial values of the state variables for computation at a parameter set were their steady-state values at a resting V m (see Supplementary Table S3 for the control conditions), which were perturbed by the current stimulus; the last values of the state variables in computation were used as initial conditions for the next computation at a new parameter set. The minimum V m during AP phase 4 (V min ) and the maximum V m during early phase 2 before emergence of an EAD (V max ), as well as APD at 90% repolarization (APD 90 ), were determined for individual APs or AP sets. Steady-state APs for the first parameter set were obtained by numerical integration for 30 min; subsequent numerical integration with each parameter set was continued until the differences in V min , V max and APD 90 between the newly calculated AP and the preceding one became <1 × 10 −3 of their preceding values.

Detection of EADs
EADs were detected as transient V m oscillations which emerged during late AP phase 2 (200 ms or later from the AP peak) and eventually led to AP repolarization to a resting V m . All the local minimum (EAD min ) and maximum (EAD max ) of V m oscillations during EAD formation, as well as a set of V min , V max and APD 90 , were determined for one AP cycle. When APs with EADs were irregular (arrhythmic), all the potential extrema (V min , V max , EAD min , and EAD max ) and APD 90 values were sampled for APs evoked by the last 10 stimuli.

Stability and Bifurcation Analyses for HVMs
We performed bifurcation analysis to explore how dynamical properties of the HVM model cell systems alter with changes in parameters. Detailed procedures for bifurcation analyses, i.e., locating equilibrium points (EPs) and limit cycles (LCs), detecting bifurcation points by determination of their stabilities, were provided in our previous articles (Kurata et al., 2008(Kurata et al., , 2012(Kurata et al., , 2013Tsumoto et al., 2017), as well as in textbooks (Guckenheimer and Holmes, 1983;Parker and Chua, 1989;Kuznetsov, 2003). In the present study, one-and twoparameter bifurcation diagrams for the non-paced cell model, as well as phase diagrams for the paced model cell, were constructed as functions of parameters, including (1) g Ks , g Kr , and g CaL , (2) scaling factor for I NCX , (3) P up , and (4) pacing CL. The maximum conductance of the ionic channel currents and P up were expressed as normalized values, i.e., ratios to the control values. Mechanisms of the initiation and termination of EADs were further examined by the slowfast decomposition analysis, in which stability and bifurcations of a fast subsystem are determined as functions of a slow variable, i.e., the gating variable xs for I Ks activation or Ca SR Qu et al., 2013;Xie et al., 2014). Basic concepts of bifurcation analysis, types of bifurcations, and methods for constructions of bifurcation/phase diagrams and slow-fast decomposition analysis are briefly described in Supplementary Materials.

Validation and Characterization of the mTP06 Models for LQTS HVMs
We first determined whether the mTP06a/b models can mimic the electrophysiological properties of I Kr -reduced LQTS type 2 (LQT2) and I Ks -reduced LQTS type 1 (LQT1) HVMs, in which EADs occur mainly at lower heart rates (bradycardia), and under β-AS, e.g., during exercise (tachycardia), respectively.

Decreases in I Kr and/or I Ks Accelerated EAD Formation in the mTP06 Model
The mTP06b model, but not the original TP06 or mTP06a model, exhibited an AP with EADs when I Kr was inhibited during 0.5-Hz pacing (Figure 1). Similarly, when g Kr was reduced by 40% during 0.2-Hz pacing, we could observe the AP with EADs in the g Kr -reduced mTP06b model, as shown in Figure 2A-a. This simulated AP with EADs was accompanied by oscillatory reactivation of I CaL , and was terminated (i.e., V m went back to the resting V m ) as I Ks increased (blue arrows for I Ks in Figure 2A-a). Further reducing g Kr by 60% caused another type of EADs with I Ks saturated before applying the second stimulus (Figure 2Ab); just before applying this second stimulus, the transient depolarization in plateau phase originating from the spontaneous SR Ca 2+ release and resulting activation of inward I NCX (red arrows in Figure 2A-b). In this case, repolarization did not occur without the next stimulus, i.e., repolarization failure occurred in the non-paced model cell after the cessation of pacing (Figure 2A-c). Figure 2B-i shows switching of AP dynamics when g Kr was gradually reduced during 0.2-Hz pacing. In this figure for the paced system, V m extrema of APs (V min /V max ) and EADs (EAD min /EAD max ), were plotted against g Kr . EADs emerged at g Kr = 0.772, i.e., with 22.8% block of I Kr (Figure 2B-i, top); the g Kr reduction led to increases in the number of EADs, resulting in the discrete increase of APD 90 values (see Figure 2B-i, bottom). The AP repolarization dynamics in the paced cell model relates to the dynamical behavior of the non-paced cell model because there is no stimulation during the AP repolarization. Therefore, we investigated the dynamical behavior of the non-paced cell model using bifurcation analysis. Figure 2B-ii shows one-parameter bifurcation diagrams as functions of g Kr , constructed for the non-paced mTP06b model (see also Supplementary Figure S1A showing those for the mTP06a model for comparison). In the non-paced mTP06b and mTP06a model, there existed three EPs as the steady states. The EP in the upper steady-state branch (V E3 in Figure 2B-ii and Supplementary Figure S1A-ii) was always unstable at positive g Kr values, while stable at negative g Kr values in the mTP06a model. When g Kr markedly reduced to a large negative value (out of range in Figure 2B-ii), the unstable EP (V E3 ) underwent the supercritical Hopf bifurcation (HB), which changed it to a stable EP and led to a generation of LC oscillation. The LCs spawned from the HB point were always unstable in the positive g Kr range (see gray lines in Figure 2B-ii and Supplementary Figure S1A-ii). On the one hand, we found small-amplitude spontaneous V m oscillations (SOs) that occurred at depolarized V m (red and blue lines in Figure 2B-ii, top) in the vicinity of the unstable LCs, as exemplified in Figure 2A-c. Just before the disappearance of SOs with increasing g Kr , the period of unstable LC markedly prolonged (see the gray zigzag trace in Figure 2B-ii, bottom). This marked prolongation of LC periods and the emergence of SOs correlated with very long APD (long-lasting EADs) and in the g Kr -reduced model cells paced at 0.2 Hz, illustrating two types of EADs. Temporal behaviors of the variables were computed for 30 min; the behaviors elicited by additional 3-4 stimuli are shown. APs with EADs were terminated (V m repolarized) by gradual increases in I Ks (a) or the next stimulus as indicated by the blue arrow (b). Spontaneous SR Ca 2+ releases and resulting increases of inward I NCX to provoke EADs occurred at g Kr = 0.4, as indicated by the red arrows (b). With the smaller g Kr , sustained V m oscillations driven by spontaneous SR Ca 2+ releases were observed after cessation of pacing, i.e., in the non-paced model cell (c). (B) Potential extrema of APs and EADs, and the AP duration (APD) measured at 90% repolarization (APD 90 ) are plotted as functions of normalized g Kr for the paced model cell (i). One-parameter bifurcation diagrams depicting steady-state V m at equilibrium points (EPs), the extrema of limit cycles (LCs) and spontaneous oscillations (SOs), and the periods of LCs and SOs plotted against g Kr are also shown for the non-paced model cell (ii). In the panel (i), AP dynamics during 0.2-Hz pacing were computed for 1 min at each g Kr value, which was reduced from 1.0 to −0.2 at an interval of 0.001. The minimum V m during AP phase 4 (V min ) and the maximum V m during AP phase 2 before EAD formation (V max ) are indicated by black dots for rhythmic APs, and by orange (V min ) and light green (V max ) dots for arrhythmic APs. When EADs appeared, their local potential minimum (EAD min ) and maximum (EAD max ) were plotted by blue and red dots, respectively. In the panel for APD 90 , the black, blue, and magenta dots represent APD 90 values for regular APs without EAD, regular APs with EADs, and arrhythmic APs with EADs, respectively (no-EAD: APs without EAD, + EAD: APs with EADs). The points labeled as "a" and "b" indicate g Kr values for which AP behaviors are shown in Panel (A). In the panel (ii), the steady-state branches as loci of V m at EPs (V E1-3 ), periodic branches as the potential minimum (LC min ) and maximum (LC max ) of LCs, and potential extrema of SOs (SO min , SO max ), as well as the periods of LCs and SOs, are plotted for the non-paced model cell. The steady-state branch V E1 is stable (green solid lines), while V E2 and V E3 unstable (black dashed lines). The periodic branches represented by gray solid lines are always unstable. The point labeled as "c" indicates the g Kr values for which SOs are shown in A-c.
irregularity of the repolarization time in the paced cell model (compare the dotted ranges in Figures 2B-i,ii).
In contrast, V E3 in the I Ks -eliminated (g Ks = 0) non-paced model cells was stabilized via an occurrence of the subcritical HB when g Kr was reduced (solid green lines to the left of the label "H" in Supplementary Figure S1B-ii). Thus, I Ks inhibition caused drastic shift of HB points toward higher g Kr values. Unstable LCs emerged via the subcritical HB, not changing their stability in the g Kr range tested. In the I Ks -eliminated paced mTP06a/b models (Supplementary Figure S1B-i), decreasing g Kr did not yield EADs, but abruptly changed APs without EADs to local responses in the depolarized V m range during pacing, i.e., arrest at stable EPs (V E3 ) without pacing.
To evaluate the dependencies of EAD formation on g Kr and g Ks , we performed AP simulations using the mTP06b model with various sets of g Kr and g Ks . Figure 3A-i shows a phase diagram of AP behaviors for changes in g Ks and g Kr values with 0.2-Hz pacing. By characterizing AP behaviors observed in the paced mTP06b model, the g Ks -g Kr parameter plane was divided into three regions: (1) AP without EAD, (2) AP with EADs (colored regions; see examples of Figure 3B for the points "d" and "e" in Figure 3A-i), and (3) local response (dotted region; FIGURE 3 | I Kr /I Ks -dependent EAD formation, dynamics and bifurcations of the mTP06b model. (A) A phase diagram indicating the region of EAD formation (and local responses) in the paced model cell (i) and two-parameter bifurcation diagrams for the non-paced model cell (ii) on the g Ks -g Kr parameter plane. In the diagram for the paced model cell (i), the thick red solid, black dashed and thin black solid lines, respectively, indicate parameter sets of critical points at which short-term EADs (APD 90 < 5 s, as in Figure 2A-a and [B-(i)], long-term or sustained EADs (APD 90 > 5 s, as in Figure 2A-b and [B-(ii)], and a local response (C) emerged; parameter regions in which short-term EADs, long-term or sustained EADs, and local responses occur are shown as the light-gray region labeled as "fR" (fast repolarization), blue region labeled as "RF" (repolarization failure), and dotted region, respectively. In the two-parameter bifurcation diagram for the non-paced model cell (ii), H, SO 1 , and SO 2 indicate parameter sets of HB points, critical points at which SOs emerged, and critical points at which SOs switched into quiescence, respectively. The parameter regions in which SOs and convergence to the steady state (V E3 ), i.e., arrest, can be observed are indicated as the shaded and orange regions, respectively. The labels "sEP" and "uEP" indicate the areas of stable and unstable EPs, respectively, divided by the HB curve. The panel (iii) is the diagram for which the phase diagram (i) is superimposed upon the two-parameter bifurcation diagram (ii). The points labeled as "N", "LQT1" and "LQT2" denote the normal, LQT1, and LQT2 conditions, respectively. The points "a"-"f" indicate parameter sets for which AP behaviors are shown in Figures 2A -a-c, and [B (d,e)] [C (f)]. (B) Representative behaviors of APs with EADs during 0.2-Hz pacing at the points labeled as "d" and "e" in (A), which are classified into the fast repolarization type (d) and repolarization failure type (e) behaviors, respectively. (C) An example of the local response during 0.2-Hz pacing at the point "f" in (A).
see an example of Figure 3C for the point "f " in Figure 3Ai). We further separated the region of the AP with EADs into two regions based on characteristics of the repolarization time in an AP with EADs: During 0.2-Hz pacing, further decreases in g Kr (and/or g Ks ) in the EAD region altered an AP with shorter APD 90 of ≤5 s that repolarizes before the next stimulus to an AP with longer APD 90 of >5 s that is repolarized by the next stimulus, as shown in Figure 2A; then, the APs with EADs were defined as "fast repolarization (fR)" type for the former and "repolarization failure (RF)" type for the latter, which are exemplified in Figures 2A-a,b, respectively. The fR and RF types were distinguished by AP behaviors after an extra stimulus following the last test stimulus to cause AP repolarization, as illustrated in Figure 3B: The fR-type AP repolarized to resting V m within 5 s ( Figure 3B-i), while the RF-type one did not ( Figure 3B-ii); in this case, APD 90 values of the RF-type AP were almost always more than 10 min. Decreases in g Kr and/or g Ks required for EAD formation were much smaller in the mTP06b model than in the mTP06a model (compare Figure 3A-i and Supplementary Figure S2). The borderline of EAD initiation (the red solid line in Figure 3Ai) shifted in a g Ks -dependent manner, with the g Kr region of EADs broadening as g Ks increased. Furthermore, two-parameter bifurcation analysis for the non-paced cell model (Figure 3A-ii) determined three areas with different behaviors: (1) quiescence at a stable EP (resting state; V E1 ) with no stable EP or LC at depolarized V m , (2) co-existence of quiescence at a stable EP (resting state) and a stable LC or SO at depolarized V m (shaded area labeled as "Spontaneous Oscillation"), and (3) co-existence of two stable EPs at V E1 and depolarized V m (V E3 ), i.e., the arrest at depolarized V m (colored area labeled as "Arrest (stable V E3 )").
To clarify relationships between AP responses observed in the paced cell model and bifurcations occurred in the non-paced cell model, we superimposed the phase diagram on the twoparameter bifurcation diagram (Figure 3A-iii). Most of the SO region in which SOs can be observed in the non-paced cell model was included in the RF region, suggesting the relation of spontaneous SR Ca 2+ release-mediated sustained EADs to SOs (Figure 3B-ii). The borderline between local response and AP with EADs corresponded to the HB set in the non-paced cell model, indicating that V m in the paced cell model converges to the stable EP (V E3 ) in the area of local response.

Slow and Rapid Pacing Facilitated EAD Formation in the mTP06b Model
To further validate the mTP06b model as a LQT2 model, we next determined whether EAD formation in the g Kr -reduced mTP06b model is facilitated at lower pacing rates (in bradycardia). Rate effects on EAD formation are shown in the diagrams depicting the g Kr regions of EADs as functions of the pacing cycle length ( Figure 4A). EAD formation in the g Kr -reduced system was promoted at lower pacing rates in the Na i -variable system, while prevented in the Na i -fixed system. As in the K05 and O11 models , the facilitation of EAD formation at lower pacing rates in the Na i -variable mTP06b model was accompanied by the decrease in Na i , which resulted in the decrease of outward I NaK leading to delays in AP repolarization and EAD formation ( Figure 4B-i). In the Na i -fixed mTP06b model, the inhibition of EAD formation at lower pacing rates accompanied marked outward shift of I NCX resulting from diminished Ca i transients ( Figure 4B-ii). In Supplementary Figure S3, two-parameter bifurcation diagrams on the g Ks -g Kr parameter plane are also shown for the Na i -variable and Na i -fixed mTP06b model cells paced at 0.2 and 1 Hz. In the Na i -variable system (Supplementary Figure S3A), slower pacing promoted EAD formation during decreases of g Kr and/or g Ks and broadened the parameter region of EADs; in the Na i -fixed system (Supplementary Figure  S3B), however, the rate-dependent changes in the onset and region of EADs were opposite to those in the Na i -variable system (compare the gray and blue areas in each panel of Supplementary Figure S3).
The rate dependence shown in Figure 4 and Supplementary Figure S3 was determined by quasi steady-state dynamics for each parameter set. In LQT2 patients, however, EADs and TdP may often be induced by abrupt pause or transient slowing of heartbeats (bradycardia) at rest or during sleep. Thus, we also determined how EADs emerge after sudden reductions of pacing rates (Supplementary Figure S4). When a pacing CL was increased from 1 s to 3, 4, and 5 s in the I Kr -reduced mTP06b model (g Kr = 0.721), EADs were first induced by 172nd, 74th, and 51st stimulus, respectively, after the reductions in pacing rates; pause-induced EAD or early onset of EADs after the increment in pacing CL was not observed, but long bradycardiac periods of more than 4 min were needed for EAD formation in this model cell.
In the mTP06b model, rapid pacing (CL < 1 s) also facilitated EAD formation via increases in Ca i and Ca SR , and resulting increases of inward I NCX (data not shown). EADs are known to be inhibited at higher pacing rates by accumulation of slowly deactivating I Ks as well as reductions in I CaL due to slow recovery.
Cumulative I Ks increments and I CaL reductions were certainly observed at the higher pacing rates in the mTP06b model as well; however, the enhanced inward I NCX appeared to cause APD prolongation and EAD formation in this model cell.

Spontaneous SR Ca 2+ Release-Medicated EADs Occurred During β-AS
To validate the mTP06b model as a LQT1 model, i.e., to determine whether the I Ks -reduced model cell can exhibit EADs under the conditions of β-AS, we examined susceptibilities to EAD generation during β-AS of the normal and LQT1 versions of the mTP06b model. For the LQT1 model cell, g Ks was reduced by 50% and 75%, following the reports for the KCNQ1 mutations M437V and A590W, respectively (see Sogo et al., 2016). Figure 5A shows simulated APs of the normal and LQT1 versions of the mTP06b model under the basal condition and conditions of β-AS with g CaL increased to 140, 150, 160, and 180% of the control value. The LQT1 model cells exhibited longer APDs under the basal condition (APD 90 of 334 ms with the normal g Ks vs. 348 ms with 50% g Ks and 356 ms with 25% g Ks ) and EADs under β-AS with g CaL increased by 60% or more for 50% g Ks and 40% or more for 25% g Ks , whereas the normal cell did not exhibit EAD. By constructing a twoparameter bifurcation diagram on the g Ks -g CaL plane for the Na i -variable model cell paced at 1 Hz (Figure 5B), we could explain their EAD formation under the conditions of β-AS. As in the K05 model , EAD formation during g CaL increases under β-AS could be inhibited by concomitant g Ks increases more effectively in the normal mTP06b model than in the LQT1 models: The LQT1 model cells entered the area of EAD formation with smaller increases in g CaL (50.8% or more for the 50% g Ks reduction and 31.0% or more for the 75% g Ks reduction), while the normal cell with more than 87.7% increases in g CaL . Thus, the mTP06b model could recapitulate EAD formation via enhancement of I CaL during β-AS in the LQT1 cardiomyocyte. Under β-AS with higher g CaL and P up , spontaneous SR Ca 2+ releases as evidenced by abrupt falls in Ca SR without I CaL reactivation often occurred, leading to Ca i elevations (see Supplementary Figure S5), increments of inward I NCX , and resultant EADs (or spontaneous V m oscillations), as indicated by the dots in Figure 5A.

Influences of SR Ca 2+ Cycling, I NCX and I CaL on EAD Formation
Following the finding of spontaneous SR Ca 2+ releases which occurred especially under β-AS with enhanced I CaL and SR Ca 2+ uptake, we next examined how SR Ca 2+ uptake/release (intracellular Ca 2+ dynamics) and I NCX regulated by the intracellular Ca 2+ , as well as I CaL regulated by the subspace Ca 2+ , affect EAD formation and bifurcations of dynamical behaviors in the mTP06b model by changing P up , the scaling factor for I NCX , or g CaL . Figure 6A shows phase diagrams on the P upg Kr and P up -g CaL parameter planes. P up values were varied from zero to 2-times the control value, assuming the effects of SR Ca 2+ pump inhibitors and β-AS (Maltsev and Lakatta, 2010;Briston et al., 2014). The region of EADs shrank with reducing P up as in the K05 and O11 models , while   Supplementary  Table S2). The LQT1 model cell was assumed to have reduced g Ks of 50% or 25% of the control value. The points of the control (basal) conditions for cardiomyocytes with the normal and reduced g Ks are labeled as "N" and "LQT1", respectively. The arrows indicate the parameter shifts from the basal condition to the conditions of β-AS with g Ks doubled and g CaL increased to 140, 150, 160 and 180% of the control value.
broadening at higher P up ; however, for the emergence of EADs (red solid lines in Figure 6A), the critical g Kr value was not decreased but slightly increased (the critical g CaL value was not increased but slightly decreased) as P up reduced. The facilitated EAD formation at smaller P up was associated with increased Ca i and resulting inward shift in I NCX as well as slight increases in I CaL during AP late phase 2 (see Supplementary Figure  S6A). The two-parameter bifurcation analysis offered further information on how the region of EADs depends on P up and g Kr . As shown in Supplementary Figure S7, the critical set of the emergence of EADs and the HB set were mostly parallel to the P up and g Kr axes, respectively, suggesting that alterations in P up contributed not to EAD formation but to rather stability changes of EP (V E3 ) in the non-paced mTP06b model; HB points disappeared with the emergence of spontaneous V m and Ca 2+ oscillations at higher P up , indicating that the SR Ca 2+ uptake/release machinery destabilizes EPs and thereby induces spontaneous oscillations. When g Kr was markedly reduced in the mTP06b model, spontaneous SR Ca 2+ releases to cause transient increases in intracellular Ca 2+ concentrations (Ca ss and Ca i ) and resulting activation of inward I NCX occasionally occurred with prolonged APD (Figure 6B-a). Increasing P up shortened the time to the emergence of the first spontaneous Ca 2+ release and raised the incidence and frequency of spontaneous Ca 2+ oscillations to yield EADs or V m oscillations (Figure 6B-b).
Contributions of I NCX to bifurcations and EAD formation were also explored in relation to those of intracellular Ca 2+ dynamics, SR Ca 2+ cycling, and I CaL . The scaling factor of I NCX was varied from 0.1 to 10, within the range of experimental changes in Na + /Ca 2+ exchanger densities or I NCX (Milberg et al., 2008(Milberg et al., , 2012bPott et al., 2012). On the I NCX -g Kr and I NCX -g CaL parameter planes (Supplementary Figure S8), enhancement of I NCX yielded the upward shift in the critical g Kr and downward shift in the critical g CaL for EAD formation (see red curves in Supplementary Figure S8). However, the TP06b model did not exhibit a significant shift in the critical g Kr or g CaL for EAD formation when I NCX was reduced; only small (20-30%) inhibition of I NCX was effective in shifting the critical points FIGURE 6 | Influences of SR Ca 2+ handling on EAD generation in the mTP06b model. (A) Phase diagrams on the P up -g Kr (i) and P up -g CaL (ii) parameter planes, depicting displacements of critical points for the occurrence of short-term EADs (red solid lines) and long-term or sustained EADs (dashed lines) as well as the emergence of local responses (black solid lines). By the parameter sets of these critical points, the parameter planes are divided into the areas of APs with short-term EADs (fR), APs with long-term or sustained EADs (RF) and local response, as described for Figure 3A. toward the prevention of EADs, with further inhibition resulting in the promotion of EADs. Whether EADs emerge or not depended mainly on the amplitude of inward I NCX and I CaL during the AP late phase 2: As exemplified in Supplementary Figure S6B, disappearance of EADs with lower I NCX density was accompanied by a decrease of inward I NCX and a slight reduction of I CaL with increased inactivation during the preconditioning phase just before initiation of the first EAD (see the ellipses and inset in Supplementary Figure S6B).
We finally examined effects of I CaL on EAD formation in the paced mTP06b model and bifurcations of dynamical behaviors in the non-paced mTP06b model by changing g CaL . g CaL -dependent changes in AP dynamics observed in the paced model cell when I Ks was normal (g Ks = 1) and one-parameter bifurcation diagrams as functions of g CaL for the non-paced cell model are shown in Supplementary Figures S9A-i,ii, respectively. The one-parameter bifurcation diagrams for g CaL (Supplementary Figure S9A-ii) suggest the scenario of EAD formation during enhancement of I CaL , which is different from those in the K05 and O11 models : Increments of g CaL yielded unstable EPs via a saddle-node bifurcation (SNB) of EPs and unstable LCs via a SNB of LCs. With normal I Ks (g Ks = 1), an enhanced g CaL of 1.298-fold the control value was high enough for EAD formation in the mTP06b model (Supplementary Figure S9A-i), whereas unrealistically large increases in g CaL (to 4.248-fold the control value) were required in the mTP06a model (Supplementary Figure S10A). Supplementary Figure S9B shows a phase diagram of AP behaviors in the paced model cell (Supplementary Figure S9Bi) and a two-parameter bifurcation diagram for the non-paced model cell (Supplementary Figure S9B-ii), as well as the merged diagram (Supplementary Figure S9B-iii), on the g CaLg Ks parameter planes. Decreasing g Ks shifted the critical g CaL value for EAD generation toward lower values and enlarged the g CaL region of EADs (RF) in the mTP06b model. Larger g CaL (going into the RF region in Supplementary Figures S9B-i,iii) led to the AP behavior classified into the RF type with smallamplitude spontaneous V m oscillations around unstable LCs. EPs (V E3 ) in the mTP06 models were unstable independent of g CaL unless g Ks was extremely low or high; no HB occurred for moderate variation of g Ks value and consequently stability changes of the EP did not occur (see also Supplementary  Figure S9B-ii). In the I Ks -eliminated mTP06a/b models (g Ks = 0), an EP (V E3 ) was stabilized via supercritical HBs at relatively small g CaL ; stable LCs emerging from the HBs were immediately destabilized via a period-doubling bifurcation (PDB) or Neimark-Sacker bifurcation (NSB) (Supplementary Figure S10B). EAD did not occur at g Ks = 0; larger I CaL caused repolarization failure, in this case, local response.

I Ks Activation-Dependent Bifurcations of the Fast Subsystem Associated With EAD Formation
To clarify the dynamical mechanisms of EAD formation in the I Kr -reduced LQT2-type mTP06b model and why EADs emerge at larger g Kr in the mTP06b model than in the mTP06a model (compare Figure 2 and Supplementary Figure S1), we further performed the slow-fast decomposition analysis Qu et al., 2013;Xie et al., 2014). The I Ks activation gating variable xs or I Ks channel open probability (xs 2 ) appears to be a slow variable yielding the termination of EADs (Figure 2Aa, the second from top). Thus, bifurcation diagrams for the fast subsystem composed of the state variables other than the slow variables xs, Na i and Ca SR were first constructed as functions of xs 2 , with Na i and Ca SR fixed at constant values (Figure 7, left); then, trajectories of the full system (with fixed Na i and Ca SR ) were superimposed on the diagrams (Figure 7, right). The quasi-EP (qEP), defined as a steady state of the fast subsystem, at depolarized quasi-V m (qV E3 ) has possessed a property of spiral sink in the mTP06b model ( Figure 7A) but spiral source in the mTP06a model ( Figure 7B) at xs 2 = 0. Stable qEP in the former was destabilized via an HB as xs 2 increased. The g Kr reduction led to broadening of the xs 2 region of stable qEPs (compare green traces of qV E3 in Figures 7A-i,ii, left). The g Kr reduction-induced broadening of the xs 2 range of stable qEPs yielded a transient trapping of the full system trajectory in the attractor basin of the stable qEP ( Figure 7A-ii, right). This trapping of the full system trajectory around the stable qEP as spiral sink sustained until the trajectory came across the steady-state xs 2 curve. This trapping phenomenon was not observed in the mTP06a model ( Figure 7B, right) or the I Kr -normal mTP06b model (Figure 7Ai, right), because the full system trajectories did not intersect with the stable steady-state branch (qV E3 ) before intersecting the steady-state xs 2 curve. These results indicate that an acceleration of the voltage-dependent I CaL inactivation to form the mTP06b model from the mTP06a model plays a critical role in the stabilization of qV E3 , consequently leading to the trapping of the full system trajectory.

Release-Mediated EAD
To clarify the dynamical mechanisms of spontaneous SR Ca 2+ release-mediated EAD formation, we further examined the stability, dynamics and bifurcations of the voltage-clamped mTP06 model. Ca 2+ dynamics during a train of 1-s depolarizing test pulses to −10 mV (from the holding potential of −85 mV) applied at 2-s intervals to mimic APs evoked by 0.5 Hz pacing were first determined for the mTP06b model with different P up (Figure 8A). Spontaneous SR Ca 2+ releases occurred when Ca SR increased at higher P up , as indicated by the dots in Figure 8A; as P up increased, the time to the first Ca 2+ release and period of spontaneous Ca 2+ releases shortened, and their frequency increased. Figure 8B shows one-parameter bifurcation diagrams of the steady-state stability and dynamics of Ca i as functions of the clamped-V m in the voltage-clamped mTP06b model. Steadystate intracellular Ca 2+ concentrations (EPs) in the voltageclamped model cell were stable at hyperpolarized and depolarized V m (green traces in the right and middle panels of Figure 8B) but became unstable via supercritical HBs in the V m range of AP phase 2 and early phase 3 (dashed traces in Figure 8B, middle). LCs emerging from the HB points were first stable but were destabilized via NSBs after small changes in V m ; spontaneous FIGURE 7 | Dynamical mechanisms of EAD initiation and termination determined by the slow-fast decomposition analysis for the mTP06 models. Shown are one-parameter bifurcation diagrams of quasi-equilibrium points (qEPs) and quasi-limit cycles (qLCs), where the steady-state branches as loci of V m at qEPs (qV E1-3 ) and periodic branches as the potential minimum (qLC min ) and maximum (qLC max ) of qLCs are depicted as functions of the square of the I Ks activation gating variable (xs 2 ), i.e., I Ks channel open probability for the fast subsystems of the g Kr -normal [A-(i), left] and g Kr -reduced [A-(ii), left] mTP06b model and g Kr -reduced mTP06a model (B, left). Other slow variables, Na i and Ca SR , were fixed at constant values: Na i = 6 mM for all cases; Ca SR was fixed at the value which was reached just before occurrence of the first EAD or the maximum values during AP phase 2 (when no EAD occurred), i.e., at 0.5 mM and 1.5 mM for the normal and g Kr -reduced mTP06b model, respectively, and at 0.5 mM for the g Kr -reduced mTP06a model. The steady-state branches consist of the stable (green solid lines) and unstable (black dashed lines) segments. The periodic branches (gray solid lines) are all unstable. The blue lines indicate the steady-state xs 2 curve. Trajectories of the full system (with the fixed Ca SR and Na i ) are superimposed on the bifurcation diagrams for the fast subsystems (red lines in each right panel). The arrows indicate the directions of changes in the state variables. H, Hopf bifurcation; hom, homoclinic bifurcation.

FIGURE 8 | Stability and bifurcations of intracellular Ca 2+ dynamics in voltage-clamped mTP06b model cell. (A)
Dynamics of intracellular Ca 2+ concentrations as well as Ca 2+ -dependent sarcolemmal currents in the Na i -fixed voltage-clamped model cell under the control condition (P up = 1) and the conditions of β-AS (P up = 1.41, 1.67, and 2.00). Temporal behaviors of the voltage-clamped model cell during a train of 1-s step depolarization from −85 mV to −10 mV at 0.5 Hz were computed for 10 min; Ca i , Ca ss , Ca SR , I CaL and I NCX for the last 12 s (6 pulses) are shown as steady-state dynamics. Spontaneous SR Ca 2+ releases, as evidenced by abrupt falls of Ca SR and increases in Ca i , Ca ss and inward I NCX with attenuated I CaL , occurred at higher P up (indicated by the dots). (B) One-parameter bifurcation diagrams of the equilibrium point (EP) and extrema of limit cycles (LC min/max ) and spontaneous Ca i oscillations (CaO min/max ) as functions of V m for the voltage-clamped model cell (left and middle). The middle panel shows an enlarged diagram of the rectangular area in the left panel. The periods of spontaneous Ca 2+ oscillations (CaO) and limit cycles (LC) are also plotted against V m (right). H 1-2 , Hopf bifurcations of the EP; NS 1-2 , Neimark-Sacker bifurcations of the LC. (C) Two-parameter diagrams on the V m -P up (left), V m -I NCX (middle) and V m -g CaL (right) planes, indicating how the unstable V m range changed depending on P up , I NCX and I CaL . HB values, i.e., the critical V m at which an EP is (de)stabilized are plotted as functions of P up , I NCX and I CaL for the Na i -fixed mTP06b model; the HB points were very close to the Neimark-Sacker bifurcation points at which LCs were destabilized with the emergence of CaOs.
Ca 2+ oscillations occurred in the V m range of unstable LCs, i.e., between NS 1 and NS 2 (gray traces labeled as LC min and LC max for the minimum and maximum Ca i during LC oscillations in Figure 8B). As shown in Figure 8C, the unstable V m region (uEP) was enlarged by increasing P up (see Figure 8C, left), decreasing I NCX activity (Figure 8C, middle), and/or enhancing I CaL (Figure 8C, right), all of which led to increases in Ca SR .
To further clarify the Ca SR -dependent mechanism of spontaneous SR Ca 2+ releases in the P up -increased mTP06b model and why SR Ca 2+ release-mediated EADs emerge more frequently at larger P up , we also performed the slow-fast decomposition analysis for the slow variable Ca SR . Bifurcation diagrams were constructed as functions of Ca SR for the voltageclamped fast subsystem composed of the voltage-independent state variables f CaL (Ca 2+ -dependent inactivation gate for I CaL ), R (proportion of closed SR Ca 2+ release channels), Ca ss , and Ca i (Figure 9). Trajectories of the voltage-clamped full system dynamics as shown in Figure 8A for the normal (1) and enhanced (1.67) P up were superimposed on the diagrams. The steady states of the fast subsystem, stable at lower Ca SR (green traces in the middle and right panels of Figure 9), became unstable via an HB at higher Ca SR (dashed traces in the middle and right panels of Figure 9). In the P up -enhanced system, spontaneous SR Ca 2+ releases as shown in blue trajectories in the middle and right panels of Figure 9B occurred when the full system trajectory, moving along the stable steady-state branch, passed through the HB point, i.e., when Ca SR exceeded the HB value. In contrast, the P up -normal system did not exhibit spontaneous SR Ca 2+ release, because an increment of Ca SR (Ca 2+ refilling of the SR) during Ca 2+ transient decay was too slow for the full system trajectory to reach the HB point for Ca SR before V m repolarization (Figure 9A, middle and right).

DISCUSSION
In this study, we theoretically investigated dynamical mechanisms of EAD formation in the TP06 model for HVMs, which has often been used for simulations and theoretical analyses of reentrant arrhythmias, automaticity, multi-stability and EAD formation in HVMs, in relation to the model cell dynamical behaviors and their bifurcations. In summary, EAD formation and its dynamics in the paced (non-autonomous) mTP06 model cell basically depended on stability and bifurcations of the non-paced (autonomous) model cell. Bifurcation phenomena and dynamical mechanisms of EAD formation in the mTP06 model were different from those in the K05 and O11 models tested previously  in several respects (see also Supplementary Materials for additional discussions).

Validation of the mTP06 Model for EAD Reproducibility in LQT1 and LQT2 Conditions
EAD Formation in LQT1 and LQT2 Conditions (mTP06a vs. mTP06b) Like the K05 model, the TP06b model with accelerated I CaL inactivation could recapitulate EAD formation in the I Ks -reduced LQT1-type and I Kr -reduced LQT2-type HVMs. The mTP06b model was much more vulnerable to EAD formation than the mTP06a model, consistent with the previous experimental finding that slowing I CaL inactivation eliminated EADs (Qu et al., 2013). As demonstrated by the slow-fast decomposition analysis (Figure 7), higher susceptibility of the mTP06b model to EAD development is attributable to the stabilization of qEPs at depolarized V m close to the plateau V m in the xs-parameterized fast subsystem by accelerating I CaL inactivation.

Rate Dependence of EAD Formation (Validation for LQT2 Model)
In LQT2 patients, fatal cardiac events often occur during sleep or at rest, i.e., in bradycardia (Shimizu and Horie, 2011;Shimizu, 2013). The K05 and O11 models could partially reproduce the bradycardia-related EADs . Like the other HVM models, the Na i -variable mTP06b model could partly reproduce the rate-dependent EAD generation in FIGURE 9 | Slow-fast decomposition analysis for the slow variable Ca SR of the voltage-clamped mTP06b model. Normal (middle) and expanded (right) scale views of one-parameter bifurcation diagrams, where the steady-state branch as a locus of Ca ss at EPs, as well as periodic branches as the minimum (LC min ) and maximum (LC max ) of Ca ss limit cycles (LCs) that were almost always unstable, are plotted as functions of Ca SR for the fast subsystems of the voltage-clamped model cells with the normal (A) and increased (B) P up . The green solid and black dashed lines represent stable and unstable steady-state branch, respectively. Subcritical HB (H) with the emergence of an unstable LC, NSB from a stable LC to an unstable one following a SNB of LCs (NS), and homoclinic bifurcation (hom) are located. Trajectories projected onto the Ca SR -Ca ss plane of the full system during a 1-s step depolarization to −10 mV (from −85 mV) and subsequent 1-s repolarization to −85 mV are superimposed on the diagrams; Ca ss and Ca SR dynamics in the full system during this voltage-clamp pulse are shown on the left. The labels "1"-"9" for the arrows on each diagram (middle, right) indicate the time course of changes in the state variables (Ca ss and Ca SR ) along the trajectories, with the labels also given for the Ca ss and Ca SR dynamics (left).
LQT2 patients (Figure 4). In the I Kr -reduced mTP06b model, however, EADs appeared only when pacing CLs increased to 3 s or more (i.e., pacing rates decreased to 20 beats/min or less) and extreme bradycardia continued for more than 4 min (Supplementary Figure S4). Although LQTS patients are known to exhibit sinus arrest or severe bradycardia due to coexisting sick sinus syndrome or atrio-ventricular block (Roden et al., 1996;Chiang and Roden, 2000), such long-lasting extreme bradycardia may be unlikely to occur often in LQT2 patients. Like the Luo-Rudy model for guinea-pig ventricular myocytes (Viswanathan and Rudy, 1999), the I Kr -reduced O11 model could reproduce pause-induced EAD formation on increasing a pacing CL from 1 s to 2 s, which was attributable to a decrease in I Ks and Ca i increase-mediated enhancement of inward I NCX at the lower pacing rate (data not shown). In contrast, the mTP06 model did not exhibit EADs by a single pause or transient bradycardia, which may be a limitation of this model cell.
In our previous study for the Na i -variable K05 and O11 models , the facilitation of EAD formation during lower rate pacing was accompanied by the decrease in Na i and resulting reductions in outward I NaK and inward shift of I NCX . This study also demonstrated for the mTP06b model that the facilitation of EAD formation at lower pacing rates was mainly due to the decrease in Na i and resultant changes in I NaK ( Figure 4B). Thus, the major mechanism for bradycardia-related EADs in the mTP06b model is essentially the same as that in the K05 model. Bradycardia-induced EADs are believed to be ascribable to a reduction of I Ks (and increment of I CaL ). In the mTP06 model, however, I Ks reduction (or I CaL increment) did not occur when a pacing CL increased from 1 s to 2-5 s; deactivation of I Ks was fast enough to complete before the next stimulus during 1-Hz pacing (Wang et al., 1994;Virág et al., 2001;Kurata et al., 2005;Jost et al., 2007).

EAD Formation During β-AS (Validation for LQT1 Model)
In LQT1 patients with smaller I Ks , fatal cardiac events are exercise-induced (tachycardia-related), because adrenergic enhancement of I CaL is no longer counterbalanced by the concomitant stimulation of I Ks ; the smaller increase in I Ks leads to the occurrence of EADs that trigger ventricular tachyarrhythmia. The K05 model, but not the O11 model, could reproduce this I Ks reduction-related EAD formation as a cause of ventricular tachycardia in LQT1 patients during β-AS (Sogo et al., 2016;Kurata et al., 2017). The mTP06b model was also capable of reproducing β-AS-related EAD formation in LQT1 cardiomyocytes, which could clearly be accounted for by the I CaL -and I Ks -dependent bifurcation properties of the model cell ( Figure 5). I CaL /I Ks -dependent properties of EAD formation in the mTP06b model were very similar to those in the K05 model , whereas I CaL /I Ks -dependent HB properties of the mTP06 models were totally different ( Supplementary  Figures S9B, S10). However, EAD formation in the mTP06b model during β-AS was different in a mechanism from that in the K05 model: EADs in the mTP06b model involved spontaneous SR Ca 2+ releases and resulting increments of inward I NCX (Figure 5A and Supplementary Figure S4), while those in the K05 model solely I CaL -reactivation dependent, not involving spontaneous SR Ca 2+ releases . It is likely that the spontaneous SR Ca 2+ release is implicated in EAD formation during β-AS; delayed afterdepolarizations (DADs), known to be induced by spontaneous SR Ca 2+ releases, often accompanied EADs (Priori and Corr, 1990;Volders et al., 2000;Zhao et al., 2012). Amplitudes of spontaneous SR Ca 2+ releases in the mTP06b model appeared to be comparable to or slightly larger than those observed in experimental studies during β-AS (Volders et al., 2000;Zhao et al., 2012) or for LQT2 cardiomyocytes (Choi et al., 2002;Němec et al., 2010Němec et al., , 2016Parikh et al., 2012) unless AP phase 2 was extremely long (Figure 8A and Supplementary Figure S5), while much larger than those reproduced by another HVM model (Trenor et al., 2013) and rabbit ventricular myocyte models (Milberg et al., 2012b;Song et al., 2015;Wilson et al., 2017;Zhong et al., 2018). More sophisticated HVM models incorporating β-ASrelated modulating factors, like those developed by Saucerman et al. (2003) and Kuzumoto et al. (2008), are required for further investigations of the mechanisms of exercise-induced EAD formation in LQT1 HVMs.

Comparisons With Other HVM Models for Bifurcation Phenomena and EAD Mechanisms
EAD Initiation Mechanisms (Roles of I CaL , I NCX , and SR Ca 2+ Release) At least two mechanisms appeared to underlie the initiation of phase-2 EADs in the mTP06b model: (1) I CaL reactivationdependent mechanism which operates and causes EADs even in the absence of spontaneous SR Ca 2+ releases at lower P up and lower pacing rates, and (2) spontaneous SR Ca 2+ release-mediated mechanism activating inward I NCX at higher P up and higher pacing rates (Figures 2A-b, 6B). Coexistence of these two distinct mechanisms for EAD formation have been demonstrated experimentally as well (Zhao et al., 2012).
The major contribution of I CaL to EAD formation was suggested in many previous experimental and theoretical studies for ventricular myocytes (January and Riddle, 1989;Ming et al., 1994;Yamada et al., 2008;Xie et al., 2010;Corrias et al., 2011;Madhvani et al., 2011;Chang et al., 2012a,b;Milberg et al., 2012a;Qu and Chung, 2012;Qu et al., 2013). EAD formation in the mTP06b model with lower P up is also attributable to I CaL reactivation in that reactivated I CaL contributes to V m depolarization (Figures 1C, 2A-a, and  Supplementary Figure S6). In the K05 and O11 models, EADs often emerged in the vicinity of the critical point at which a stable LC appeared during I CaL increases , suggesting that EAD formation depends on I CaL responsible for the instability of EPs and generation of stable LCs. EADs also occurred in the I CaL -enhanced mTP06b model when unstable LCs emerged, but stable LCs were not detected (Supplementary Figure S9A). The slow-fast decomposition analysis of the guinea-pig ventricular myocyte model have suggested that the I CaL -dependent destabilization of a qEP and formation of a stable quasi-LC (qLC) via an HB in the fast subsystem is required for EAD generation in the full system Qu et al., 2013;Song et al., 2015). In the mTP06b model, however, transient trapping of the full system trajectory occurred around the stable and unstable qEPs without forming a stable qLC, indicating that the emergence of a stable qLC is not necessarily needed for EAD formation. This scenario for EAD formation is essentially the same as that in a two-current three-variable AP model (Xie et al., 2014). Nevertheless, the absence of a stable qLC may result in EADs of relatively small amplitudes, as demonstrated for the mTP06b model. As mentioned above, the initiation of EADs in the mTP06b model is attributable to the stabilization of qEPs at depolarized V m in the xs-parameterized fast subsystem, which causes transient trapping of the full system trajectory around the stable qEP; I Kr reduction promotes EAD formation by broadening the region of stable qEPs at depolarized V m (Figure 7).
As another possible mechanism for EAD initiation, many recent experimental and simulation studies have strongly suggested the spontaneous SR Ca 2+ release causing Ca i oscillations, oscillatory increases in inward I NCX , and resulting V m depolarization during β-AS (Choi et al., 2002;Volders et al., 2003;Zhao et al., 2012;Song et al., 2015;Wilson et al., 2017;Zhong et al., 2018) and in I Kr -reduced LQT2 cardiomyocytes (Choi et al., 2002;Kim et al., 2005;Němec et al., 2010Němec et al., , 2016, which is similar to the mechanism for DADs induced by spontaneous SR Ca 2+ releases under Ca 2+ overload conditions or β-AS (e.g., Volders et al., 2003;Zhao et al., 2012) and the Ca 2+ clock mechanism for sinoatrial node cell pacemaking (Maltsev and Lakatta, 2009). The K05 or O11 model could not reproduce the spontaneous SR Ca 2+ release as a cause of phase-2 EADs . In contrast, the mTP06b model could clearly replicate this scenario in a Ca SR -dependent manner (Figures 2Ab, 5A, 6B, 9), while it was not found in the previous study using a modified TP06 model (Vandersickel et al., 2014). Such SR Ca 2+ release-mediated EADs under β-AS conditions (Figure 5) have also been reproduced by the rabbit ventricular myocyte model (Volders et al., 2000;Song et al., 2015;Zhong et al., 2018).
One of the prominent properties of the mTP06 model is the instability of steady-state intracellular Ca 2+ concentrations resulting in the spontaneous SR Ca 2+ release at higher P up to increase Ca SR (Figures 8, 9). The K05 or O11 model did not exhibit spontaneous SR Ca 2+ releases even at higher P up , because steady-state intracellular Ca 2+ concentrations were always stable independently of V m ; although Ca i oscillations occurred during EADs in the K05 and O11 models, these Ca i oscillations were not induced by spontaneous SR Ca 2+ releases but by oscillatory reactivation of I CaL . Thus, this study newly suggests that the occurrence of spontaneous SR Ca 2+ releases and Ca 2+ oscillations as a cause of phase-2 EADs are attributable to instability of intracellular Ca 2+ concentrations in a steady state, destabilization of which leads to spontaneous Ca 2+ oscillations (Figures 8, 9). This scenario, i.e., steadystate destabilization for spontaneous Ca 2+ oscillations involving ryanodine or IP 3 receptors, has previously been suggested by bifurcation analyses for cardiac myocytes (Keizer and Levine, 1996;Tveito et al., 2012) and for other cells (Schuster et al., 2002;Higgins et al., 2006;Kusters et al., 2007). However, the Ca 2+ oscillations reported in these previous studies were much longer in period than those observed in the mTP06b model, not relating to EAD formation. Wilson et al. (2017) demonstrated spontaneous SR Ca 2+ releases and sustained Ca 2+ oscillations in a voltage-clamped rabbit ventricular myocyte model, suggesting instability of intracellular Ca 2+ dynamics; however, dynamical mechanisms for the Ca 2+ oscillation were not clarified by bifurcation analysis. To the best of our knowledge, this is the first report demonstrating instability of steady-state intracellular Ca 2+ concentrations and resulting spontaneous SR Ca 2+ releases that cause EADs in the HVM model (Kurata et al., 2019). In the mTP06b model, enhanced I CaL further contributed to spontaneous SR Ca 2+ releases via the increment of SR Ca 2+ contents and resultant enhancement of the instability of intracellular Ca 2+ dynamics ( Figure 8C). Elevations of Ca ss by spontaneous SR Ca 2+ releases caused transient I CaL reductions due to Ca 2+ -dependent inactivation (Figures 2A-b, 6B, 8A), which may be regarded as a negative feedback mechanism leading to inhibition of EADs.
EAD Termination Mechanisms (Roles of I Ks , I Kr , and I CaL ) The mTP06 model requires I Ks for EAD termination, i.e., repolarization failure occurred abruptly during I Kr inhibition or I CaL enhancement when I Ks was absent or small, whereas I Ks was not necessarily needed in the K05 or O11 model. Thus, EADs in the mTP06b model during pacing appeared to terminate in an I Ks activation-dependent (or stimulus-dependent) manner: The open probability of I Ks channels (xs 2 ) increased progressively in the model cell with relatively small I Kr , i.e., LQT2-like cells (Figure 2A-a). Tran et al. (2009) suggested the major role of the slow I Ks activation for the guinea-pig ventricular myocyte model by the slow-fast decomposition analysis in which the slow I Ks activation gating variable was assumed to be a parameter for the fast subsystem. In the diagram for the slow gating variable-parameterized fast subsystem with a superimposed full system trajectory, gradual increases in the slow variable led the full system trajectory slowly across the stable steady-state branch of qEP and then into the region of the stable periodic branch of qLC through an HB point, resulting in the termination of EADs via a homoclinic bifurcation of qLC. This scenario is known as the Hopf-homoclinic bifurcation mechanism Qu et al., 2013;Song et al., 2015;Huang et al., 2018). Consistent with these previous reports, the mTP06b model exhibited slow I Ks activationdependent EADs in the xs 2 regions of stable and unstable qEPs ( Figure 7A); however, a stable qLC region or homoclinic bifurcation to yield EAD termination was not detected for the fast subsystem of the mTP06b model. EAD terminated simply via the destabilization of a qEP in the I Kr -reduced mTP06b model, suggesting that the Hopf-homoclinic bifurcation scenario is not necessarily applicable.
In the I Ks -eliminated system, EADs would not terminate unless there exist other slow components or factors, such as the slowly inactivating I CaL or late I Na and intracellular Na + accumulation to increase outward I NaK gently. Our previous study using the K05 and O11 models indicated that EAD termination might occur in a slow I CaL inactivation-dependent manner when I Ks was relatively small (Figure 3 in Kurata et al., 2017). However, I CaL inactivation-dependent EAD termination was not clearly detected in the TP06 model. Other candidates for slow variables to cause EAD termination include the slow inactivation of late I Na (Horvath et al., 2013;Trenor et al., 2013;Asakura et al., 2014) not incorporated into the TP06 model and gradual increases in Na i (Chang et al., 2012a;Xie et al., 2015). After cessation of pacing, the I Kr -reduced and/or I CaL -enhanced mTP06b model could exhibit long-term EAD bursts the termination of which was induced by slow elevation of Na i and resulting enhancement of outward I NaK (data not shown). This Na i -dependent mechanism has previously been demonstrated for a rabbit ventricular AP model as well (Chang et al., 2012a). Nevertheless, the slow Na i elevation (intracellular Na + accumulation) is unlikely as a termination mechanism for short-term EADs during pacing at 0.2-2 Hz in the mTP06b model.
As another termination mechanism, stimulus-induced repolarization was observed when APD became very long with stable AP phase 2 (Figures 2A-b, 3B-ii, 5, 6). In the mTP06b model, it was yielded mainly by the prolonged decrease (inactivation) of I CaL due to its slow recovery and resulting outward shift in the total membrane current after the stimulus off. In terms of bifurcation theory, this phenomenon is related to bistability (co-existence of two stable EPs at resting V m and depolarized V m close to AP phase 2) and a transition between the two stable EPs by the stimulus; the transition occurs when following an application of the stimulus current a trajectory of a system starting from one stable EP at the depolarized V m goes outside the attractor basin of the stable EP and enters the attractor basin of the other stable EP at the resting V m (Vinet and Roberge, 1990). Bistability and the stimulus-induced transition between two stable states were found in other cardiomyocyte models (Landau et al., 1990;Vinet and Roberge, 1990). We could not find any experimental evidence for the AP repolarization induced by a stimulus current during stable AP phase 2, but it is theoretically possible. Experimental studies for wider ranges of channel conductance or other parameters may verify that the stimulus-induced repolarization really occurs.

Limitations and Perspectives of Study
As summarized in our preceding article , bifurcation analyses have been used for elucidating the dynamical mechanisms of sinoatrial node pacemaking, abnormal automaticity in ventricular myocytes, generation of biological pacemaker activity, and EAD formation in ventricular myocytes. These theoretical studies have clearly demonstrated the significance of bifurcation analyses for general understanding and systematic description of the dynamical mechanisms of normal and abnormal oscillatory behaviors.
There are many limitations of our study including incompleteness of the model and inconsistency between model predictions and experimental observations, as well as the lack of experimental evidence for bifurcation phenomena in real HVMs. The aim of this study was not to refine but to validate the TP06 model. Nevertheless, more sophisticated HVM models have to be used or developed for more detailed theoretical investigations. As mentioned above, simulated Ca 2+ transients induced by spontaneous SR Ca 2+ releases during AP phase 2 were larger than those observed in many experimental studies. The larger Ca 2+ transients in the model cell may be due to the one compartment SR with weak Ca 2+ leak, which results in higher SR Ca 2+ load and greater Ca 2+ releases during AP phase 2; a two-compartment SR model may be required for reproducing experimentally observed smaller Ca 2+ releases, as suggested previously (Wilson et al., 2017). Moreover, incorporation of more elaborate schemes for the mechanisms of SR Ca 2+ release and intra-SR Ca 2+ transfer (Laver, 2007(Laver, , 2009Chen et al., 2014;Song et al., 2015;Zhong et al., 2018) would also be crucial. Our preceding Tsumoto et al., 2017) and present studies have demonstrated that EAD mechanisms are different depending on models and parameter values. Therefore, we have to test as many models as possible for providing more profound understanding of EAD mechanisms.
In this study, bifurcation analysis was limited to a single cell model. However, EAD-related arrhythmias are suggested to be induced by synchronization of EADs in multiple cells Xie et al., 2010) and also influenced by heterogeneity of ventricular myocytes; because of electrotonic interactions, EAD formation in multicellular or tissue models including epicardial, endocardial and M cell models may be very different in conditions from that in single cell models (Gibb et al., 1994;Huelsing et al., 2000;Weiss et al., 2010;Corrias et al., 2011). Therefore, we need investigations of the mechanisms for EAD formation and for triggering arrhythmias in human ventricles in vivo, which require multicellular (tissue) models, like those used in previous simulation studies de Lange et al., 2012;Vandersickel et al., 2014;Chang et al., 2015;Liu et al., 2018). Despite many limitations, our studies provide significant insights into the dynamical mechanisms of EAD generation in LQT1 and LQT2 HVMs by utilizing recently developed HVM models.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
YaK conception and design of the research and drafted the manuscript. YaK and KT performed the programing, simulations, bifurcation analyses, and analyzed the data. YaK, KT, and KH interpreted the results. YaK, KT, MT, and YuK prepared the figures. YaK, KT, KH, and IH edited and revised the manuscript. YaK, KT, KH, IH, MT, and YuK approved the final version of the manuscript.