Determinants of early afterdepolarization properties in ventricular myocyte models

Early afterdepolarizations (EADs) are spontaneous depolarizations during the repolarization phase of an action potential in cardiac myocytes. It is widely known that EADs are promoted by increasing inward currents and/or decreasing outward currents, a condition called reduced repolarization reserve. Recent studies based on bifurcation theories show that EADs are caused by a dual Hopf-homoclinic bifurcation, bringing in further mechanistic insights into the genesis and dynamics of EADs. In this study, we investigated the EAD properties, such as the EAD amplitude, the inter-EAD interval, and the latency of the first EAD, and their major determinants. We first made predictions based on the bifurcation theory and then validated them in physiologically more detailed action potential models. These properties were investigated by varying one parameter at a time or using parameter sets randomly drawn from assigned intervals. The theoretical and simulation results were compared with experimental data from the literature. Our major findings are that the EAD amplitude and takeoff potential exhibit a negative linear correlation; the inter-EAD interval is insensitive to the maximum ionic current conductance but mainly determined by the kinetics of ICa,L and the dual Hopf-homoclinic bifurcation; and both inter-EAD interval and latency vary largely from model to model. Most of the model results generally agree with experimental observations in isolated ventricular myocytes. However, a major discrepancy between modeling results and experimental observations is that the inter-EAD intervals observed in experiments are mainly between 200 and 500 ms, irrespective of species, while those of the mathematical models exhibit a much wider range with some models exhibiting inter-EAD intervals less than 100 ms. Our simulations show that the cause of this discrepancy is likely due to the difference in ICa,L recovery properties in different mathematical models, which needs to be addressed in future action potential model development.


Introduction
Under diseased conditions or influence of drugs, cardiac myocytes can exhibit early afterdepolarizations (EADs) [1][2][3]. EADs are depolarization events during the repolarizing phase of an action potential (AP), which are known to be arrhythmogenic [4][5][6][7]. Many experimental and computational studies have been carried out, which have greatly improved our understanding of the causes and mechanisms of the genesis of EADs. It is well known that EADs can occur in an AP when inward currents are increased and/or outward currents are reduced, a condition called reduced repolarization reserve [8]. Under this condition, L-type calcium (Ca 2+ ) current (I Ca,L ) can be reactivated to cause depolarizations in the repolarization phase of the AP to manifest as EADs. The importance of I Ca,L reactivation for EAD genesis has been widely demonstrated in experiments [1,9] and computer simulations [10,11]. Recent studies [12][13][14][15] using bifurcation theories have brought in additional mechanistic insights into the genesis of EADs, which show that EADs are oscillations originating via a supercritical or subcritical Hopf bifurcation and terminating via a homoclinic bifurcation, or via an unstable manifold of a saddle focus fixed point in the full AP dynamics [14]. In other words, irrespective to the specific ionic causes, EADs are oscillations caused by a Hopf bifurcation [16,17] after which the quasi-equilibrium state becomes unstable and the system oscillates around the unstable equilibrium. As the slow outward currents grow gradually, the oscillation amplitude increases until a new bifurcation point-the homoclinic bifurcation at which the oscillation stops. A detailed discussion on the links between the ionic causes and nonlinear dynamics for the genesis of EADs was presented in our previous review article [18].
Despite the wide experimental and computational studies on the ionic causes and dynamical mechanisms of EAD genesis, less attention has been paid on EAD properties, such as the EAD latency (the time from the upstroke of the AP to the upstroke of the first EAD), the inter-EAD interval, and the EAD amplitude. Although it is well-known that reactivation of I Ca,L is required for EAD genesis, since EADs are a collective behavior arising from the interactions of many ionic currents, it is unclear how these ionic currents affect the EAD properties and what are the major determinants. For example, since I Ca,L plays a critical role in the genesis of EADs, one would intuitively expect that increasing the maximum conduction of I Ca,L might increase the amplitude of EADs, but as we show in this study that this is not the case. On the other hand, understanding the EAD properties and their determinants is important for understanding the mechanisms of EAD-related arrhythmogenesis. For example, in an early experimental study [19], Damiano and Rosen showed that phase-2 EADs cannot propagate as premature ventricular complexes (PVCs) while phase-3 EADs can propagate as PVCs. This was also shown in our simulation studies [20][21][22]. Therefore, understanding what determine the EAD amplitude and takeoff potential may provide insights into EAD propagation to produce PVCs. If a PVC is a direct consequence of EAD propagation, then the EAD latency may provide information for the coupling interval between a sinus beat and the following PVC. EADs are also thought to be responsible for focal arrhythmias in the heart, and if this is true, then the oscillation frequency of EADs should be the same as the excitation frequency of ventricular arrhythmias. Furthermore, understanding the EAD properties and their determinants can also be important for the development of robust mathematical AP models. For example, we observed a discrepancy in inter-EAD interval between those from some widely used AP models and the experimental data. Experimental measurements in ventricular myocytes isolated from animal and human hearts almost exclusively show that the inter-EAD intervals are greater than 200 ms with few exceptions (Table 1). However, many ventricular myocyte AP models show inter-EAD intervals much shorter than 200 ms [11,20,[23][24][25][26][27], raising a question on what ionic current properties have been missed in these models.
The previous computational studies mainly focused on the ionic causes and dynamical mechanisms of the genesis of EADs, and to our knowledge, no studies have been carried out to investigate the EAD properties and their determinants. Furthermore, a close comparison of the EAD properties from mathematical models with those from experimental recordings has not been done. The objective of this study is to use bifurcation theories and computer simulations to systematically investigate the EAD properties and their major determinants. We first made theoretical predictions of the EAD properties using the 1991 Luo and Rudy (LR1) model [28] based on our previous bifurcation theory of EADs [12]. We then carried out computer simulations using physiologically more detailed ventricular AP models [20,22,26,[29][30][31] to verify the theoretical predictions. In computer simulations, we also used parameter sets randomly drawn from assigned intervals so that large parameter ranges are explored to ensure generality of the simulation results. Theoretical and simulation results were compared with experimental results, and potential caveats of the current AP models were discussed.

Action potential models
Computer simulations were carried out in single ventricular myocytes. The governing equation of the transmembrane voltage (V) for the single cell is where I ion is the total ionic current density and I sti the stimulus current density. C m is the membrane capacitance which was set as C m = 1 μF/cm 2 . We simulated six ventricular AP models: the 1991 Luo and Rudy (LR1) guinea pig model [28]; the 1994 Luo and Rudy (LRd) guinea pig model in a modified version [29]; the UCLA (H UCLA ) rabbit model with modifications by Huang et al [22]; the 2004 ten Tusscher et al (TP04) human model [30]; the Grandi et al (GB) human model [26]; and the O'Hara et al (ORd) human model [31].

Simulation methods
The differential equations were numerically solved using a first-order Euler method and the Rush and Larson method [32] for the gating variables with a fixed time step Δt = 0.01 ms.

Control parameters and parameter variations
For each model, a set of control parameters was used. The control parameter set is not the parameter set of the original model but a set we used for the AP to exhibit EADs. The major changes of parameters from the original models are either by increasing the maximum conductance of both I Ca,L (was called slow inward current in the LR1 model, denoted as I si ) and the slow component of the delayed rectifier potassium current (I Ks ) or by increasing the maximum conductance of I Ca,L but decreasing I Ks or I Kr (the rapid component of the delayed rectifier potassium current). The former corresponds to a normal myocyte (the original model) under isoproterenol while the later corresponds to the condition of long QT syndrome with isoproterenol. The specific changes of each model are detailed in S1 Text and the control APs exhibiting EADs are shown in Fig A in S1 Text.
To explore the effects of ionic current conductance on EAD properties in a wide parameter range, we varied the parameters in two ways: 1) We varied one parameter incrementally at a time but kept other parameters in their control values. The fold change of a specific parameter  (Fig 2(b)) Human -15 25~700~1500 [37] (Fig 6) where p c is the control value of p; and 2) We randomly selected parameter sets, with each parameter drawn randomly from a uniform distribution in the interval (0.4p c , 1.6p c ). The rationale for choosing such a parameter interval is that this interval can cover the range from no EAD to many EADs in an AP for most of the parameters in the models simulated in this study (see the figures in S1 Text).

Defining the EAD properties
We investigated three EAD properties in the AP models-amplitude, inter-EAD interval, and latency. As illustrated in Fig 1A, the EAD amplitude (A EAD ) is defined as the difference between the takeoff voltage (V takeoff ) and the peak voltage (V peak ) of an EAD. The inter-EAD interval (T EAD ) is defined as the time interval between the peaks of two consecutive EADs. The latency (L EAD ) is defined as the time interval between the AP upstroke and the time when the 1st EAD takes off. Table 1 summarizes the EAD properties observed experimentally in isolated ventricular myocytes from literature survey , which includes V takeoff , A EAD , T EAD , and L EAD . Based on this literature survey, we found that the V takeoff is always above -50 mV except some of the mouse [53] and guinea pig [48,49] experiments. However, the EADs in the guinea pig experiments [48,49] were induced by injection of a constant inward current, which may behave differently from the ones occurring intrinsically. Therefore, the observed EADs in isolated ventricular myocytes are mainly phase-2 EADs. The maximum amplitude of the phase-2 Definitions of A EAD , T EAD , L EAD , V peak , and V takeoff . B. V peak versus V takeoff from a short segment of sheep Purkinje fiber [1], a short segment of canine Purkinje fiber [1], and a human-induced pluripotent stem cell-derived cardiomyocyte [2]. The dashed straight lines were added as references for the slopes of the negative linear correlations of the data.

Experimentally observed EAD properties
EADs can be as large as 70 mV. It has also been observed in experiments that V peak (thus A EAD ) exhibits a negative linear correlation with V takeoff (Fig 1B). As shown in Table 1, T EAD ranges from 200 ms to 500 ms except some of the mouse experiments [53]. L EAD varies in a wide range, from 30 ms to 3 seconds.

Determinants of EAD amplitude and takeoff potentials
EAD amplitude in the LR1 model: Theoretical predictions and simulation results. We first investigated theoretically the relationship between A EAD and V takeoff based on bifurcation theories of EADs. In a previous theoretical study [12], we showed that phase-2 EADs are caused by a supercritical Hopf bifurcation and terminate at a homoclinic bifurcation in the fast subsystem of the model as the slow variable grows. Namely, in the LR1 model, the Na + current and Ca 2+ current activate and inactivate much faster than the time-dependent K + current, forming the fast subsystem with voltage. The slow variable is the X-gating variable of the timedependent K + current. The equilibrium states (black lines in Fig 2A) of the fast subsystem were obtained by treating X as a parameter, and stability analysis revealed a Hopf bifurcation (the upper black line, labeled as "Q", changes from solid to dashed) leading to oscillations. As X increases, the oscillation amplitude increases, bounded by a bell-shaped envelope plotted as Bifurcations in the fast subsystem when treating the slow subsystem (X) as a parameter. Q, S, and R are the three equilibria in the fast subsystem. The green bell-shaped envelope is the steady state oscillation amplitude (from V takeoff to V peak ) from the Hopf bifurcation point to the homoclinic bifurcation point. The lowest possible takeoff potential is at the homoclinic bifurcation point, which is around -41.5 mV. The red trace (arrows indicate the time course) is an AP from the whole system where X is a variable. Open red circle is the resting state. Note: the bifurcations and simulations shown in this panel and panel B were done without the presence of I Na . B. The steady-state oscillation period from the Hopf bifurcation point to the homoclinic bifurcation point. C. A EAD versus V takeoff . The parameters were randomly drawn from the assigned intervals as described in Methods. For each parameter set, the amplitudes of all EADs in the AP were included. The inset shows representative A EAD versus V takeoff from individual APs distinguished by colors. D. Same as C but for shifted I si kinetics (d 1 and f 1 ). The green data points are for d 1 and f 1 shifted 8 mV toward more negative voltages and the blue ones are for d 1 and f 1 shifted 8 mV toward more positive voltages. The inset shows the corresponding shifts. Panels A and B were replots from Tran et al [12]. The dashed straight lines in C and D are references for slopes. https://doi.org/10.1371/journal.pcbi.1006382.g002 Early afterdepolarizations in ventricular myocytes the thick green line in Fig 2A. The oscillation period also increases as X increases ( Fig 2B). When X increases to a certain value, the takeoff potential of the oscillation meets the unstable equilibrium (labeled as "S") and the oscillation period becomes infinite. When X is greater than this value, no oscillation exists. For the whole system in which X is a variable (note: in the bifurcation analysis, X was treated as a parameter not a variable), as X grows slowly, the system (the red line with arrows indicating time progression) passes slowly through the Hopf bifurcation toward the homoclinic bifurcation with a couple of oscillations before repolarizes to the resting potential. These oscillations are embraced by the bell-shaped envelope. Defining the EAD amplitude as one can approximate V peak by the upper bound of the bell and V takeoff by the lower bound of the bell. If the oscillation is symmetric with respect to the equilibrium point (we call it the quasi-equilibrium state since it is an equilibrium state in the fast subsystem but not in the whole system), then V peak − V QES = V QES − V takeoff with V QES being the voltage of the quasiequilibrium state. In general, the oscillation is not symmetric. We assume an asymmetry factor β and rewrite the V peak and V takeoff relationship as Substituting V peak in Eq 3, we obtain: For a symmetric oscillation (β = 0), Eq 4 predicts a linear relationship between A EAD and V takeoff with a slope -2. If one plots V peak against V takeoff as done in the plots of experimental data shown in Fig 1B, then Eq 4 is rewritten as which predicts a slope -1 for an oscillation symmetric with respect to the equilibrium point. The experimental results shown in Fig 1B correspond to β = 0.3, 1, and 1.3, respectively. A non-zero β indicates that the oscilaltions manifesting as EADs in the experiments are not symmetric to the quasi-equilibrium point. Note that Eq 5 is not obtained from a rigorous derivation but just an empirical observation based on the property of the dual Hopf-homoclinic bifurcation. The Hopf bifurcation shown in Fig 2A is a supercritical Hopf bifurcation, but previous studies [14,15] showed that a subcritical Hopf bifurcation could also be responsible for EAD genesis. In the latter case, Eqs 4 and 5 still hold based on the same reasoning. We then carried out numerical simulations to investigate the relationship between A EAD and V takeoff in the LR1 model. For EADs from the same AP (the points with the same color in the inset of Fig 2C are from the same AP), their amplitudes and takeoff potentials approximately exhibit a linear relationship with slope around -4 (the slope is around -3 if we plot V peak against V takeoff as in Fig 1B). We then plotted A EAD against V takeoff for all the EADs from the randomly drawn parameter sets. The data points are much more scattered, in particular when A EAD is small. The lowest V takeoff is around -34 mV with A EAD~5 5 mV. Based on the bifurcation analysis (Fig 2A), the lowest possible V takeoff in the LR1 model is around -41.5 mV, and the numerical simulation results agree with the prediction of the bifurcation analysis. Since the EADs are caused by reactivation of I Ca,L during the plateau phase, shifting the steady-state activation curve (d 1 ) and steady-state inactivation curve (f 1 ) simultaneously toward more negative or positive voltages results in almost the same shift of the A EAD and V takeoff relationship ( Fig 2D). The maximum EAD amplitude becomes larger when the steady-state curves are shifted toward more negative voltages.
To investigate how A EAD is affected by the maximum conductance and kinetics of different ionic currents, we varied one parameter at a time while maintaining other parameters at their control values. Fig 3A shows A EAD versus the fold of the control G si (α as defined in Eq 2) for EADs in the AP (note: I Ca,L is denoted as I si in the LR1 model and G si is the maximum conductance). When α is smaller than 0.658, no EAD occurs. When α is between 0.658 and 0.74, one EAD appears in the AP. As G si increases, more EADs appear in the AP. When α is greater than 1.2, there are 10 EADs in the AP. Since I si is an inward current, one may also anticipate that increasing I si would progressively increase the EAD amplitude and APD. However, this is not the case. For example, increasing α from 0.66 to 0.74, the EAD amplitude decreases quickly from around 35 mV to around 10 mV and the APD also becomes shorter (see the APs in Fig  3B). As α is slightly above 0.74, a second EAD with a large amplitude (>30 mV) suddenly appears in the AP and also abruptly increases the APD (see the APs in Fig 3C). Note that the parameters for the blue and magenta traces differ slightly (a 0.3% difference), and the two traces are almost identical until at the very end of phase-2 where one repolarizes (exits the basin of attraction of the limit cycle) and the other depolarizes (retains in the basin of the limit cycle) to result in an EAD, a typical all-or-none behavior. The amplitude of this EAD decreases quickly and the APD also decreases as G si increases (see the APs in Fig 3D) until another new EAD appears in the AP. This process repeats as G si increases. Therefore, as shown in Fig 3A,  Early afterdepolarizations in ventricular myocytes the number of EADs in an AP increases progressively with G si , but the amplitude of a specific EAD decreases with G si until a new EAD (all-or-none) suddenly appears in the AP with a maximum amplitude. In other words, the maximum EAD amplitude always occurs when a new EAD appears in the AP. Before a new EAD occurs, the amplitudes of all EADs in an AP are relatively small (such as the EADs in the black AP in Fig 3D). The overall maximum EAD amplitude does not exhibit an apparent change against G si (remained at around 35 mV for α changed from 0.65 to 1.2).
Increasing the conductance of the time-dependent outward K + current (I K ) decreases the number of EADs in the AP, but increases the EAD amplitude (Fig 3E), which is the opposite to the effects of increasing I si . As G K increases, the amplitude of the last EAD reaches maximum immediately before its disappearance from the AP. Slowing the activation and inactivation time constants of I si reduces the number of EADs in an AP but increases the maximum EAD amplitude (Fig 3F). The dependence of EAD amplitude on diastolic interval (DI) is similar to changing a conductance (Fig 3G), i.e., decreasing DI is equivalent to increasing G K .
To more closely link the dual Hopf-homoclinic bifurcation to the behaviors of EAD amplitude shown in Fig 3, we compared the bifurcations of the fast subsystem and the EADs of the whole system in the same way as in Fig 2A. Since changing the maximum conductance of ionic currents not only changes the EAD behavior but also changes the bifurcation of the fast subsystem, to only change the EAD properties but not the bifurcations of the fast subsystem, we changed the time constant of the X-gating variable, τ X . This is because the bifurcation diagram is obtained by treating X as a parameter and thus the bifurcation will remain the same for any τ X . Fig 4A plots A EAD versus the fold change of τ X , showing that increasing τ X increases the number of EADs in the AP and gives rise to the same A EAD behavior as changing the maximum conductance shown in Fig 3. Fig 4B shows three APs with α indicated by the colored arrows in Fig 4A. At the red arrow (the red AP trace in Fig 4B), there are four EADs in the AP. Increasing τ X slightly (the blue arrow in Fig 4A and the blue AP trace in Fig 4B), a new EAD (the 5 th EAD) with a much larger amplitude appears in the AP. As shown in Fig 4C, the new EAD in the blue race takes off right before the homoclinic bifurcation while the red trace just passes the homoclinic bifurcation point during its 4 th EAD. The reason is that X grows slightly slower due to a slightly larger τ X for the blue trace so that the 4 th EAD ends just right before the homoclinic bifurcation, allowing a new EAD to take off. Increasing τ X further quickly reduces the EAD amplitude of the 5 th EAD because its takeoff point moves away from the homoclinic bifurcation point (open arrows in Fig 4D) due to a slower X growth caused by a larger τ X .

EAD amplitude in physiologically more detailed models
To further assess the theoretical predictions and simulation results from the LR1 model, we carried out simulations using physiologically more detailed AP models. Fig 5 shows A EAD versus V takeoff for the five physiologically detailed models we simulated. The negative linear correlation holds roughly for all AP models while the slopes vary from -2 to -5 (corresponding to slopes ranging from -1 to -4 in plots of V peak versus V takeoff ). Shifting the steady-state activation and inactivation curves of I Ca,L results in roughly the same shift in the A EAD and V takeoff relationship, indicating that the I Ca,L reactivation is responsible for EADs in all the models. The maximum EAD amplitude and the lowest detectable takeoff potential also vary largely from model to model. The maximum EAD amplitude (without the shifts in I Ca,L kinetics) of the TP04 model is the smallest (<35 mV) while that of the ORd model is the largest (~90 mV) among the five models. The lowest V takeoff of the TP04 model is the highest (~-20 mV) while that of the H UCLA model is the lowest (~-45 mV) among the five models. Although the EAD amplitude properties vary largely from model to model, they generally agree with the experimental data shown in Fig 1B and Table 1.
We investigated the effects of the maximum conductance of the major ionic currents on A EAD for all five models, including I Ca,L (Fig B in S1 Text), I Ks (Fig C in S1 Text), I Kr (Fig D in  S1 Text), I K1 (Fig E in S1 Text), I NCX (Fig F in S1 Text), I NaK (Fig G in S1 Text), I NaL (Fig H in  S1 Text), and I to (Fig I in S1 Text). The effects of the inactivation time constant of I Ca,L were also shown (Fig J in S1 Text). The general observations are the same as those from the LR1 model, i.e., increasing an inward current increases the number of EADs in the AP, decreases A EAD until a new EAD suddenly appears in the AP at which A EAD becomes maximum. Increasing an outward current decreases the number of EADs in the AP, and increases A EAD and A EAD of the last EAD in the AP reaches maximum before it suddenly disappears from the AP. However, there are some exceptions. For example, increasing the maximum conductance of I to can either promote or suppress EADs (Fig I in S1 Text). In both the LRd model and the ORd model, increasing I to promotes EADs (more number of EADs in the AP), but suppresses EADs in the TP04 model and GB model. For the H UCLA model, increasing I to,s promotes EADs, while increasing I to,f first promotes EADs but then suppresses EADs. Since I to is an outward current, it is generally known that it suppresses EADs [54]. However, recent studies have demonstrated that I to can also promote EADs [45,55]. Whether I to promotes or suppresses EADs depends on its magnitude, speed of inactivation, and its pedestal component. Slow inactivation or large pedestal current tends to suppress EADs [55]. Another exception is I NCX .
Increasing I NCX promotes EADs in all the models except in the TP04 model in which EADs may also be suppressed by increasing I NCX (see Fig F in S1 Text).

Determinants of inter-EAD interval
As shown in Table 1, T EAD in isolated ventricular myocytes is typically from 200 ms to 500 ms. In this section, we investigate T EAD and its determinants in the AP models. We first used the LR1 model for theoretical treatments and then used the physiologically more detailed models to verify the theory.
Inter-EAD interval in the LR1 model: Theoretical predictions and simulation results. Based on Tran et al [12], EADs arises from a Hopf bifurcation in the fast subsystem of the LR1 model. Here we first calculate the period of oscillation at the Hopf bifurcation point analytically and then compare the period of oscillation to the inter-EAD interval. The eigenvalue λ of the Jacobian (Eq 2 in Tran et al [12]) satisfies the following equation: where τ d and τ f are the time constant of the d-gate and the f-gate of I si , respectively; s d and s f are the slopes of the steady-state activation and inactivation curves, respectively. a ¼ @F @V ; b ¼ @F @d , and c ¼ @F @f with F being the total ionic current as a function of voltage (V), gating variable d, and gating variable f (see Tran et al [12] for more details). These quantities are their values at the steady state. At the Hopf bifurcation point, λ 1,2 = ±iω, satisfying (λ − iω)(λ + iω)(λ − λ 3 ) = 0, i.e.,

Early afterdepolarizations in ventricular myocytes
Comparing the coefficients of the Eqs (6) and (7), one obtains: To calculate ω, one needs to define the Hopf bifurcation point. In the LR1 model, the slow subsystem is the X-gate which is treated as a parameter for the stability analysis of the fast subsystem. For a given X value, one can easily obtain the steady-state voltage V QES and thus the steady-state values of all other gating variables. From Eqs (6) and (7), one also has: . This leads to: Defining then at the Hopf bifurcation point, h = 0. Moreover, at the vicinity of the Hopf bifurcation point, λ 1,2 = ε ± iω, satisfying (λ − ε − iω)(λ − ε + iω)(λ − λ 3 ) = 0. Following the same analysis above, one obtains: Since ε<0 before the Hopf bifurcation and ε>0 after the Hopf bifurcation, thus h>0 before the Hopf bifurcation and h<0 after the bifurcation. Therefore, by calculating h using Eq 10 and scanning X from 0 to 1, one can determine the Hopf bifurcation point. After the Hopf bifurcation point is determined, one can then calculate the oscillation period of the limit cycle of the fast subsystem at the Hopf bifurcation using Eq 8. Fig 6 shows the oscillation period versus different parameters using the theoretical approach (open red symbols). The oscillation period decreases slowly with increasing G si (Fig 6A), has almost no change with G K (Fig 6C) and G K1 (Fig 6D), but increases quickly with slowing the inactivation time constant of I si (a 2-fold change in τ f resulted in roughly a 2-fold change in the oscillation period, Fig 6E). For comparison, we carried out simulations of the LR1 model using the same parameters. We plotted all T EAD in an AP as we did for A EAD with the ordering of T EAD in an AP shown in Fig 6B. Similar to A EAD , when a new EAD suddenly appears in or disappears from the AP, the appearing or disappearing T EAD is the longest. As shown in Fig 6A, as G si increases, an EAD suddenly appears in the AP (green arrow in Fig 6A and green AP in Fig 6B) which results in a long T EAD . As G si increases further, this T EAD decreases quickly toward the oscillation period of the limit cycle in the fast subsystem predicted from the theory, and becomes insensitive to G si . This T EAD behavior repeats when a new EAD occurs in the AP as G si increases. The fast decaying phase is a result of the homoclinic bifurcation, in which the oscillation period of the limit cycle decreases quickly as the system is away from the homoclinic bifurcation point (see Fig 2B). In other words, when a new EAD first appears in an AP, it is always the one closest to the homoclinic bifurcation, and thus the A EAD is the largest and T EAD the longest (see Fig 4). As the takeoff of this EAD is further away from the homoclinic point, the corresponding T EAD (e.g., from the blue AP to the magenta AP in Fig 4B) quickly approaches to those of the EADs prior to it, roughly the oscillation period of the limit cycle of the fast subsystem at the Hopf bifurcation. Therefore, changing G si has only a small effect on T EAD except when an EAD appears in or disappears from the AP. Changing the maximum conductance of an outward current exhibits the same behavior (Fig 6C and 6D). However, T EAD increases with the time constant τ f of the I si inactivation gate more sensitively than with the maximum conductance of the ionic currents (Fig 6E), which also agrees with the theoretical prediction. This indicates that besides the dual Hopf-homoclinic bifurcation, I si inactivation and recovery (although τ f is the inactivation time constant, it also determines the recovery of the channel in the LR1 model) is a major parameter determining T EAD .  To more systematically explore the dependence of T EAD on ion channel conductance, we carried out simulations by randomly drawn the maximum conductance of different ionic currents and measured all T EAD . The data is presented as a histogram in Fig 6F. The distribution shows that T EAD is mainly between 250 ms and 500 ms (peak~280 ms), which is the same range seen in Fig 6A-6D (the range is narrower comparing to Fig 6E since the control τ f was used for all the random parameter sets).
Therefore, based on the theoretical predictions and simulation results of the LR1 model, we conclude that T EAD is mainly determined by the time constant of I Ca,L (I si in the LR1 model), which is the oscillation period at the Hopf bifurcation. The change in an ionic current conductance exhibits small effects on T EAD until it causes an EAD to disappear or appear at which T EAD is mainly influenced by the homoclinic bifurcation.

Inter-EAD interval in physiologically more detailed models
The simulation results from the physiologically more detailed models show similar characteristics of T EAD dependence on different parameters (Fig 7). Increasing an inward current decreases T EAD while increasing an outward current (except I to ) increases T EAD . In the LRd and H UCLA models, increasing I to conductance decreases T EAD . This is because I to promotes EADs in these models. Similar to the LR1 model, the change in T EAD is not very sensitive to a change in the maximum conductance of an ionic current until it is close to the transition from two EADs to one EAD (or from one EAD to two EADs) in the AP at which T EAD changes rapidly. The predominant T EAD for the 5 models are roughly: LRd-90 ms; H UCLA -85 ms; TP04-270 ms; ORd-275 ms; and GB-125 ms. As shown in Fig 7, T EAD in the LRd, H UCLA , and GB models are shorter than 200 ms (between 50 ms to 200 ms), while those in the TP04 and ORd models range from 250 ms to 500 ms. In all the models, τ f exhibits a stronger effect on T EAD than the other parameters we explored.
In Fig 7, we only show the first inter-EAD interval versus a specific parameter. To explore wider ranges of T EAD in these models, we used random parameter sets (see Methods) and measured all T EAD in an AP. Fig 8 shows the T EAD distributions for the AP models. The T EAD ranges are: LRd-from 75 ms to 125 ms (peak~90 ms); H UCLA -from 50 ms to 150 ms (peak 85 ms); TP04-from 225 ms to 350 ms (peak~250 ms); ORd-from 200ms to 500 ms (peak 250 ms); and GB-from 100 ms to 200 ms (peak~125 ms). Therefore, the T EAD of the TP04 and ORd model is in the experimentally observed range, while those of the GB, LRd, and H UCLA models are too short comparing to the experimental recordings. Note that the T EAD range using the random parameter sets is similar to the range seen in Fig 7 for each model, indicating that the T EAD range of an AP model is not sensitive to ionic current conductance but mainly determined by the period variation in the dual Hopf-homoclinic bifurcation.
Based on the theoretical predictions that τ f is the main determinant of the inter-EAD interval, we changed τ f functions in the LRd and H UCLA models from upward bell-shaped functions to downward bell-shaped functions (see Fig 8C) to lengthen τ f in the plateau phase, we can effectively shift the T EAD distributions toward the longer periods (red histograms in Fig 8A  and 8B). The inactivation time constants of I Ca,L in the plateau voltage for the TP04 model and the ORd model are much longer, and thus inter-EAD intervals of these two models are also much longer. The GB model has a much shorter inactivation time constant of I Ca,L in the plateau, similar to those of LRd and H UCLA , and thus the inter-EAD interval is also short.

Determinants of EAD latency
As shown in Table 1, L EAD varied in a large range, from less than 100 ms to a couple of seconds. Based on the bifurcation theory of EADs [12,18], for EADs to occur, besides the instability leading to oscillations, the voltage needs to decay into the window voltage range of I Ca,L activation and the LCCs need to be recovered by a certain amount so that there are enough LCCs available for re-opening.
To reveal how L EAD is determined by the ion channel properties, we started with simulations of the LR1 model (Fig 9A). We varied four parameters: G si , G K , G K1 , and τ f . Increasing G si first decreases L EAD quickly and then increases L EAD slowly (red curve in Fig 9A). The longest L EAD occurs when the first EAD appears in the AP (green arrow in Fig 9A). When a new EAD first appears in an AP, its takeoff potential is the lowest which is close to the potential of  the homoclinic bifurcation point (see Fig 2A), and it takes a longer time for the EAD depolarization to occur (the same reason that the T EAD is the longest and A EAD is the largest when a new EAD appears). As G si increases, the takeoff potential is higher and thus L EAD shortens. However, increasing G si also slows the decay of voltage into the window range of LCC reactivation, and thus lengthens L EAD . Increasing G K does the opposite (blue curve in Fig 9A) for the same reasons mentioned for G si . Changing G K1 has no effect since I K1 is almost negligible in early phase-2 of the AP. τ f has a big effect on L EAD , which can vary L EAD in a much wider range than the conductance. τ f affects L EAD by two ways: 1) slowing τ f causes a slower inactivation of I si which delays the voltage decay to the window range; 2) slowing τ f delays recovery of LCCs, which then delays the depolarization of the first EAD.
Similar behaviors occur in all other models (Fig 9B-9F): changing a conductance usually has a small effect on L EAD until it causes the only EAD in the AP to disappear, at which L EAD changes steeply; and the inactivation time constant of I Ca,L is the most sensitive parameter for L EAD . The L EAD varies largely from model to model. In Fig 10, we show L EAD distributions from random parameter sets for all models, and the L EAD ranges are: LR1-from 600 ms to 1000 ms (peak~750 ms); LRd-from 240 ms to 550 ms (peak~300 ms); H UCLA -from 135 ms to 180 ms (peak~150 ms); TP04-from 640 ms to 720 ms (peak~650 ms); ORd-from 360 ms to 600 ms (peak~400 ms); GB-from 240 ms to 500 ms (peak~280 ms).

Discussion
In this study, we used theoretical analyses and computer simulations to investigate EAD properties and their major determinants in AP models of ventricular myocytes. Our major observations and conclusions are summarized and discussed below.

EAD takeoff potential and amplitude
In all the models we simulated, the EAD takeoff potentials are usually above -40 mV (see the 0 mV shift cases in Fig 5), which is in agreement with experimental data from isolated Early afterdepolarizations in ventricular myocytes ventricular myocytes (Table 1). A negative linear correlation between EAD amplitude and takeoff potential, which has been shown in experimental recordings [1,2,9], can be implied from the dual Hopf-homoclinic bifurcation and is shown in simulations of all models. The slopes of the negative linear correlations, the lowest takeoff potentials, and the maximum EAD amplitudes vary substantially from model to model. The ORd model exhibits the largest maximum EAD amplitude and the steepest slope of the linear correlation between EAD amplitude and takeoff potential.
Although EADs are promoted by increasing inward currents or decreasing outward currents, once EADs occur in an AP, increasing the maximum conductance of an inward current or decreasing that of an outward current does not increase the amplitudes of the EADs. Increasing the maximum conductance of an inward current causes more EADs in the AP, and the maximum EAD amplitude occurs when a new EAD appears in the AP. The amplitude of this new EAD then decreases as the conductance increases. Increasing the maximum conductance of an outward current decreases the number of EADs in the AP but increases the EAD amplitude. The maximum EAD amplitude occurs before an EAD disappears from the AP. Based on the bifurcation analysis [12][13][14][15], the EAD amplitude grows as the system evolves from the Hopf bifurcation point to the homoclinic bifurcation point (Fig 2A and Fig 4). This behavior should hold for both supercritical [12,15] and subcritical [14,15] Hopf bifurcation. At the homoclinic bifurcation point, the takeoff potential is the lowest and the EAD amplitude becomes the maximum. Therefore, the amplitude of the last EAD in an AP depends on how far away the takeoff potential is from the homoclinic bifurcation point. Born of a new EAD or death of an existing EAD always occurs when the EAD takes off near the homoclinic bifurcation point.
Note that in our simulations of assessing EAD amplitude and takeoff potential, we used random parameter sets to explore a wide range of parameters for each model. In these simulations, only phase-2 EADs in the ventricular myocyte models were observed. A recent simulation study [15] using the ORd model and the Kurata et al model [56] also showed that only phase-2 EADs could be observed. This indicates that varying maximum ionic current conductance may not produce phase-3 EADs. This agrees with the experimental data that EADs observed in isolated ventricular myocytes are mainly phase-2 EADs (Table 1), while phase-3 EADs are rarely observed in isolated cells, except under Ca 2+ overload [57] or by external current injection [49]. On the other hand, phase-3 EADs are observed in Purkinje fibers [19] and cardiac tissue [58][59][60]. Previous computer simulations showed that phase-3 EADs could occur in single cells with a strong I nsCa under elevated intracellular Ca 2+ concentration [20,21,61] or in tissue with repolarization heterogeneity induced dynamical instabilities [22,58]. Therefore, phase-3 EADs can be caused either by strong Ca 2+ overload in isolated myocytes or by repolarization heterogeneities in tissue with reduced repolarization reserve, while phase-2 EADs are mainly due to reduced repolarization reserve and reactivation of I Ca,L [9,10,18,62,63].
Since phase-2 EADs cannot propagate into PVCs in tissue [19][20][21][22], this raises question on how are EADs linked to arrhythmias under LQTS and many other diseased conditions where Ca 2+ may not be overloaded. In recent studies [22,64], we demonstrated how phase-2 EADs and tissue-scale dynamical instabilities interact to result in PVCs and arrhythmias under LQTS, linking mechanistically phase-2 EADs to arrhythmogenesis.

Inter-EAD interval
Based on the bifurcation analysis, the inter-EAD interval is governed by the period of the limit cycle oscillation between the Hopf bifurcation and the homoclinic bifurcation. Similar to A EAD , the inter-EAD interval increases as the system evolves from the Hopf bifurcation point to the homoclinic bifurcation point. This behavior has been demonstrated in experimental recordings previously [65]. Our theoretical analysis and simulation of the AP models showed that the inter-EAD intervals (except the last one in an AP), which are mainly determined by the period of the limit cycle oscillation at the Hopf bifurcation, is insensitive to the change of maximum ionic current conductance but more sensitive to the recovery of LCCs (Figs 7 and 8, and Eq 7). The T EAD range of an AP model is determined by the oscillation period between the Hopf bifurcation and the homoclinic bifurcation. However, the inter-EAD interval from different models exhibits different ranges, differing several folds. On the other hand, the inter-EAD intervals observed in isolated ventricular myocyte experiments mostly are in the range from 200 ms to 500 ms, irrespective of species (Table 1). In the AP models simulated in this study, the inter-EAD intervals of LR1, TP04, and ORd are in the same range as observed experimentally, but other models exhibit much faster inter-EAD intervals, indicating that caveats may exist in these models. Our simulation indicates that the major caveat may lie in the formulation of I Ca,L , in particular the recovery time of I Ca,L during the plateau phase (see more detailed discussion below).

EAD latency
EAD latency is determined by many factors. In term of biophysics, since EADs are caused by reactivation of LCCs, the voltage needs to decay into the window range for reactivation of LCCs, which depends on the speed of activation of the outward currents (namely, I Ks , I Kr , and I to ) and inactivation of the inward currents (namely I Ca,L ). Then there are enough LCCs recovered for re-opening, which depends on how fast the LCCs recovers. In more general term of nonlinear dynamics as indicated by the bifurcation theory, the voltage and other variables need to enter the basin or the vicinity of the basin of attraction of the limit cycle. This requires not only the voltage but also the other state variables to reach their proper values. For example, in the LR1 model, the X-gating variable needs to grow to a certain value to engage the Hopf bifurcation as shown in Fig 2A. If X grows too slowly, the system may stay at the quasi-equilibrium state for a long time with no oscillations until reaching the Hopf bifurcation point, such as the cases shown in Song et al [51]. Note that transient oscillations around a stable focus can occur before the Hopf bifurcation, and thus EADs can occur before the Hopf bifurcation (see Fig 2A, Fig 4C and 4D in this study, and Figs 2 and 4 in the study by Kügler [14]). Therefore, the EAD latency can be very variable, explaining the experimental observation that EAD latency varies in a wide range, from less than 100 ms to several seconds (Table 1).

Implications to mathematical modeling of cardiac APs
It is obvious that understanding the EAD properties and nonlinear dynamics is of great importance for understanding arrhythmogenesis in cardiac diseases [66], but it also provides important information for cardiac AP modeling. Previous AP modeling has mainly considered AP morphology, APD, as well as APD restitution, but not the EAD properties. For example, the inter-EAD intervals of the guinea pig ventricular myocyte models [11,23,24] are much shorter than what have been observed in isolated guinea pig ventricular myocytes (Table 1). This is also true for the rabbit ventricular myocyte models. As shown in Fig 8A and 8B, we can effectively increase the inter-EAD interval by increasing the inactivation time constant τ f of I Ca,L . However, for both the guinea pig model [67] and the rabbit model [68], the original inactivation time constants were based on experimental measurements. Since in the Hodgkin-Huxley formulation of I Ca,L , the f-gate is a voltage-dependent inactivation gate, but it also governs the recovery of I Ca,L . Therefore, one would conclude that experimentally-based τ f might be a correct constant but the recovery properties of I Ca,L in these models are incorrect, which gives rise to the discrepancy in inter-EAD intervals between mathematical models and experimental measurements. On the other hand, τ f is large in the LR1, TP04, and ORd models, which gives rise to the right recovery times to result in inter-EAD intervals in the ranges as observed in experiments. That also does not mean that I Ca,L models are completely correct in these AP models since we know that the one for the LR1 model gives rise to a too slow inactivation of I Ca,L . Therefore, the EAD properties provide additional important information for AP modeling, which need to be considered in future AP model development.

Limitations
A major limitation is reliable experimental data curation. First, most of the values in Table 1 were estimated from the published figures, which is difficult to be accurate and unbiased. Second, experimental data of EADs recorded from isolated ventricular is not abundant. Moreover, to calculate T EAD , we have to select APs with two or more EADs, which further limited our data sources. Third, most of the experimental plots do not have coordinates but indicated by scale bars. Sometimes, these scale bars may not be correctly labeled due to different reasons. For example, we confirmed with Dr. Li that the time scale bar in Fig 2C of Ref. [33] is 300 ms instead of 150 ms. In computer simulations, we only explored the maximum conductance of the ion currents and the voltage-dependent inactivation time constant of I Ca,L , but it is obvious any parameter that affects repolarization will have an effect on EAD genesis and EAD properties. For example, the time constants of ion channel activation and inactivation, the Ca 2 + -dependent inaction of I Ca,L [69], the intracellular Ca 2+ transient, mitochondrial metabolism and Ca 2+ cycling [61,70], as well as spatial distribution of the ion channels will impact the EAD behaviors, which need to be investigated in future studies. Another limitation of the simulations is that our conclusions may depend on the setting of control parameter and the assigned intervals for random parameters. However, as we show in this study, despite certain distinct difference between models, such as the inter-EAD interval, the general conclusions are not model dependent, and thus, not likely to be affected by the choice of control parameter sets and the assigned intervals. We only investigated the effects of voltage-dependent inactivation of I Ca,L on EAD properties. It is shown that Ca 2+ -dependent inactivation also play important roles in EAD genesis [69], one would anticipate that the Ca 2+ -dependent inactivation might have a large effect on EAD properties. However, the Ca 2+ -dependent inactivation is modeled very differently in different models, we do not have a unique way to alter a parameter to study the effects in the models. For example, in most models, Ca 2+ -dependent inactivation was modeled by an instantaneous function of Ca 2+ , and thus it is not clear how to change the time constant of Ca 2+ -dependent inactivation in these models. One major caveat of the current study is that in the models we simulated, the EADs are caused by reactivation of I Ca,L , however, experimental studies showed EADs can also be caused by spontaneous Ca 2+ release [57,71,72]. In a recent modeling study by Wilson et al [73], the authors showed that spontaneous Ca 2 + oscillation can lead to EADs. We performed the same analyses of this model as we did for other models and showed the results in Fig K in S1 Text. The EAD behaviors and their dependence on the ionic conductances are similar to other models, but the model indeed exhibits some differences from the other models. For example, T EAD exhibits a much wider range (from 100 ms to 1200 ms) than those of other models and the T EAD histogram shows characteristic distributions indicating that there are two mechanisms of EADs involved. However, further investigation is needed to pinpoint the underlying mechanisms of EADs in this model and compare model results with experimental data in future studies.