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

Abstract

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 (I_Kr) and slow (I_Ks) components of delayed-rectifier K⁺ channel currents, L-type Ca²⁺ channel current (I_CaL), Na⁺/Ca²⁺ exchanger current (I_NCX), and intracellular Ca²⁺ 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 I_CaL 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, I_CaL reactivation–dependent and spontaneous SR Ca²⁺ release–mediated EADs, were detected. (3) Termination of EADs (AP repolarization) during pacing depended on the slow activation of I_Ks. (4) Spontaneous SR Ca²⁺ releases occurred at higher Ca²⁺ uptake rates, attributable to the instability of steady-state intracellular Ca²⁺ 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) (Weiss et al., 2010; 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²⁺ 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; Guo D. et al., 2007; 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²⁺ release from the sarcoplasmic reticulum (SR) (Volders et al., 2000; Choi et al., 2002; Zhao et al., 2012). In our recent theoretical study (Kurata et al., 2017) 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²⁺ release-mediated EADs. With respect to the termination of EADs (AP repolarization), theoretical studies (Tran et al., 2009; 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 (Kurata et al., 2017) 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 (Kurata et al., 2017; 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 parameter-dependent qualitative changes in dynamical behaviors, in non-linear 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 non-paced 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²⁺ 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²⁺ 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 EAD-triggered reentrant arrhythmias in the LQTS human ventricle, as well as for prevention and treatments of life-threatening arrhythmias, like TdPs, in LQTS.

Materials and Methods

Mathematical Modeling for HVMs

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²⁺ 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²⁺ 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, 2004; Guo T. et al., 2007). 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²⁺ uptake rate (P_up) to reduce the Ca²⁺ 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 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²⁺ uptake/release rates for the modified versions, as well as for the original version, are given in Supplementary Table S1.

FIGURE 1

The TP06 model for the normal activity of single HVMs is described as a non-linear dynamical system of 19 first-order ordinary differential equations. The membrane current system includes the Na⁺ channel current (I_Na), I_CaL, I_Kr, I_Ks, 4-aminopyridine-sensitive transient outward current (I_to), inward-rectifier K⁺ channel current (I_K1), background K⁺ (I_pK), Na⁺ (I_bNa) and Ca²⁺ (I_bCa) currents, Na⁺-K⁺ pump current (I_NaK), Na⁺/Ca²⁺ exchanger current (I_NCX), and Ca²⁺ pump current (I_pCa). Time-dependent changes in the membrane potential (V_m) are described by the equation,

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²⁺ (Ca_i), and subspace Ca²⁺ (Ca_ss), while external concentrations of K⁺, Na⁺ and Ca²⁺ 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¹. In addition, the original TP06, mTP06a, and mTP06b models have been implemented in PhysioDesigner as XML-based Physiological Hierarchy Markup Language (PHML)² models. These models can be referred from PHML database (ID938 to 940)³, and simulations of their temporal behaviors can be performed using the software, Flint⁴.

Modeling LQTS Cardiomyocytes With Simulated EADs

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 (Kurata et al., 2017). 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.

Numerical Methods for Dynamic Simulations

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₉₀), 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₉₀ 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₉₀, 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₉₀ 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, 2012, 2013, 2017; Tsumoto et al., 2017), as well as in textbooks (Guckenheimer and Holmes, 1983; Parker and Chua, 1989; Kuznetsov, 2003). In the present study, one- and two-parameter 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 slow-fast 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 (Tran et al., 2009; 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.

Results

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 2A-b); just before applying this second stimulus, the transient depolarization in plateau phase originating from the spontaneous SR Ca²⁺ 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 2

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₉₀ 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 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; see an example of Figure 3C for the point “f” in Figure 3A-i). 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₉₀ of ≤5 s that repolarizes before the next stimulus to an AP with longer APD₉₀ 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₉₀ values of the RF-type AP were almost always more than 10 min.

FIGURE 3

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 3A-i) 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 two-parameter 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²⁺ 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 (Kurata et al., 2017), 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).

FIGURE 4

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²⁺ 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₉₀ 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 two-parameter 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 (Kurata et al., 2017), 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²⁺ 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.

FIGURE 5

Influences of SR Ca²⁺ Cycling, I_NCX and I_CaL on EAD Formation

Following the finding of spontaneous SR Ca²⁺ releases which occurred especially under β-AS with enhanced I_CaL and SR Ca²⁺ uptake, we next examined how SR Ca²⁺ uptake/release (intracellular Ca²⁺ dynamics) and I_NCX regulated by the intracellular Ca²⁺, as well as I_CaL regulated by the subspace Ca²⁺, 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_up–g_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²⁺ 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 (Kurata et al., 2017), while 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²⁺ oscillations at higher P_up, indicating that the SR Ca²⁺ uptake/release machinery destabilizes EPs and thereby induces spontaneous oscillations. When g_Kr was markedly reduced in the mTP06b model, spontaneous SR Ca²⁺ releases to cause transient increases in intracellular Ca²⁺ 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²⁺ release and raised the incidence and frequency of spontaneous Ca²⁺ oscillations to yield EADs or V_m oscillations (Figure 6B-b).

FIGURE 6

Contributions of I_NCX to bifurcations and EAD formation were also explored in relation to those of intracellular Ca²⁺ dynamics, SR Ca²⁺ 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²⁺ exchanger densities or I_NCX (Milberg et al., 2008, 2012b; Pott 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 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 (Kurata et al., 2017): 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 S9B-i) 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_CaL–g_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 small-amplitude 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.

Dynamical Mechanisms for Initiation and Termination of EADs in the mTP06 Model

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 (Tran et al., 2009; Qu et al., 2013; Xie et al., 2014). The I_Ks activation gating variable xs or I_Ks channel open probability (xs²) appears to be a slow variable yielding the termination of EADs (Figure 2A-a, 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², 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² = 0. Stable qEP in the former was destabilized via an HB as xs² increased. The g_Kr reduction led to broadening of the xs² 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² 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² curve. This trapping phenomenon was not observed in the mTP06a model (Figure 7B, right) or the I_Kr-normal mTP06b model (Figure 7A-i, right), because the full system trajectories did not intersect with the stable steady-state branch (qV_E3) before intersecting the steady-state xs² 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.

FIGURE 7

Dynamical Mechanisms of Spontaneous SR Ca²⁺ Release-Mediated EAD

To clarify the dynamical mechanisms of spontaneous SR Ca²⁺ release-mediated EAD formation, we further examined the stability, dynamics and bifurcations of the voltage-clamped mTP06 model. Ca²⁺ 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²⁺ 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²⁺ release and period of spontaneous Ca²⁺ 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. Steady-state intracellular Ca²⁺ concentrations (EPs) in the voltage-clamped 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 Ca²⁺ oscillations occurred in the V_m range of unstable LCs, i.e., between NS₁ and NS₂ (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.

FIGURE 8

To further clarify the Ca_SR-dependent mechanism of spontaneous SR Ca²⁺ releases in the P_up-increased mTP06b model and why SR Ca²⁺ 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 voltage-clamped fast subsystem composed of the voltage-independent state variables f_CaL (Ca²⁺-dependent inactivation gate for I_CaL), R (proportion of closed SR Ca²⁺ 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²⁺ 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²⁺ release, because an increment of Ca_SR (Ca²⁺ refilling of the SR) during Ca²⁺ 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).

FIGURE 9

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 (Kurata et al., 2017) 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.

EAD amplitudes in the mTP06b model (LQT1/2 versions) during pacing (∼30 mV) were smaller than those in the K05 and O11 models (Zimik et al., 2015; Kurata et al., 2017) as well as those in rabbit and guinea-pig ventricular myocyte models (Song et al., 2015; Zhong et al., 2018); however, they were comparable to those in many experimental reports for isolated HVMs (Verkerk et al., 2000; Veldkamp et al., 2001) and human iPS cell-derived LQT2 cardiomyocytes (Itzhaki et al., 2011) as well as for ferret, rabbit, and mouse ventricular myocytes (Marban et al., 1986; Liu et al., 2012; Edwards et al., 2014). Periods of EADs (∼200 ms) were shorter than those in the other HVM models, but comparable to the experimental data from HVMs (Verkerk et al., 2000).

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 (Kurata et al., 2017). Like the other HVM models, the Na_i-variable mTP06b model could partly reproduce the rate-dependent EAD generation in 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 (Kurata et al., 2017), 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 (Kurata et al., 2017), 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²⁺ 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²⁺ releases (Kurata et al., 2017). It is likely that the spontaneous SR Ca²⁺ release is implicated in EAD formation during β-AS; delayed afterdepolarizations (DADs), known to be induced by spontaneous SR Ca²⁺ releases, often accompanied EADs (Priori and Corr, 1990; Volders et al., 2000; Zhao et al., 2012). Amplitudes of spontaneous SR Ca²⁺ 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., 2010, 2016; Parikh 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 β-AS-related 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²⁺ Release)

At least two mechanisms appeared to underlie the initiation of phase-2 EADs in the mTP06b model: (1) I_CaL reactivation-dependent mechanism which operates and causes EADs even in the absence of spontaneous SR Ca²⁺ releases at lower P_up and lower pacing rates, and (2) spontaneous SR Ca²⁺ 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; Guo D. et al., 2007; 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 (Kurata et al., 2017), 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 (Tran et al., 2009; 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²⁺ 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., 2010, 2016), which is similar to the mechanism for DADs induced by spontaneous SR Ca²⁺ releases under Ca²⁺ overload conditions or β-AS (e.g., Volders et al., 2003; Zhao et al., 2012) and the Ca²⁺ clock mechanism for sinoatrial node cell pacemaking (Maltsev and Lakatta, 2009). The K05 or O11 model could not reproduce the spontaneous SR Ca²⁺ release as a cause of phase-2 EADs (Kurata et al., 2017). In contrast, the mTP06b model could clearly replicate this scenario in a Ca_SR-dependent manner (Figures 2A-b, 5A, 6B, 9), while it was not found in the previous study using a modified TP06 model (Vandersickel et al., 2014). Such SR Ca²⁺ 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²⁺ concentrations resulting in the spontaneous SR Ca²⁺ release at higher P_up to increase Ca_SR (Figures 8, 9). The K05 or O11 model did not exhibit spontaneous SR Ca²⁺ releases even at higher P_up, because steady-state intracellular Ca²⁺ 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²⁺ releases but by oscillatory reactivation of I_CaL (Kurata et al., 2017). Thus, this study newly suggests that the occurrence of spontaneous SR Ca²⁺ releases and Ca²⁺ oscillations as a cause of phase-2 EADs are attributable to instability of intracellular Ca²⁺ concentrations in a steady state, destabilization of which leads to spontaneous Ca²⁺ oscillations (Figures 8, 9). This scenario, i.e., steady-state destabilization for spontaneous Ca²⁺ oscillations involving ryanodine or IP₃ 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²⁺ 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²⁺ releases and sustained Ca²⁺ oscillations in a voltage-clamped rabbit ventricular myocyte model, suggesting instability of intracellular Ca²⁺ dynamics; however, dynamical mechanisms for the Ca²⁺ 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²⁺ concentrations and resulting spontaneous SR Ca²⁺ 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²⁺ releases via the increment of SR Ca²⁺ contents and resultant enhancement of the instability of intracellular Ca²⁺ dynamics (Figure 8C). Elevations of Ca_ss by spontaneous SR Ca²⁺ releases caused transient I_CaL reductions due to Ca²⁺-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²) 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 (Tran et al., 2009; Qu et al., 2013; Song et al., 2015; Huang et al., 2018). Consistent with these previous reports, the mTP06b model exhibited slow I_Ks activation-dependent EADs in the xs² 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 (Kurata et al., 2017), 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²⁺ transients induced by spontaneous SR Ca²⁺ releases during AP phase 2 were larger than those observed in many experimental studies. The larger Ca²⁺ transients in the model cell may be due to the one compartment SR with weak Ca²⁺ leak, which results in higher SR Ca²⁺ load and greater Ca²⁺ releases during AP phase 2; a two-compartment SR model may be required for reproducing experimentally observed smaller Ca²⁺ releases, as suggested previously (Wilson et al., 2017). Moreover, incorporation of more elaborate schemes for the mechanisms of SR Ca²⁺ release and intra-SR Ca²⁺ transfer (Laver, 2007, 2009; Chen et al., 2014; Song et al., 2015; Zhong et al., 2018) would also be crucial. Our preceding (Kurata et al., 2017; 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 (Sato et al., 2009; 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 (Weiss et al., 2010; 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.

References

• 1

  AdeniranI.El HarchiA.HancoxJ. C.ZhangH. (2012). Proarrhythmia in KCNJ2-linked short QT syndrome: insights from modelling.Cardiovasc. Res.9466–76. 10.1093/cvr/cvs082

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  AndersonC. L.DelisleB. P.AnsonB. D.KilbyJ. A.WillM. L.TesterD. J.et al (2006). Most LQT2 mutations reduce Kv11.1 (hERG) current by a class 2 (trafficking-deficient) mechanism.Circulation113365–373. 10.1161/CIRCULATIONAHA.105.570200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  AntzelevitchC.ShimizuW.YanG. X.SicouriS.WeissenburgerJ.NesterenkoV. V.et al (1999). The M cell: its contribution to the ECG and to normal and abnormal electrical function of the heart.J. Cardiovasc. Electrophysiol.101124–1152. 10.1111/j.1540-8167.1999.tb00287.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  AsakuraK.ChaC. Y.YamaokaH.HorikawaY.MemidaH.PowellT.et al (2014). EAD and DAD mechanisms analyzed by developing a new human ventricular cell model.Prog. Biophys. Mol. Biol.11611–24. 10.1016/j.pbiomolbio.2014.08.008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  BristonS. J.DibbK. M.SolaroR. J.EisnerD. A.TraffordA. W. (2014). Balanced changes in Ca buffering by SERCA and troponin contribute to Ca handling during β-adrenergic stimulation in cardiac myocytes.Cardiovasc. Res.104347–354. 10.1093/cvr/cvu201

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  CarmelietE. (1992). Voltage- and time-dependent block of the delayed K+ current in cardiac myocytes by dofetilide.J. Pharmacol. Exp. Ther.262809–817.

  • Pubmed Abstract
  • Google Scholar

• 7

  ChangM. G.ChangC. Y.de LangeE.XuL.O’RourkeB.KaragueuzianH. S.et al (2012a). Dynamics of early afterdepolarization-mediated triggered activity in cardiac monolayers.Biophys. J.1022706–2714. 10.1016/j.bpj.2012.05.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  ChangM. G.SatoD.de LangeE.LeeJ. H.KaragueuzianH. S.GarfinkelA.et al (2012b). Bi-stable wave propagation and early afterdepolarization-mediated cardiac arrhythmias.Heart Rhythm9115–122. 10.1016/j.hrthm.2011.08.014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  ChangP. C.WoH. T.LeeH. L.LinS. F.WenM. S.ChuY.et al (2015). Role of sarcoplasmic reticulum calcium in development of secondary calcium rise and early afterdepolarizations in long QT syndrome rabbit model.PLoS One10:e0123868. 10.1371/journal.pone.0123868

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  ChenW.WangR.ChenB.ZhongX.KongH.BaiY.et al (2014). The ryanodine receptor store-sensing gate controls Ca2+ waves and Ca2+-triggered arrhythmias.Nat. Med.20184–192. 10.1038/nm.3440

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  ChiangC. E.RodenD. M. (2000). The long QT syndromes: genetic basis and clinical implications.J. Am. Coll. Cardiol.361–12. 10.1016/s0735-1097(00)00716-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  ChoiB. R.BurtonF.SalamaG. (2002). Cytosolic Ca2+ triggers early afterdepolarizations and torsade de pointes in rabbit hearts with type 2 long QT syndrome.J. Physiol.543615–631. 10.1113/jphysiol.2002.024570

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  ChouabeC.NeyroudN.GuicheneyP.LazdunskiM.RomeyG.BarhaninJ. (1997). Properties of KvLQT1 K+ channel mutations in romano-ward and jervell and lange-nielsen inherited cardiac arrhythmias.EMBO J.165472–5479. 10.1093/emboj/16.17.5472

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  CorriasA.GilesW.RodriguezB. (2011). Ionic mechanisms of electrophysiological properties and repolarization abnormalities in rabbit Purkinje fibers.Am. J. Physiol. Heart Circ. Physiol.300H1806–H1813. 10.1113/jphysiol.2002.024570

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  de LangeE.XieY.QuZ. (2012). Synchronization of early afterdepolarizations and arrhythmogenesis in heterogeneous cardiac tissue models.Biophys. J.103365–373. 10.1016/j.bpj.2012.06.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  EdwardsA. G.GrandiE.HakeJ. E.PatelS.LiP.MiyamotoS.et al (2014). Nonequilibrium reactivation of Na+ current drives early afterdepolarizations in mouse ventricle.Circ. Arrhythm. Electrophysiol.71205–1213. 10.1161/CIRCEP.113.001666

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  GibbW. J.WagnerM. B.LeshM. D. (1994). Effects of simulated potassium blockade on the dynamics of triggered cardiac activity.J. Theor. Biol.168245–257. 10.1006/jtbi.1994.1106

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  GuckenheimerJ.HolmesP. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields.New York, NY: Springer-Verlag.

  • Google Scholar

• 19

  GuoD.ZhaoX.WuY.LiuT.KoweyP. R.YanG. X. (2007). L-type calcium current reactivation contributes to arrhythmogenesis associated with action potential triangulation.J. Cardiovasc. Electrophysiol.18196–203. 10.1111/j.1540-8167.2006.00698.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  GuoT.AiX.ShannonT. R.PogwizdS. M.BersD. M. (2007). Intra-sarcoplasmic reticulum free [Ca2+] and buffering in arrhythmogenic failing rabbit heart.Circ. Res.101802–810. 10.1161/CIRCRESAHA.107.152140

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  HigginsE. R.CannellM. B.SneydJ. (2006). A buffering SERCA pump in models of calcium dynamics.Biophys. J.91151–163. 10.1529/biophysj.105.075747

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  HimenoY.SaraiN.MatsuokaS.NomaA. (2008). Ionic mechanisms underlying the positive chronotropy induced by beta1-adrenergic stimulation in guinea pig sinoatrial node cells: a simulation study.J. Physiol. Sci.5853–65. 10.2170/physiolsci.RP015207

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  HorvathB.BanyaszT.JianZ.HegyiB.KistamasK.NanasiP. P.et al (2013). Dynamics of the late Na+ current during cardiac action potential and its contribution to afterdepolarizations.J. Mol. Cell. Cardiol.6459–68. 10.1016/j.yjmcc.2013.08.010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  HuangX.SongZ.QuZ. (2018). Determinants of early afterdepolarization properties in ventricular myocyte models.PLoS Comput. Biol.14:e1006382. 10.1371/journal.pcbi.1006382

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  HuelsingD. J.SpitzerK. W.PollardA. E. (2000). Electrotonic suppression of early afterdepolarizations in isolated rabbit purkinje myocytes.Am. J. Physiol. Heart Circ. Physiol.279H250–H259. 10.1152/ajpheart.2000.279.1.H250

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  ItzhakiI.MaizelsL.HuberI.Zwi-DantsisL.CaspiO.WintersternA.et al (2011). Modelling the long QT syndrome with induced pluripotent stem cells.Nature471225–229. 10.1038/nature09747

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  JanuaryC. T.RiddleJ. M. (1989). Early afterdepolarizations: mechanism of induction and block; a role for L-type Ca2+ current.Circ. Res.64977–990. 10.1161/01.res.64.5.977

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  JanuaryC. T.RiddleJ. M.SalataJ. J. (1988). A model for early afterdepolarizations: induction with the Ca2+ channel agonist Bay K 8644.Circ. Res.62563–571. 10.1161/01.res.62.3.563

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  JostN.PappJ. G.VarróA. (2007). Slow delayed rectifier potassium current (IKs) and the repolarization reserve.Ann. Noninvasive Electrocardiol.1264–78. 10.1111/j.1542-474X.2007.00140.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  JostN.VirágL.BitayM.TakácsJ.LengyelC.BiliczkiP.et al (2005). Restricting excessive cardiac action potential and QT prolongation: a vital role for IKs in human ventricular muscle.Circulation1121392–1399. 10.1161/CIRCULATIONAHA.105.550111

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  KazbanovI. V.ten TusscherK. H.PanfilovA. V. (2016). Effects of heterogeneous diffuse fibrosis on arrhythmia dynamics and mechanism.Sci. Rep.6:20835. 10.1038/srep20835

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  KeizerJ.LevineL. (1996). Ryanodine receptor adaptation and Ca2+-induced Ca2+ release-dependent Ca2+ oscillations.Biophys. J.713477–3487. 10.1016/S0006-3495(96)79543-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  KimJ. J.NěmecJ.LiQ.SalamaG. (2005). Synchronous systolic subcellular Ca2+-elevations underlie ventricular arrhythmia in drug-induced long QT type 2.Circ. Arrhythm. Electrophysiol.8703–712. 10.1161/CIRCEP.114.002214

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  KondoT.HisatomeI.YoshimuraS.MahatiE.NotsuT.LiP.et al (2016). Characterization of the novel mutant A78T-HERG from a long QT syndrome type 2 patient: instability of the mutant protein and stabilization by heat shock factor 1.J. Arrhythm.32433–440. 10.1016/j.joa.2015.10.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  Krogh-MadsenT.SchafferP.SkriverA. D.TaylorL. K.PelzmannB.KoidlB.et al (2005). An ionic model for rhythmic activity in small clusters of embryonic chick ventricular cells.Am. J. Physiol. Heart Circ. Physiol.289H398–H413. 10.1152/ajpheart.00683.2004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  KurataY.HisatomeI.MatsudaH.ShibamotoT. (2005). Dynamical mechanisms of pacemaker generation in IK1-downregulated human ventricular myocytes: insights from bifurcation analyses of a mathematical model.Biophys. J.892865–2887. 10.1529/biophysj.105.060830

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  KurataY.HisatomeI.ShibamotoT. (2012). Roles of sarcoplasmic reticulum Ca2+ cycling and Na+/Ca2+ exchanger in sinoatrial node pacemaking: insights from bifurcation analysis of mathematical models.Am. J. Physiol. Heart Circ. Physiol.302H2285–H2300. 10.1152/ajpheart.00221.2011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  KurataY.HisatomeI.TanidaM.ShibamotoT. (2013). Effect of hyperpolarization-activated current If on robustness of sinoatrial node pacemaking: theoretical study on influence of intracellular Na+ concentration.Am. J. Physiol. Heart Circ. Physiol.304H1337–H1351. 10.1152/ajpheart.00777.2012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  KurataY.MatsudaH.HisatomeI.ShibamotoT. (2008). Regional difference in dynamical property of sinoatrial node pacemaking: role of Na+ channel current.Biophys. J.95951–977. 10.1529/biophysj.107.112854

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  KurataY.TsumotoK.HayashiK.HisatomeI.KudaY.TanidaM. (2019). Multiple dynamical mechanisms of phase-2 early afterdepolarizations in a human ventricular myocyte model: involvement of spontaneous SR Ca2+ release.BioRxiv[Preprint]. Available at: https://www.biorxiv.org/content/10.1101/613182v1(accessed June 12, 2019).

  • Google Scholar

• 41

  KurataY.TsumotoK.HayashiK.HisatomeI.TanidaM.KudaY.et al (2017). 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.312H106–H127. 10.1152/ajpheart.00115.2016

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  KustersJ. M.CortesJ. M.van MeerwijkW. P.YpeyD. L.TheuvenetA. P.GielenC. C. (2007). Hysteresis and bistability in a realistic cell model for calcium oscillations and action potential firing.Phys. Rev. Lett.98:098107. 10.1103/PhysRevLett.98.098107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  KuznetsovY. A. (2003). Elements of Applied Bifurcation Theory, 3rd Edn, New York, NY: Springer-Verlag.

  • Google Scholar

• 44

  KuzumotoM.TakeuchiA.NakaiH.OkaC.NomaA.MatsuokaS. (2008). Simulation analysis of intracellular Na+ and Cl- homeostasis during β1-adrenergic stimulation of cardiac myocyte.Prog. Biophys. Mol. Biol.96171–186. 10.1016/j.pbiomolbio.2007.07.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  LandauM.LorenteP.MichaelsD.JalifeJ. (1990). Bistabilities and annihilation phenomena in electrophysiological cardiac models.Circ. Res.661658–1672. 10.1161/01.res.66.6.1658

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  LaverD. R. (2007). Ca2+ stores regulate ryanodine receptor Ca2+ release channels via luminal and cytosolic Ca2+ sites.Biophys. J.923541–3555. 10.1529/biophysj.106.099028

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  LaverD. R. (2009). Luminal Ca2+ activation of cardiac ryanodine receptors by luminal and cytoplasmic domains.Eur. Biophys. J.3919–26. 10.1007/s00249-009-0417-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  LiuG. X.ChoiB. R.ZivO.LiW.de LangeE.QuZ.et al (2012). Differential conditions for early after-depolarizations and triggered activity in cardiomyocytes derived from transgenic LQT1 and LQT2 rabbits.J. Physiol.5901171–1180. 10.1113/jphysiol.2011.218164

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  LiuW.KimT. Y.HuangX.LiuM. B.KorenG.ChoiB. R.et al (2018). Mechanisms linking T-wave alternans to spontaneous initiation of ventricular arrhythmias in rabbit models of long QT syndrome.J. Physiol.5961341–1355. 10.1113/JP275492

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  LuoC.RudyY. (1991). A model of the ventricular cardiac action potential: depolarization, repolarization, and their interaction.Circ. Res.681501–1526. 10.1161/01.res.68.6.1501

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  MadhvaniR. V.XieY.PantazisA.GarfinkelA.QuZ.WeissJ. N.et al (2011). Shaping a new Ca²⁺ conductance to suppress early afterdepolarizations in cardiac myocytes.J. Physiol.5896081–6092. 10.1113/jphysiol.2011.219600

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  MaltsevV. A.LakattaE. G. (2009). Synergism of coupled subsarcolemmal Ca²⁺ clocks and sarcolemmal voltage clocks confers robust and flexible pacemaker function in a novel pacemaker cell model.Am. J. Physiol. Heart Circ. Physiol.296H594–H615. 10.1152/ajpheart.01118.2008

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  MaltsevV. A.LakattaE. G. (2010). A novel quantitative explanation for the autonomic modulation of cardiac pacemaker cell automaticity via a dynamic system of sarcolemmal and intracellular proteins.Am. J. Physiol. Heart Circ. Physiol.298H2010–H2023. 10.1152/ajpheart.00783.2009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  MarbanE.RobinsonS. W.WierW. G. (1986). Mechanisms of arrhythmogenic delayed and early afterdepolarizations in ferret ventricular muscle.J. Clin. Invest.781185–1192. 10.1172/JCI112701

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  MilbergP.FinkM.PottC.FrommeyerG.BiertzJ.OsadaN.et al (2012a). Blockade of ICa suppresses early afterdepolarizations and reduces transmural dispersion of repolarization in a whole heart model of chronic heart failure.Br. J. Pharmacol.166557–568. 10.1111/j.1476-5381.2011.01721.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  MilbergP.PottC.FrommeyerG.FinkM.RuheM.MatsudaT.et al (2012b). Acute inhibition of the Na+/Ca2+ exchanger reduces proarrhythmia in an experimental model of chronic heart failure.Heart Rhythm9570–578. 10.1016/j.hrthm.2011.11.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  MilbergP.PottC.FinkM.FrommeyerG.MatsudaT.BabaA.et al (2008). Inhibition of the Na+/Ca2+ exchanger suppresses torsades de pointes in an intact heart model of long QT syndrome-2 and long QT syndrome-3.Heart Rhythm51444–1452. 10.1016/j.hrthm.2008.06.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  MingZ.NordinC.AronsonR. S. (1994). Role of L-type calcium channel window current in generating current-induced early afterdepolarization.J. Cardiovasc. Electrophysiol.5323–334. 10.1111/j.1540-8167.1994.tb01169.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  NěmecJ.KimJ. J.GabrisB.SalamaG. (2010). Calcium oscillations and T-wave lability precede ventricular arrhythmias in acquired long QT type 2.Heart Rhythm71686–1694. 10.1016/j.hrthm.2010.06.032

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  NěmecJ.KimJ. J.SalamaG. (2016). The link between abnormal calcium handling and electrical instability in acquired long QT syndrome–Does calcium precipitate arrhythmic storms?Prog. Biophys. Mol. Biol.120210–221. 10.1016/j.pbiomolbio.2015.11.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  O’HaraT.RudyY. (2012). Quantitative comparison of cardiac myocyte electrophysiology and response to drugs in human and nonhuman species.Am. J. Physiol. Heart Circ. Physiol.302H1023–H1030. 10.1152/ajpheart.00785.2011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  O’HaraT.VirágL.VarróA.RudyY. (2011). Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation.PLoS Comput. Biol.7:e1002061. 10.1371/journal.pcbi.1002061

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  ParikhA.MantravadiR.KozhevnikovD.RocheM. A.YeY.OwenL. J.et al (2012). Ranolazine stabilizes cardiac ryanodine receptors: a novel mechanism for the suppression of early afterdepolarization and torsades de pointes in long QT type 2.Heart Rhythm9953–960. 10.1016/j.hrthm.2012.01.010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  ParkerT. S.ChuaL. O. (1989). Practical Numerical Algorithms for Chaotic Systems.New York, NY: Springer-Verlag.

  • Google Scholar

• 65

  PottC.MuszynskiA.RuheM.BögeholzN.SchulteJ. S.MilbergP.et al (2012). Proarrhythmia in a non-failing murine model of cardiac-specific Na+/Ca2+ exchanger overexpression: whole heart and cellular mechanisms.Basic Res. Cardiol.1071–13. 10.1007/s00395-012-0247-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  PrioriS. G.CorrP. B. (1990). Mechanisms underlying early and delayed afterdepolarizations induced by catecholamines.Am. J. Physiol.258H1796–H1805. 10.1152/ajpheart.1990.258.6.H1796

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  QuZ.ChungD. (2012). Mechanisms and determinants of ultralong action potential duration and slow rate-dependence in cardiac myocytes.PLoS One7:e43587. 10.1371/journal.pone.0043587

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  QuZ.XieL. H.OlceseR.KaragueuzianH. S.ChenP. S.GarfinkelA.et al (2013). Early afterdepolarizations in cardiac myocytes: beyond reduced repolarization reserve.Cardiovasc. Res.996–15. 10.1093/cvr/cvt104

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  RodenD. M.LazzaraR.RosenM.SchwartzP. J.TowbinJ.VincentG. M. (1996). Multiple mechanisms in the long-QT syndrome. Current knowledge, gaps, and future directions. the SADS foundation task force on LQTS.Circulation941996–2012. 10.1161/01.cir.94.8.1996

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 70

  SatoD.XieL. H.SovariA. A.TranD. X.MoritaN.XieF.et al (2009). Synchronization of chaotic early afterdepolarizations in the genesis of cardiac arrhythmias.Proc. Natl. Acad. Sci. U.S.A.1062983–2988. 10.1073/pnas.0809148106

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  SaucermanJ. J.BruntonL. L.MichailovaA. P.McCullochA. D. (2003). Modeling β-adrenergic control of cardiac myocytes contractility in silico.J. Biol. Chem.27847997–48003. 10.1074/jbc.M308362200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  SchusterS.MarhlM.HöferT. (2002). Modelling of simple and complex calcium oscillations. From single-cell responses to intercellular signalling.Eur. J. Biochem.2691333–1355. 10.1046/j.0014-2956.2001.02720.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  ShannonT. R.GuoT.BersD. M. (2003). Ca2+ scraps: local depletions of free [Ca2+] in cardiac sarcoplasmic reticulum during contractions leave substantial Ca2+ reserve.Circ. Res.9340–45. 10.1161/01.RES.0000079967.11815.19

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 74

  ShannonT. R.WangF.PuglisiJ.WeberC.BersD. M. (2004). A mathematical treatment of integrated Ca dynamics within the ventricular myocyte.Biophys. J.873351–3371. 10.1529/biophysj.104.047449

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  ShimizuW. (2013). Update of diagnosis and management inherited cardiac arrhythmias.Circ. J.772867–2872. 10.1253/circj.cj-13-1217

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  ShimizuW.HorieM. (2011). Phenotypic manifestations of mutations in genes encoding subunits of cardiac potassium channels.Circ. Res.10997–109. 10.1161/CIRCRESAHA.110.224600

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  SogoT.MorikawaK.KurataY.LiP.IchinoseT.YuasaS.et al (2016). Electrophysiological properties of iPS cell-derived cardiomyocytes from a patient with long QT syndrome type 1 harboring the novel mutation M437V of KCNQ1.Regener. Ther.49–17. 10.1016/j.reth.2015.12.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 78

  SongZ.KoC. Y.NivalaM.WeissJ. N.QuZ. (2015). Calcium-voltage coupling in the genesis of early and delayed afterdepolarizations in cardiac myocytes.Biophys. J.1081908–1921. 10.1016/j.bpj.2015.03.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  ten TusscherK. H.HrenR.PanfilovA. V. (2007). Organization of ventricular fibrillation in the human heart.Circ. Res.100e87–e101. 10.1161/CIRCRESAHA.107.150730

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 80

  ten TusscherK. H.PanfilovA. V. (2006). Alternans and spiral breakup in a human ventricular tissue model.Am. J. Physiol. Heart Circ. Physiol.291H1088–H1100. 10.1152/ajpheart.00109.2006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 81

  TranD. X.SatoD.YochelisA.WeissJ. N.GarfinkelA.QuZ. (2009). Bifurcation and chaos in a model of cardiac early afterdepolarizations.Phys. Rev. Lett.102:258103. 10.1103/PhysRevLett.102.258103

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  TrenorB.CardonaK.SaizJ.RajamaniS.BelardinelliL.GilesW. R. (2013). Carbon monoxide effects on human ventricular action potential assessed by mathematical simulation.Front. Physiol.4:282. 10.3389/fphys.2013.00282

  • CrossRef
  • Google Scholar

• 83

  TsumotoK.KurataY.FurutaniK.KurachiY. (2017). Hysteretic dynamics of multi-stable early afterdepolarisations with repolarisation reserve attenuation: a potential dynamical mechanism for cardiac arrhythmias.Sci. Rep.7:10771. 10.1038/s41598-017-11355-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 84

  TveitoA.LinesG. T.HakeJ.EdwardsA. G. (2012). Instabilities of the resting state in a mathematical model of calcium handling in cardiac myocytes.Math. Biosci.23697–107. 10.1016/j.mbs.2012.02.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 85

  VandersickelN.KazbanovI. V.NuitermansA.WeiseL. D.PanditR.PanfilovA. V. (2014). A study of early afterdepolarizations in a model for human ventricular tissue.PLoS One9:e84595. 10.1371/journal.pone.0084595

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 86

  VeldkampM. W.VerkerkA. O.van GinnekenA. C.BaartscheerA.SchumacherC.de JongeN.et al (2001). Norepinephrine induces action potential prolongation and early afterdepolarizations in ventricular myocytes isolated from human end-stage failing hearts.Eur. Heart J.22955–963. 10.1053/euhj.2000.2499

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 87

  VerkerkA. O.VeldkampM. W.de JongeN.WildersR.van GinnekenA. C. (2000). Injury current modulates afterdepolarizations in single human ventricular cells.Cardiovasc. Res.47124–132. 10.1016/s0008-6363(00)00064-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 88

  VinetA.RobergeF. A. (1990). A model study of stability and oscillations in the myocardial cell membrane.J. Theor. Biol.147377–412. 10.1016/s0022-5193(05)80495-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 89

  VirágL.IostN.OpincariuM.SzolnokyJ.SzécsiJ.BogátsG.et al (2001). The slow component of the delayed rectifier potassium current in undiseased human ventricular myocytes.Cardiovasc. Res.49790–797. 10.1016/s0008-6363(00)00306-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 90

  ViswanathanP. C.RudyY. (1999). Pause induced early afterdepolarizations in the long QT syndrome: a simulation study.Cardiovasc. Res.42530–542. 10.1016/s0008-6363(99)00035-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 91

  VoldersP. G.StenglM.van OpstalJ. M.GerlachU.SpatjensR. L.BeekmanJ. D.et al (2003). Probing the contribution of IKs to canine ventricular repolarization: key role for beta-adrenergic receptor stimulation.Circulation1072753–2760. 10.1161/01.CIR.0000068344.54010.B3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 92

  VoldersP. G.VosM. A.SzaboB.SipidoK. R.de GrootS. H.GorgelsA. P.et al (2000). Progress in the understanding of cardiac early afterdepolarizations and torsades de pointes: time to revise current concepts.Cardiovasc. Res.46376–392. 10.1016/s0008-6363(00)00022-5

  • CrossRef
  • Google Scholar

• 93

  WangZ.FerminiB.NattelS. (1994). Rapid and slow components of delayed rectifier current in human atrial myocytes.Cardiovasc. Res.281540–1546. 10.1093/cvr/28.10.1540

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 94

  WeissJ. N.GarfinkelA.KaragueuzianH. S.ChenP. S.QuZ. (2010). Early afterdepolarizations and cardiac arrhythmias.Heart Rhythm71891–1899. 10.1016/j.hrthm.2010.09.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 95

  WienerR.HaitinY.ShamgarL.Fernández-AlonsoM. C.MartosA.Chomsky-HechtO.et al (2008). The KCNQ1 (Kv7.1) COOH terminus, a multitiered scaffold for subunit assembly and protein interaction.J. Biol. Chem.2835815–5830. 10.1074/jbc.M707541200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 96

  WilsonD.ErmentroutB.NěmecJ.SalamaG. (2017). A model of cardiac ryanodine receptor gating predicts experimental Ca2+-dynamics and Ca2+-triggered arrhythmia in the long QT syndrome.Chaos27:093940. 10.1063/1.5000711

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 97

  XieY.IzuL. T.BersD. M.SatoD. (2014). Arrhythmogenic transient dynamics in cardiac myocytes.Biophys. J.1061391–1397. 10.1016/j.bpj.2013.12.050

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 98

  XieY.SatoD.GarfinkelA.QuZ.WeissJ. N. (2010). So little source, so much sink: requirements for afterdepolarizations to propagate in tissue.Biophys. J.991408–1415. 10.1016/j.bpj.2010.06.042

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 99

  XieY.LiaoZ.GrandiE.ShiferawY.BersD. M. (2015). Slow [Na]_i changes and positive feedback between membrane potential and [Ca]_i underlie intermittent early afterdepolarizations and arrhythmias.Circ. Arrhythm. Electrophysiol.81472–1480. 10.1161/CIRCEP.115.003085

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 100

  YamadaM.OhtaK.NiwaA.TsujinoN.NakadaT.HiroseM. (2008). Contribution of L-Type Ca2+ channels to early afterdepolarizations induced by IKr and IKs channel suppression in guinea pig ventricular myocytes.J. Membr. Biol.222151–166. 10.1007/s00232-008-9113-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 101

  ZengJ.RudyY. (1995). Early afterdepolarizations in cardiac myocytes: mechanism and rate dependence.Biophys. J.68949–964. 10.1016/S0006-3495(95)80271-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 102

  ZhaoZ.WenH.FefelovaN.AllenC.BabaA.MatsudaT.et al (2012). Revisiting the ionic mechanisms of early afterdepolarizations in cardiomyocytes: predominant by Ca waves or Ca currents?Am. J. Physiol. Heart Circ. Physiol.302H1636–H1644. 10.1152/ajpheart.00742.2011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 103

  ZhongM.ReesC. M.TerentyevD.ChoiB. R.KorenG.KarmaA. (2018). NCX-mediated subcellular Ca2+ dynamics underlying early afterdepolarizations in LQT2 cardiomyocytes.Biophys. J.1151019–1032. 10.1016/j.bpj.2018.08.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 104

  ZimikS.VandersickelN.NayakA. R.PanfilovA. V.PanditR. (2015). A comparative study of early afterdepolarization-mediated fibrillation in two mathematical models for human ventricular cells.PLoS One10:e0130632. 10.1371/journal.pone.0130632

  • Pubmed Abstract
  • CrossRef
  • Google Scholar