Bifurcations and Slow-Fast Analysis in a Cardiac Cell Model for Investigation of Early Afterdepolarizations

Abstract

In this study, we teased out the dynamical mechanisms underlying the generation of arrhythmogenic early afterdepolarizations (EADs) in a three-variable model of a mammalian ventricular cell. Based on recently published studies, we consider a 1-fast, 2-slow variable decomposition of the system describing the cellular action potential. We use sweeping techniques, such as the spike-counting method, and bifurcation and continuation methods to identify parametric regions with EADs. We show the existence of isolas of periodic orbits organizing the different EAD patterns and we provide a preliminary classification of our fast–slow decomposition according to the involved dynamical phenomena. This investigation represents a basis for further studies into the organization of EAD patterns in the parameter space and the involved bifurcations.

1. Introduction

In a healthy heart, the sinoatrial node sends out an electrical impulse, which spreads cell to cell throughout the heart, activating all cardiomyocytes to produce an electrical response called the action potential (AP). Upon being stimulated, cardiomyocytes in the lower chambers of the heart, the ventricles, suddenly increase their transmembrane voltage (depolarization). This increase is followed by a small partial voltage decrease (transient repolarization) and a prolonged plateau phase where voltage remains approximately constant. In the final part of the AP, transmembrane voltage decreases (repolarization) while returning to the resting potential level, which is maintained until receiving the next stimulus.

Under some circumstances, the normal sequence of AP phases can be reversed due to inward currents raising the transmembrane voltage during the plateau or repolarization phases of the AP, producing so-called phase-2 or phase-3 early afterdepolarizations (EADs). Drug side effects, ion channel dysfunction or oxidative stress, among others, can lead to the genesis of EADs [1,2,3]. In heart failure and genetic syndromes, EADs have been documented to be an important cause for lethal ventricular arrhythmias [4,5,6], but further knowledge is required to understand the mechanisms underlying EAD generation and the relationship between these cellular abnormalities and the occurrence of arrhythmias at tissue and whole-organ level that could eventually lead to sudden cardiac death.

Computational models of cardiac electrical activity have greatly contributed to shed light into varied cardiac phenomena, including EADs. Multiple models of mammalian ventricular cardiomyocytes have been used in the literature. Some of these are highly detailed, complex models with tens of state variables and hundreds of parameters, while there are other simpler models with just a few variables and parameters. Whereas complex (high-dimensional) models allow for greater realism in reproducing experimental observations and facilitating biophysical interpretation, simple (low-dimensional) models aid in isolating the dynamical mechanisms underlying a particular phenomenon and in performing a comprehensive theoretical study. In the present work, we use a low-dimensional approach based on a reduced three-variable cardiomyocyte model.

In 1991, Luo and Rudy [7] introduced a mathematical model of the membrane AP of a mammalian ventricular cell based on experimental data recently available at that time. This model, differently from subsequent models developed by the same authors, is called “passive” because all ionic concentrations, except for intracellular calcium, remain unchanged rather than varying dynamically over time. Due to its simple formulation and its ability to represent phenomena involving both depolarization and repolarization of ventricular cell, this model, either as such or with some simplifications, has been extensively used. In [8], a reduction in the number of variables of the Luo–Rudy model was proposed to investigate the mechanisms of EADs when decreasing the stimulation frequency. Initially, a model able to produce EADs was considered that contained three ionic currents, namely a fast $Na$ current, an intermediate time-scale $Ca$ current, and a slow K current. Later, the fast $Na$ current was discarded due to it having little effect on EAD generation, as it is activated mostly during the AP upstroke and is practically null during the plateau and repolarization AP phases. In our analysis, we use this reduced model to simplify the number of state variables at its maximum and, thus, facilitate a theoretical study.

From a mathematical point of view, this is a fast–slow dynamical system with multi-timescale phenomena. A fast–slow analysis of a reduced version of the Luo–Rudy model is presented in [8], considering a system with one slow and two or three fast variables. In that paper, the presence of a subcritical Poincaré–Andronov–Hopf bifurcation is shown as a signature of pseudo-plateau bursting. Similar analyses were published based on mathematical neuron models [9,10,11], where the classification of bursting models and the generation of new oscillations (spikes) were related, among others, with Poincaré–Andronov–Hopf and homoclinic bifurcations in the fast subsystem of two variables (dimension 2). While this approach has been successfully used in a variety of cardiac studies to investigate the causes for the presence or absence of EADs, it has recently been shown to fail in explaining the lack of certain types of EADs [12]. Considering the reduced three-variable version of the Luo–Rudy model presented in [8], a 1-fast, 2-slow decomposition is proposed [12,13] to provide more insight into the facilitation or inhibition of EADs as a function of the pacing frequency or pharmacological interventions. In this paper, we use this approach for the fast–slow decomposition and we characterize dynamical behaviors by introducing a sweeping technique, namely the spike-counting method. Besides, by using continuation techniques, we describe some bifurcations of the system and we show, for the first time to the best of our knowledge, the presence of isolas of families of periodic orbits in the reduced Luo–Rudy model.

The paper is organized as follows. In Section 2, the reduced model used to describe the electrical behavior of a mammalian ventricular cardiomyocyte is presented. In Section 3, the dynamical analysis performed to identify parametric regions with EADs and generate associated bifurcation diagrams is described. In Section 4, the discussion and conclusions of the study as well as recapitulating classification figures regarding EAD generation are presented.

2. The Reduced Luo–Rudy Mammalian Ventricular Cell Model

Following a Hodgkin and Huxley formalism [14], Luo and Rudy proposed a mathematical model of a mammalian ventricular cell (LR91) [7]. The rate of change of the transmembrane potential (V) is given by

$$dV/dt=-(1/C_m)(I_i+I_{sti})$$

where $C_m$ is the membrane capacitance, in this study set at $0.5\mathsf{\mu }\mathrm{F}/{\mathrm{cm}}^2$, as in [13], and subsequently varied to investigate its role in facilitating the generation of EADs. In the above equation, $I_{sti}$ is an external stimulus current and $I_i$ is the sum of the ionic currents in the cell. All ionic currents were computed for a membrane area of $1{\mathrm{cm}}^2$. Six ionic currents were defined in the LR91 model: $I_{Na}$, a fast sodium current; $I_{Ca}$, a slow inward current; $I_K$, a time-dependent potassium current; $I_{K1}$, a time-independent potassium current; $I_{Kp}$, a plateau potassium current; and $I_b$, a time-independent background current.

Since here we focus on EAD generation, we follow the approach proposed in [8] and we discard the fast $I_{Na}$ current for the reduced version of the LR91 model. Although some studies have described a role for the sodium current in the generation of EADs [2,3], it should be noted that this current has two components, namely the fast sodium current and the late sodium current. While the late sodium current flows throughout the AP plateau and its involvement in EAD generation has been well documented, the fast sodium current contributes to the AP depolarization and has a more limited contribution to EADs. In the reduced model used in this work, only the fast sodium current is included and, thus, we discard it based on its reduced contribution to the investigated phenomenon. The other two ionic currents in the reduced model are: $I_{Ca}=G_{Ca}·d·f·(V-E_{Ca})$, which is a calcium current with an activation gating variable d and an inactivation gating variable f; and $I_K=G_K·x·x_i·(V-E_K)$, which is a time-dependent potassium current with a time-dependent activation gating variable x and a time-independent inactivation gating variable $x_i$ set to one for simplification. The reversal potential of calcium was set at $E_{Ca}=100\mathrm{mV}$, rather than being time-dependent as in the LR91 model.

The values of the gating variables used to define the ionic currents are obtained as the solution of a coupled system of nonlinear ordinary differential equations (ODEs) of the form $dy/dt=(y_{\infty }\left(V\right)-y)/{\tau }_y\left(V\right)$, where y represents any gating variable, ${\tau }_y$ is its time constant, and $y_{\infty }$ is the steady-state value of y [7].

An additional simplification to the four-variable system (V, d, f, and x) proposed in [8] was later described by Kügler [15]. The gating variable d was replaced with its steady-state function $d_{\infty }\left(V\right)$ and the time-constant functions ${\tau }_f\left(V\right)$ and ${\tau }_x\left(V\right)$ were assumed to be constant, thus being represented by ${\tau }_f$ and ${\tau }_x$.

With all these simplifications, the three-variable model used here [12] was described by the following ODE system:

$$\left\{\begin{array}{c}{C}_m{\displaystyle \frac{dV}{dt}}=-(I_{Ca}+I_K)+I_{sti}\equiv h(V,f,x),\hfill \\ {\displaystyle \frac{df}{dt}}={\displaystyle \frac{f_{\infty }\left(V\right)-f}{{\tau }_f}}\equiv g_1(V,f),\hfill \\ {\displaystyle \frac{dx}{dt}}={\displaystyle \frac{x_{\infty }\left(V\right)-x}{{\tau }_x}}\equiv g_2(V,x),\hfill \end{array}\right.$$

(1)

with the inward ionic calcium current $I_{Ca}$ (with the calcium channel conductance $G_{Ca}$ and the dynamic inactivation variable f) and the outward ionic potassium current $I_K$ (with the potassium channel conductance $G_K$ and the dynamic activation variable x) defined by:

$$I_{Ca}=G_{Ca}·d_{\infty }\left(V\right)·f·(V-E_{Ca});I_K=G_K·x·(V-E_K).$$

(2)

The steady-state functions were given by:

$${\displaystyle f_{\infty }\left(V\right)=\frac{1}{1+e^{\frac{V-E_f}{{\theta }_f}}},x_{\infty }\left(V\right)=\frac{1}{1+e^{\frac{E_x-V}{{\theta }_x}}},d_{\infty }\left(V\right)=\frac{1}{1+e^{\frac{E_d-V}{{\theta }_d}}}.}$$

The fixed parameter values used in this study, unless otherwise stated, were:

$$\begin{array}{ccccc}{C}_m=0.5\mathsf{\mu }{\mathrm{F}/\mathrm{cm}}^2,\hfill & {\tau }_f=80\mathrm{ms},\hfill & {\tau }_x=300\mathrm{ms},\hfill & G_{Ca}=0.025{\mathrm{mS}/\mathrm{cm}}^2,\hfill & E_{Ca}=100\mathrm{mV},\hfill \\ G_K=0.04{\mathrm{mS}/\mathrm{cm}}^2,\hfill & E_K=-80\mathrm{mV},\hfill & E_f=-20\mathrm{mV},\hfill & {\theta }_f=8.6\mathrm{mV},\hfill & E_x=-40\mathrm{mV},\hfill \\ {\theta }_x=5\mathrm{mV},\hfill & E_d=-35\mathrm{mV},\hfill & {\theta }_d=6.24\mathrm{mV}.\hfill \end{array}$$

It should to be noted that, whereas the LR91 model well represents the behavior of a ventricular cardiomyocyte with constant resting membrane potential during Phase 4, the reduced model described above presents a Phase 4 with slowly increasing transmembrane potential. Thus, an external stimulus is required to depolarize the AP in the LR91 model, but the AP depolarizes spontaneously when transmembrane potential reaches a threshold in the case of the reduced model. Consequently, $I_{sti}=0$ in the above equations for the three-variable model.

3. Dynamical Study

We showed that the reduced version of the LR91 model described above presents a behavior comparable to that of typical spiking-bursting activity [9,10,11]. To achieve a better understanding of the model dynamics, we performed a detailed numerical analysis using both a spike-counting technique and numerical continuation to detect the main bifurcations. Numerical simulations were performed using a variable-stepsize variable-order Taylor series method (software TIDES (https://sourceforge.net/projects/tidesodes/) [16,17]), which provides a highly accurate numerical ODE solver, using $TOL={10}^{-12}$ as error tolerance.

3.1. Spike-Counting Sweeping

When EADs are generated, extra spikes can be seen in the plateau of the AP. Figure 1 shows a normal beat with no EADs, a beat with one EAD and a beat with several EADs (more than 2 spikes).

The spike-counting technique consists of detecting the number of spikes of the limit cycles of the system, when these exist [18]. Regions of parameter values with the same number of spikes per beat are represented with the same color. This allows appreciating different bands in the parameter space, which characterize the structure of the model. The detection of the number of spikes per orbit is performed by computing all relative maxima and minima and by counting the spikes where the relative difference between the maxima and minima in the voltage variable is higher than a threshold. Here, the threshold is set at ${10}^{-2}$ so as not to include very small oscillations that occur in action potentials such as the one represented on the right panel of Figure 1. The number of counted spikes is clearly sensitive to the selection of the threshold value, but the global behavior in terms of the occurrence of EADs is generally well captured, as all main voltage changes during the AP plateau are detected.

In our analysis, we used the conductances of the calcium and potassium currents, $G_K$ and $G_{Ca}$, as the main set of parameters to investigate their involvement in EAD occurrence. We also investigated variations in the cell capacitance, $C_m$. This represents a further step from previous studies that considered variations only in $G_{Ca}$ [8] or both $G_{Ca}$ and $G_K$ [12] to determine the presence of EADs. While using a reduced model facilitates theoretical understanding on the mechanisms for EAD generation, future works should be conducted with more complex and more realistic cardiomyocyte models that allow meaningful interpretations regarding the parameters involved in the occurrence of EADs.

Figure 2 shows the results of the spike-counting technique for the biparametric plane ($G_K$, $G_{Ca}$), with the upper left inset presenting a magnification of a specific area. Colors indicate different cellular behaviors. Darkest blue denotes cases where no beats are generated. Slightly lighter blue corresponds to normal beats without EADs, an example of which is illustrated in Figure 3b for the parameter values of Point I. In the region of Point II, the periodic orbit shows one EAD, as exemplified in Figure 3d. When the parameters move to the lightest blue and red regions, more EADs are present, as shown in Figure 3f for Point III and in Figure 4c,f for Points IV and V, respectively. In particular, Point III corresponds to the default values $G_{Ca}=0.025$ and $G_K=0.04$, for which EADs are present. Based on the different cellular behaviors indicated by the color band structure in Figure 2, further analyses were performed for the default value of $G_{Ca}=0.025$ while varying $G_K$ along the horizontal red line, as described in the following sections.

Figure 5 shows additional spike-counting results for the biparametric ($G_K$, $C_m$) plane, with the horizontal red line corresponding to the default value $C_m=0.5$ (and $G_{Ca}=0.025$). In this case as well, color bands can be observed, rendering similar structures to those shown for the ($G_K$, $G_{Ca}$) plane in Figure 2. Based on the results shown in the two figures, bounded regions in the global parameter space can be expected in the characterization of EAD dynamics.

3.2. Fast–Slow Analysis

The reduced version of the LR91 model used in this study is a fast–slow system, since its state variables change at different time scales. When analyzing a fast–slow system, it is of major relevance to investigate the bifurcation diagrams of the limit cases (when the parameters responsible for the difference in time scales are considered to be zero) [10]. This approach uses geometric singular perturbation methods and allows (partially) explaining the dynamics when such parameters are small enough (see [19] for a review and [9] for the basic theory and for classification of bursting mechanisms).

The default values of the time constant parameters characterizing the state variables f and x in the reduced LR91 model are ${\tau }_f=80\mathrm{ms}$ and ${\tau }_x=300\mathrm{ms}$. An estimate of the time constant for the state variable V can be obtained by ${\tau }_V=\frac{C_m}{G_{Ca}+G_k}=7.7\mathrm{ms}$ [12]. This means that V is the fastest state variable and x is the slowest one, while the state variable f can be considered as a fast or a slow one. In [8], the four-variable (V, d, f, and x) system was decomposed into a fast three-variable (V, d, and f) subsystem and a slow one-variable (x) subsystem, thus including f as a fast variable. In [15], a simpler three-variable (V, f, x) model was considered, as d was replaced with its steady-state function $d_{\infty }\left(V\right)$. The fast subsystem was defined to contain two variables (V and f), thus again taking f as a fast state variable. On this basis, a fast two-variable and a slow one-variable subsystems were analyzed and their information was combined to characterize the dynamics of the full three-variable system. By fixing x, the model analysis provides two invariant objects: a curve of equilibrium points and a manifold of limit cycles. The main issue with this approach is that it does not explain some changes in the system, as illustrated in [12]. To solve this issue, some studies have decomposed the same three-variable system into a fast one-variable (V) and a slow two-variable (f and x) subsystems [12,13]. A systematic analysis of the geometric singular perturbation methods for 1-fast, 2-slow variables is given in [20,21].

In this study, we considered, on the one hand, the reduced LR91 model (Equation (1)) using the fast time scale $\tau =t/C_m$:

$${\displaystyle \frac{dV}{d\tau }}=h(V,f,x),{\displaystyle \frac{df}{d\tau }}=\epsilon ·g_1(V,f),{\displaystyle \frac{dx}{d\tau }}=\epsilon ·g_2(V,x).$$

(3)

where $\epsilon$ is a parameter that takes small values ($C_m/{\tau }_f,C_m/{\tau }_x$) and $h(V,f,x)=-(I_{Ca}+I_K)$, as no external stimulus was considered ($I_{sti}=0$), as in [12,13]. When $\epsilon$ decreases to zero ($\epsilon ⟶0$), the system described in Equation (3) defines orbits that converge during fast dynamics to solutions of the fast subsystem or layer equations given by:

$${\displaystyle \frac{dV}{d\tau }}=h(V,f,x),{\displaystyle \frac{df}{d\tau }}=0,{\displaystyle \frac{dx}{d\tau }}=0.$$

(4)

On the other hand, when the orbits have slow dynamics, the fast variable moves so rapidly that it can be considered to have already reached steady-state. Thus, changing to the time scale $\tilde{t}=\tau /\epsilon$ and taking the limit case, the orbits converge to solutions of the differential algebraic equation (DAE) system, called the slow-flow system, given by:

$$0=h(V,f,x),{\displaystyle \frac{df}{d\tilde{t}}}=g_1(V,f),{\displaystyle \frac{dx}{d\tilde{t}}}=g_2(V,x).$$

(5)

The solutions of Equation (5) evolve on the manifold given by $h(V,f,x)=0$, which is called the critical manifold and is denoted by ${\mathcal{M}}_{\mathrm{critical}}$. Besides, this gives the manifold where the equilibria of the fast subsystem are. It follows from Fenichel theory [22] that this manifold perturbs to invariant manifolds that exist for small enough $\epsilon$ in the full system.

In the following, a detailed investigation of the geometry of the bursting orbit dynamics for the reduced LR91 model, and in particular the orbits for the set of parameters corresponding to Points I–V in Figure 2, is performed. Taking the model equations defined by Equation (1), the critical manifold is given by the cubic-shaped surface:

$${\mathcal{M}}_{\mathrm{critical}}:=\left\{(V,f,x)\left|h\right(V,f,x)=0\right\}\to {\mathcal{M}}_{\mathrm{critical}}=\left\{(V,f,x)|f=-\frac{G_{K}x(V-E_K)}{G_{Ca}{d}_{\infty }\left(V\right)(V-E_{Ca})}\right\}.$$

(6)

Figure 3 shows the outer sheets of the surface, which are attracting ones, and the middle sheet, which is a repelling one. The different sheets are separated by curves $F^{\pm }$ corresponding to fold bifurcations of the fast subsystem:

$$F^{\pm }=\left\{(V,f,x)\in {\mathcal{M}}_{\mathrm{critical}}|\frac{\partial h(V,f,x)}{\partial V}=0,\frac{{\partial }^{2}h(V,f,x)}{{\partial }^{2}V}\ne 0\right\}.$$

(7)

From the analysis of limit cases and Fenichel theory, the solutions of the global system were found to evolve in the slow epochs close to one of the attracting sheets of the critical manifold. The evolution on the attracting sheets was found to give the stable hyperpolarized and stable depolarized steady-states. Close to the fold bifurcations, the orbits quickly fell down to the other attracting sheet. The rapid transitions between the attracting sheets were approximated by solutions of the fast one-variable subsystem in Equation (4).

As described above, the equations of the slow motion on the critical manifold are given by Equation (5), where differential equations describe the motions of the variables f and x, whereas the fast variable V is implicitly described by the first algebraic equation in Equation (5). Through several manipulations (see [12,13,19] for details), the slow subsystem can be transformed into the so-called desingularized system:

$$\begin{array}{c}{\displaystyle \frac{dV}{d{t}_d}}={\displaystyle \frac{\partial h}{\partial f}}{g}_1+{\displaystyle \frac{\partial h}{\partial x}}{g}_2,\hfill \\ {\displaystyle \frac{dx}{d{t}_d}}=-{\displaystyle \frac{\partial h}{\partial V}}{g}_2.\hfill \end{array}$$

(8)

Folded singularities are equilibria of the desingularized system only, but not of the slow subsystem in Equation (5) or the initial model in Equation (1). Thus, they lie on a fold curve and satisfy:

$$FN=\left\{(V,f,x)\in F^{\pm }|{\displaystyle \frac{\partial h}{\partial f}}{g}_1+{\displaystyle \frac{\partial h}{\partial x}}{g}_2=0\right\}.$$

(9)

A study of their linear stability suggested they are folded-node equilibria in this system. These special points give a route to the solutions of the slow subsystem to cross the folds from an attracting sheet and to move some time on the repelling sheet. Therefore, these points may generate the so-called singular canards (see [19,23] for details).

The existence of canard orbits have relevant consequences in many different systems, as they allow uncovering mechanisms of sudden changes. For instance, they were linked to the spike-adding phenomena in the Hindmarsh-Rose neuron model [24,25]. In a reduced LR91 model, canard orbits were shown to organize the first EADs generated on the orbits in phase space [12]. This is well explained in references [12,13] and therefore here we focus on the different geometric characterizations of the periodic orbits of the system showing where EADs are produced by the canard orbits and where different phenomena are also present.

In Figure 3, the fast–slow geometry is described for the orbits with labels I, II, and III in Figure 2 on the parametric line $G_{Ca}=0.025$. The critical manifold is shown in orange, with the attracting and repelling sheets separated by the fold lines, $F^{+}$ and $F^{-}$, in red and with the periodic orbits in blue. The temporal evolution of the variable V of the periodic orbit is shown on the right panels. It can be observed that in Case I for $G_K=0.05$ the orbit presents a normal beat without any EADs (Figure 3a,b). In this case, the orbit presents the “nominal behavior”, remaining most of the time on the critical manifold of the slow motion and fast transitioning between the attracting sheets close to the fold lines. Note that the orbit is far from the folded-node singular point in magenta and, therefore, there is no canard orbit. By moving forward on the line $G_{Ca}=0.025$ towards Point II corresponding to $G_K=0.045$, it can be observed that the orbit has one EAD (Figure 3c,d). In this case, the orbit passes near the folded-node singular point and this allows the orbit to make a loop on the repelling sheet before progressing to the other attracting sheet (see [13]). By further moving on the parametric line towards Point III corresponding to $G_K=0.04$, more EADs are generated by the maximal canards (Figure 3e,f). Maximal canards are given by the intersections of the depolarized attracting and the repelling sheets of the slow manifold of the system and extended by the flow (for more details, see [12,13,26]). The equilibrium point, which is a saddle-node of type (1,2) with eigenvalues ${\lambda }_1=-3.91295$, ${\lambda }_2=0.0928806$, and ${\lambda }_3=4.27207$, is attracting the orbit along the stable manifold $W^s\left(eq\right)$ of the equilibrium (in green) at the same time that the canard orbits generate new EADs until the orbit progresses to the other attracting sheet (see inset in Figure 3e). This phenomenon has been previously observed in other 1-fast, 2-slow variables systems [27,28,29].

Figure 4 shows how further variations in the bifurcation parameter $G_k$ lead to higher numbers of EADs as the orbit approaches the equilibrium point shown in black. The two orbits presented in Figure 4 contain a large number of EADs, particularly the second one. The mathematical mechanism explaining such a large number of EADs is twofold: the described canard phenomenon that gives rise to EADs when the orbit passes near the folded-node and the stable manifold $W^s\left(eq\right)$ of the equilibrium that increasingly pulls the orbit closer to that point. On the orbit of Point IV corresponding to $G_K=0.035$, there is a saddle-focus equilibria of type (1, 2) with eigenvalues ${\lambda }_1=-0.00436926$, ${\lambda }_2=0.00861114+0.035771i$, and ${\lambda }_3=0.00861114-0.035771i$. It can be observed in Figure 4 how such an orbit oscillates due to the canard orbits around the stable manifold of the equilibria $W^s\left(eq\right)$ until reaching a point close to the equilibrium point (Figure 4a,b). Consequently, the orbit remains more time active as the orbit approaches the equilibrium and it subsequently spirals outward (last spikes of the orbit). This can be appreciated in the two-dimensional local unstable manifold $W_{loc}^u\left(eq\right)$ of the equilibrium shown in Figure 4. In a previous study, this phenomenon was observed as a first approximation using the unstable linear subspace $E^u\left(eq\right)$ of the equilibrium [15]. In addition, similar observations hold for the last investigated point, i.e., Point V corresponding to $G_K=0.03345$ (Figure 4d,e). As before, the first oscillations are due to the canard phenomena and the fact of passing near the folded-node singularity, while subsequent oscillations can be explained by the orbit approaching the saddle-focus equilibria of type (1,2) with eigenvalues ${\lambda }_1=-0.0040745$, ${\lambda }_2=0.00137656+0.0399119i$ and ${\lambda }_3=0.00137656-0.0399119i$. It can be noted that in this case the ratio $|{\lambda }_1|/\mathrm{Re}{\lambda }_{2,3}=2.95991$, which makes the attracting dynamics strong and prolongs the approach towards the equilibrium point. When the orbit is already close enough to the equilibrium point, the unstable manifold takes control of its escape dynamics and the orbit follows this manifold. This is illustrated in the inset of Figure 4d, which shows how the orbit spirals the two-dimensional local unstable manifold $W_{loc}^u\left(eq\right)$. The small variations in the plateau phase of the AP closely resemble the oscillations present in pseudo-plateau bursting, whose dynamics has also been studied for 1-fast, 2-slow variables systems [21].

As the parameter $G_K$ takes decreasing values, the orbit approaches more and more the saddle-focus equilibria. As shown above, more and more small spikes are then generated. The fact of setting a threshold on voltage when counting the number of spikes may render the algorithm unable to detect a small spike, as shown in Figure 6, which can be additionally influenced by the enlargement of the orbit approaching the saddle-focus equilibria. This may generate transitions between areas of different colors, which would seem to indicate a change in the number of spikes by one or more spikes even if this is not actually the case.

3.3. Isolas of Periodic Orbits

To investigate bifurcations on the parametric line $G_{Ca}=0.025$, numerical analysis with continuation techniques was performed using the software AUTO [30,31]. The main focus of this study was to explore the presence of isolas of periodic orbits, that is, simple closed curves of families of periodic orbits. These isolas have been described for other models, such as the Koper model of chemical reactions with mixed-mode oscillations (see Figure 19 in [19]). For the reduced LR91 model used in this study, some bifurcations are identified in [12] but no isolas are shown.

Figure 7 shows the bifurcation diagram obtained with the software AUTO using $G_k$ as the continuation parameter on the line $G_{Ca}=0.025$. The bifurcation diagram shows the maximum (the highest peak of the orbit) and minimum (resting membrane potential) values of the transmembrane voltage variable V for different periodic orbits corresponding to a range of given values of the continuation parameter. The continuation of the equilibria (black thick line) and limit cycles (color lines) are shown too. Continuous lines are stable invariants while discontinuous lines are unstable ones. In the figure, presenting a large parametric interval $G_k\in [0.02,0.46]$, two Poincaré–Andronov–Hopf bifurcations can be observed, a subcritical one on the left part and a supercritical one on the right part. As $G_k$ decreases, the stable periodic orbits generated at the supercritical Hopf bifurcation move to the left towards the subinterval for $G_k$ with associated EADs.

In the subinterval with several EADs, shown in Figure 8b, families of periodic orbits can be observed, which are plotted in different colors. These families are organized in isolas, that is, closed family curves of periodic orbits with different numbers of oscillations in the form of EADs. The different isolas have a section formed by stable limit cycles placed mainly on the top (and bottom, as it is the maxima and minima) of the isola. These stable periodic orbits lose their stability at fold bifurcations (or saddle-node bifurcations of limit cycles), as indicated in [12]. However, the stable line is in fact a discontinuous one, as it is formed by the top subintervals of stable periodic orbits of the different isolas, and not a continuous line with several bifurcations. The unstable parts of the isolas shown in the figure evolve with similar maximum and minimum values, which explains why they are so closely represented in Figure 8b, while they are formed by different limit cycles, as can be appreciated in Figure 8a. Besides, there are fold and period-doubling bifurcations on the bifurcated periodic orbits from the subcritical Hopf bifurcation.

Two of the isolas in Figure 8b are represented separately in Figure 8c only for the maximum value of voltage. Now, the fold bifurcations are clearly seen as the limit of the stable periodic orbits. The temporal evolution of voltage along the hearbeat is represented for the two parametric values marked with a black dot in each of the two isolas, each showing a different number of EADs. The stable part of the isolas is located on the flat top segment, which corresponds to the observed limit cycles. These results confirm that the reduced LR91 model exhibits isolas of periodic orbits, as reported for other fast–slow models in the literature [19]. This can be the basis for further studies into the organization of EAD patterns in the parameter space and the involved bifurcations.

4. Discussion, Conclusions, and Future Work

In this study, we investigated the dynamical mechanisms for EAD generation in the reduced LR91 mammalian ventricular cell model. As in recent studies [12,13], we considered a fast–slow decomposition of the system with 1-fast and 2-slow variables rather than 2-fast and 1-slow variables [8], and we analyzed the influence of maximal canard orbits. We used sweeping techniques (the spike-counting method) as well as continuation techniques. The former allowed us to identify different parametric regions with EADs (see a summary in Figure 9), whereas the latter was used to generate bifurcation diagrams for the mentioned parametric regions. We showed the existence of isolas of periodic orbits and we performed a preliminary classification of the fast–slow decomposition in various cases according to the involved dynamical phenomena.

As a summary, based on the dynamics on the critical manifold and near the unstable equilibrium point, several regions in the ($G_K$, $G_{Ca}$) and ($G_K$, $C_m$) parametric planes associated with different dynamical behaviors can be identified. Specifically, regarding EADs, the first region corresponds to cases where the orbit has an EAD due to a maximal canard orbit. As more canards take part in the process, more EADs are generated. When the orbit approaches the equilibrium, this gives rise to additional oscillations. When the equilibrium point is of saddle-focus type, the orbit gets closer and closer to the equilibrium and spirals following its unstable manifold. These mechanisms create EADs corresponding to small oscillations of voltage during the course of the AP, commonly in association with AP prolongation. Figure 9 recapitulates our classification of the dynamical behaviors of the orbits in the ($G_K$, $G_{Ca}$) and ($G_K$, $C_m$) parametric planes: white is associated with no activity, that is, the equilibria is an attractor; blue corresponds to normal beats; yellow denotes EADs created by the maximal canards; red corresponds to the attracting behavior of the equilibria on the stable manifold direction with AP prolongation; and maroon is for very elongated APs remaining close to the saddle-focus equilibria. Color regions, associated with active cell status, are delimited by two Hopf bifurcations that make the equilibrium point unstable in the middle parametric interval. As can be noted from the two panels in Figure 9, results for both ($G_K$, $G_{Ca}$) and ($G_K$, $C_m$) parametric planes render bands that allow making a parametric description of the areas in the parameter space that lead to EAD generation. This sets the basis to investigate actions able to move the system, i.e., the mammalian ventricular cell, to parameter regions far from areas prone to arrhythmogenic EADs. Part of our future research in this line focuses on studying the effects of including an external stimulus in this model and the dynamics of other, more complex and realistic cardiomyocyte models. This would allow more meaningful interpretations regarding the parameters involved in the occurrence of EAD.