Reduced Models of Cardiomyocytes Excitability: Comparing Karma and FitzHugh–Nagumo

Since Noble adapted in 1962 the model of Hodgkin and Huxley to fit Purkinje fibres, the refinement of models for cardiomyocytes has continued. Most of these models are high-dimensional systems of coupled equations so that the possible mathematical analysis is quite limited, even numerically. This has inspired the development of reduced, phenomenological models that preserve qualitatively the main feature of cardiomyocyte’s dynamics. In this paper, we present a systematic comparison of the dynamics between two notable low-dimensional models, the FitzHugh–Nagumo model (FitzHugh in Bull Math Biophys 17:257–269, 1955, J Gen Physiol 43:867–896, 1960, Biophys J 1:445–466, 1961) as a prototype of excitable behaviour and a polynomial version of the Karma model (Karma in Phys Rev Lett 71(7):16, 1993, Chaos 4:461, 1994) which is specifically developed to fit cardiomyocyte’s behaviour well. We start by introducing the models and considering their pure ODE versions. We analyse the ODEs employing the main ideas and steps used in the setting of geometric singular perturbation theory. Next, we turn to the spatially extended models, where we focus on travelling wave solutions in 1D. Finally, we perform numerical simulations of the 1D PDE Karma model varying model parameters in order to systematically investigate the impact on wave propagation velocity and shape. In summary, our study provides a reference regarding key similarities as well as key differences of the two models.


Introduction
Excitability is a fundamental property of cardiomyocytes that generally plays a very important role in biology and medicine. Indeed, phenomena such as exchange of information between neurons, muscle contraction including cardiac arrhythmias, the emergence of organised patterns during development, all emerge as a result of excitable behaviour.
Hodgkin and Huxley proposed the first ionic model to represent excitable behaviour, namely action potentials in a nerve fibre [24]. This model is given by a 4-dimensional system of differential equations with one variable for the voltage and three gating variables for the ion channels. By adapting those equations to cardiac cells Noble developed in [37] a similar model for Purkinje cells opening up a new line in mathematical modelling focused on the heart. Since then there have been many different models either including different currents or ionic pumps (see for example further models developed by Noble et al. [34,10]) and also models for other parts of the heart instead of Purkinje fibres like e.g. ventricular cells in the Beeler-Reuter model [4]. To find a more extensive list of cardiac cell models see [16].
All the models mentioned above are quite complex with at least four variables and highly nonlinear. That makes analytical results often impossible, or at least extremely cumbersome. Furthermore, many models are even too complex for detailed numerical analysis and simulation for many parameters. This is why already in 1962 FitzHugh developed a simplified model of the Hodgkin-Huxley equations as a variation of the van der Pol oscillator [41,42], which focused on capturing the excitable properties of the system (see [17,18,19]). Almost at the same time Nagumo published the corresponding electrical circuit in [36]. Their model reduces the four-dimensional system of Hodgkin and Huxley into two equations separating the fast time-scale of the excitation and the slow time-scale of recovery. The equations are given by where v corresponds to the voltage and w to a slow gating variable, a, b are model parameters, D is the diffusion coefficient and I is an external current. Although its representation of nerve fibres or cardiac cells is not as precise as it would be with a more complex model, the FitzHugh-Nagumo (FHN) model has been studied extensively in the literature and is one of the prototype models for excitable media due to its simplicity.
In the same spirit as FitzHugh, Karma developed in [28,29] a reduced version of the Noble model preserving the fast-slow structure of FitzHugh-Nagumo. Here we present a systematic analysis of the Karma model comparing it to the FitzHugh-Nagumo model given that both models have been extensively used to model the behaviour of cardiomyocytes. Some additional references can be found in [3,35] for the Karma model and [2,5,38] for the FitzHugh-Nagumo model. In section 2 we show under what assumptions the 1993 [28] and 1994 [29] versions of the Karma model are equivalent. In section 3 we present a systematic comparison of the FitzHugh model eqrefeq:FHN and the Karma model as defined in section 2. We conclude our comparison in section 4 with numerical simulations of the full PDE systems with a focus on the Karma model.

The Karma model
As mentioned above, the Karma model introduced in [28] in 1993 is a twovariable model involving one fast and one slow variable similarly to the FitzHugh-Nagumo model but incorporating additionally important dynamic features of the Noble model for cardiac cells. The Karma model equations read (2.1) with the function θ(x) = max{x, 0}, which is also known as a rectifier and is commonly employed currently as a rectified linear unit (ReLU) in machine learning. The parameter 0 < n B < 1 controls the position of the excitable wave and the parameter M 1 controls the insensitivity of the excitable wave velocity with respect to the slow gating variable n. Furthermore the constant E * = 1.5415 has been fitted such that for some E where f E is the right-hand side of the first equation without the diffusion term.
In a follow-up paper [29] in 1994 the model was formalized in a slightly more general way as follows (2.3) One may view τ E and τ n as defining the scales for the reaction terms at which E and n respectively evolve. We therefore define ε = τ E /τ n as a single parameter separating the time scales. Next, to make sure there are exactly two stable equilibria for n fixed (corresponding to the depolarized and polarized states) a common choice [28,29] for the reaction function h is h(E) = (1 − tanh(E − 3)) E 2 2 (2.4) and the parameter E * is kept as defined above. Alternatively, a common suggestion [28,29] is a function of the form h(E) = E 2 − δE 3 which we will treat in more detail later.
To fully define the model we still have to specify the restitution function R(n) and dispersion function D(n). The former is responsible for the length A of an action potential after a diastolic or rest interval of length D. The latter function defines the relation between the dispersion velocity c of a pulse with respect to the previous diastolic interval. In theory, both functions can be chosen to fit arbitrary restitution and dispersion curves of the system to be modelled. where c(D) is the dispersion curve and C is a function that can be fitted numerically by a third order polynomial. Karma Furthermore we rescale the gating variableñ = n/n B and use the parameter transformation n B = 1 − e −Re . Again, after dropping the tildes we have which differs from the 1993 model (2.1) only in the diffusion parameter. Since the model is non-dimensional we can assume without loss of generality that τ E = 1. In particular τ E is independent of ε so we have again a slightly more general formulation of the same model. By choosing γ = ε 2 we obtain the 1993 model.
In the remainder of this paper we use the simpler form of the reaction function mentioned above. To stay as close to the function used by Karma as possible we have chosen the reaction function which is the third-order Taylor expansion of (2.4) at E = 3. Figure 1: Reaction function used in [28,29] (orange) together with the function h(E) we are going to analyse in this paper (blue).
Due to the change in the reaction function we have to adapt the parameter E * such that condition (2.2) is still satisfied, which yields E * = 1.5.
Summarizing, in the reminder of the paper we are analysing the following model equations for E * = 1.5 and δ = 0.25 where we used a mixed form of the scalings in [28] and [29].

FitzHugh-Nagumo and Karma model
In this section we proceed analysing and comparing the Karma model (2.12) to the classical FHN system (1.1). We note that both models are two-dimensional systems with a clear fast-slow structure and a diffusive term for the fast variable representing the voltage. A concise introduction into the mathematical theory we will be applying is given in Appendix A. In this paper, we are going to focus on illustrating and extracting the main geometric and analytical insights needed in the proofs of different types of dynamics, which makes it more transparent, where the similarities and differences between the two models are, and how to interpret these differences biologically.

Pure ordinary differential equations (ODE) models
We start by comparing the simplified version of both models by considering the pure ODE models, i.e. we set the diffusion coefficients equal 0. Hence we are working with the equations with 0 < b < 1 and 1 − 2 3 b < a < 1.

FitzHugh-Nagumo
The FitzHugh-Nagumo model has been analyzed extensively in the literature due to its simplicity and generality. In this paper we choose the parameter values a = 0.7 and b = 0.8 as standard configuration for cardiac cells following [19] such that for I = 0 the unique equilibrium is stable corresponding to the polarized state. For completeness we present now a short overview over the most important steps of the analysis of the ODE system when I = 0 by exploited the time-scale separation in the system. For more extensive proofs and deeper analysis of the FitzHugh-Nagumo model see [26,27,39,40].
In the FitzHugh-Nagumo model (3.2) the flow is always controlled by the third order critical manifold as we can observe in its phase plane in Figure 2. The manifold can be divided by its extrema into three branches, where the outer ones are attracting and the middle branch is repelling, therefore the flow away from the critical manifold will approach fast to one of the outer branches. When I = 0, orbits on the middle or close to the right branch follow the slow flow upwards towards the maximum where they jump fast towards the left branch. Once close to the left branch every orbit will finally converge towards the sable equilibrium. Figure 2: Phase plane of the FitzHugh-Nagumo system for I = 0. We can see the critical manifold (orange) and the w-nullcline (dashed) as well as a prototypical orbit (yellow) converging to the global equilibrium (black).
The next theorems formalize the main results. The proof is based on the decomposition of the different time scales when ε is small. For this reason we first look in Theorem 3.1 at the singular limit separating the analyse of the layer problem and the reduced system before constructing the candidate orbits. Finally in Theorem 3.2 we perturbed the candidate orbits for ε > 0.
Theorem 3.1. In the singular limit ε = 0 of the FitzHugh-Nagumo equations (3.2) with I = 0 we have a unique stable equilibrium and the third order critical manifold Every candidate orbit can be constructed as concatenation of fast segments converging to one of the outer branches of C 0 which are attracting and slow segments on the critical manifold. Eventually all orbits converge to the fixed point.
Sketch of the Proof. The layer problem defines a one-dimensional system with equilibria given by (3.3) and their Jacobian is negative on the outer branches, positive in the middle branch and 0 on both extrema. The non-zero eigenvalues divide the critical manifold into 2 attracting outer branches and a repelling inner branch while both non-hyperbolic points satisfy the conditions of generic folds.
On the other hand, the slow flow is determined by the w-nullcline which is a linear function. On the left of this line the flow moves downwards while it moves upwards on the right. Last, the w-nullcline crosses C 0 exactly once, for I = 0 the intersection occurs on the left branch of the critical manifold which results in a stable equilibrium as previously mentioned. The candidate orbits can now be constructed following first the fast flow towards one of the attracting branches. On the middle and right branch of the critical manifold we converge by the slow flow to the maximum where we switch again to a fast fiber connecting to the left branch. Finally, on the left branch the slow flow converges to the equilibrium.
To finish our analysis we show that the candidate orbits we constructed when ε = 0 correspond in fact to solutions of the FitzHugh-Nagumo model for ε > 0.
Theorem 3.2. The candidate orbits found in the singular limit ε = 0 of equations (3.2) when I = 0 can be perturbed to solution curves of the full system with ε > 0.
Sketch of the Proof. We have already seen that, excluding both extrema, that the critical manifold is normally hyperbolic. Therefore, choosing an appropriate compact subset of C 0 , all the conditions of Fenichel's theorems A.1-A.3 are satisfied, which means that the the slow flow on the critical manifold and switches from the fast fibers to the slow flow persist under a smooth perturbation. Finally, we need to transform a neighbourhood of the extrema into the normal form of a generic fold point to apply geometric blow-up as introduced in [30] (see also Appendix A.3). This method provides the persistence of the switches at the maximum and minimum.

Karma: No external current
To understand the Karma model equations (3.1) we will now perform a similar analysis. We will show that the dynamics of the Karma model are similar to FitzHugh-Nagumo since again the system is controlled by the critical manifold presenting a similar shape as shown in Figure 3. As before, we shall indicate the main geometric steps of each proof; see also the appendix for more background on the geometric view via geometric singular perturbation theory.
As before we have two attracting branches separated by a repelling one and exactly one stable equilibrium. An arbitrary orbit will either approach the right branch and then slowly ascend towards the fold point, where it jumps towards the left branch or it approaches directly the left branch where it slowly converges to the stable equilibrium at the origin. In contrast to FHN, the Karma model shows in addition to the stable equilibrium two further unstable fixed points. In general, these points do not affect the overall dynamics, however, the system (3.1) has not only two additional fixed points, which do not converge to the stable equilibrium but also a slow singular heteroclinic orbit between them. We can see the critical manifold (orange) and the w-nullcline (dashed) as well as the global equilibria (black) and a prototypical orbit (yellow) converging to the stable fixed point (0,0).
To formally analyse the dynamics, we want to exploit the different time scales as we did for FitzHugh-Nagumo and therefore consider first the limit ε = 0 analysing the layer and the reduced problem separately.
Remark. In contrast to FHN the Karma model is continuous but not smooth due to the rectifier in the second equation. Although some of the analysis techniques used for FitzHugh-Nagumo have to be modified or extended, the existence and uniqueness of solutions is still guaranteed by the Picard-Lindelöff Theorem.
Theorem 3.3. In the singular limit ε = 0 of the Karma equations (3.1) with I = 0 we have a stable, an unstable and a saddle equilibrium and the critical manifold .
Every candidate orbit can be constructed as concatenation of fast segments converging to one of the outer branches of C 0 , which are attracting and slow segments on the critical manifold. An orbit will either eventually converge to the stable equilibrium at the origin or it is one of the two unstable fixed points or a unique heteroclinic orbit between them.
Sketch of the Proof. Layer problem: Like in FHN we have the one-dimensional fast subsystem where n is a parameter. In this subsystem E = 0 is always an equilibrium and, depending on n we have either two further equilibrium points for n < 1, exactly one for n = 1, or no further equilibria when n > 1. We can calculate the derivative on these points and get which for (E, n) ∈ C 0 is negative on the outer branches, positive in the middle branch and 0 only at p F = (2, 1).
Therefore the critical manifold is normally hyperbolic everywhere except at p F . It is straightforward to check that p F satisfies all the conditions of a generic fold point.

Reduced problem:
The slow flow has a piece-wise linear nullcline which intersects C 0 exactly three times as shown in Figure 3: twice on the E-axis and once with E, n > 0. Therefore, we have three global fixed points: a stable equilibrium at the origin, a saddle at the intersection of the unstable branch of C 0 and the E-axis and an unstable fixed point on the unstable branch. The reduced problem is given bẏ this means that the slow flow on the left branch as well as on the middle branch below the unstable equilibrium points downwards while it points upwards on the right branch as well as between the unstable node and the fold point p F .
Combining this information we now want to construct the candidate orbits in the singular limit. Any orbit starting away from the critical manifold will first follow the fast flow converging to one of the attracting branches of C 0 . The orbits on the right branch follow then the slow flow upwards to p F where they jump with the fast flow to the left branch. There, all orbits follow the slow flow downwards converging to the global equilibrium (0, 0). Orbits starting on the unstable branch of the critical manifold will either converge to p F and jump to the left branch if they start above the unstable node or they will converge downwards towards the saddle if they start below the unstable fixed point.
Finally we show in the following two theorems that the Karma model for ε > 0 has an equivalent behaviour as in the singular limit. To prove this we want to apply geometric singular perturbation theory (see Appendix A.2 for more details). Nevertheless this theory requires differentiability of the system which we loose when E = 1. Since this line crosses the repelling branch of the critical manifold below the unstable node we excluded the heteroclinic connection between the unstable fixed points in Theorem 3.4. This segment will be analyzed separately in Theorem 3.5.
Theorem 3.4. Away from the heteroclinic segment of C 0 between the saddle and the unstable node, candidate orbits found in the singular limit ε = 0 of equations (3.1) with I = 0 can be perturbed to solution curves of the full system with ε > 0.
Sketch of the Proof. For this proof we first need to divide our phase space along the line E = 1 to be able to guarantee smoothness, therefore we will analyse the left and right parts of the critical manifold separately. Furthermore, Fenichel's theorems require smooth vector fields defined on R 2 . In order to satisfy this condition we extend the systems on each side to the entire real plane so that we will be working with either and the unchanged fast equation (3.5) defined in both cases for (E, n) ∈ R 2 . Since all the results of Fenichel's Theory are local around the subset of the critical manifold we are focusing on, these extensions do not change the results. Away from the fold point p F we determined above that the critical manifold is normally hyperbolic so, taking any compact subset of the left branch, we are able to apply Fenichel's first Theorem A.1 to perturb the attracting and repelling branches separately to locally invariant manifolds of the full system. Furthermore, by Fenichel's second and third Theorem A.2 and A.3 the switching between the fast and the slow flow is also preserved for ε > 0. Returning now to our original system (3.1) we can extend the fast fibres over E = 1 using the continuity of the flow. Last it remains to prove that the switching at p F is preserved as well. We know that this point is a generic fold so we can do a coordinate transformation to normal form. Krupa and Szymolyan presented in detail in [30] the analysis of the normal form of a generic fold by performing a geometric blow-up (see also Appendix A.3) which applied to our model concludes the proof.
Theorem 3.5. The heteroclinic segment of the critical manifold in the Karma model (3.1) with I = 0 can be perturbed to a heteroclinic orbit between the equilibria for ε > 0.
Sketch of the Proof. Following the proof of Theorem 3.4 we are able to perturb any compact subset of the heteroclinic segment without the point at E = 1 but we cannot directly guarantee that the left and right subsets connect. To demonstrate the existence of the expected heteroclinic orbit connecting the unstable equilibrium to the saddle we need to look directly at the system with ε > 0.
Since both unstable manifolds of the saddle converge by the analysis above to the origin we are able to define an invariant set delimited by them as shown in Figure 4. Furthermore, from the previous analysis we know that every orbit starting away from the heteroclinic segment will eventually converge to the origin so we follow that there is no periodic orbit and therefore no limit cycle in this set. Now we are able to apply the Poincaré-Bendixson Theorem in the limit t → −∞. Since the origin is unstable in backward time and there are no limit cycles the theorem shows that in fact the now unstable manifold of the saddle needs to converge to the now stable equilibrium proving the existence of the expected heteroclinic orbit for positive ε. : Phase plane (E, n) for ε > 0. In grey we have the invariant set enclosed by the unstable manifolds of the saddle (blue and orange). Furthermore we have its stable manifold for E < 1 (yellow), the nullclines (dotted) and the equilibria of the system (black) In summary, we have shown that although the general structure of the Karma ODE model is similar to the FHN ODE model, there are subtle mathematical differences, particularly in the case of using the standard variants in the literature.

Karma: External current I > 0
Next we focus on the case where the external current I > 0. In the FitzHugh-Nagumo model it is well known that an external current shifts the critical manifold upwards as shown in Figure 5. Without changing anything else, the dynamics switch from a stable resting state to an oscillatory behaviour to a stable depolarized state as the input I increases. To mathematically show this different behaviours we can perform an analogous, yet more complicated, analysis as presented in Section 3.1.1; see also [31]. Figure 5: Phase plane of the FitzHugh-Nagumo system for I = 1 (left) and I = 2 (right). We can see the critical manifold (orange) and the w-nullcline (dashed) as well as the global equilibrium (black) and a prototypical orbit (yellow) oscillating for I = 1 and converging to the stable equilibrium for I = 2.
In contrast to that, adding an external current to the Karma model results in a significant change in the shape of the critical manifold as shown in Figure  6. While we have a regime of I where the critical manifold is "S"-shaped comparable to FHN giving rise to similar relaxation oscillations, when I is big, the manifold flattens out in such a way that the curve is monotonous. In particular this means that the model does not allow any relaxation oscillations or pulses for a high input I. Furthermore, in the Karma model the stable resting state disappears when it collides with the saddle in a fold bifurcation whereas in FHN only the stability of the already unique equilibrium changes. Lastly, the change of stability of the unstable node is for the most part independent of the shape of C 0 . This means that, depending on the model parameters, we can observe bistability as well as a relaxation pulse with a stable depolarized state similar to FHN in addition to the dynamics we have already described. The next theorem formalises all the different dynamic regimes described above.
Theorem 3.6. In the singular limit ε = 0 of the Karma equations (3.1) with I > 0 we have the critical manifold .
Every candidate orbit can be constructed as concatenation of fast segments converging to one of the outer branches of C 0 which are attracting and slow segments on the critical manifold but we need to differentiate multiple parameter regimes for I giving rise to different overall dynamics. The first threshold is given by I = I 2 where the equilibrium with n > 0 changes from unstable to stable as I increases. Furthermore we have • I < I 0 ≈ 0.08718: The equations have three equilibrium points with a stable node and a saddle on the E-axis. If I < I 2 the behaviour is equivalent to the case I = 0. Otherwise, the system is bistable.
• I 0 < I < I 1 = 4 9 : The equations have a unique equilibrium point. If I < I 2 the equilibrium is unstable and the system has a stable relaxation oscillation. Otherwise, the equilibrium is globally stable although there are relaxation pulses.
• I > I 1 (> I 2 ): The middle branch of C 0 disappears such that the critical manifold is attracting everywhere and the unique equilibrium is globally stable.
Sketch of the Proof. Analogously to the previous section, we are now going to study separately the fast and slow subsystems in the singular limit in order to proof Theorem 3.6.
Layer problem: The layer problem is defined by the equation for n fixed. We can see directly that the derivative of the right-hand-side is still given by (3.6). Since by definition any equilibrium is contained in C 0 we can plug in the equality and rewrite that way the Jacobian depending on the external current I instead of n as follows To isolate any non-hyperbolic equilibrium of the fast system we set the J(E; I) = 0 and obtain after simplifying Solving the quadratic equation for E we find 2 curves of non-hyperbolic equilibria given by which connect and disappear for I ≥ I 1 := 4 9 . It is important to check for which values of I the curves E ± (I) are in fact on the critical manifold, more precisely, whether 14) The curve E + (I) satisfies this inequality for all I ∈ [0, I 1 ] but for the curve E − (I) the inequality (3.14) is only satisfied when with I 0 ≈ 0.08718.
Having isolated the non-hyperbolic equilibria we check that, similar to the previous section, when I < I 1 we have a division of the critical manifold into three separate branches where the Jacobian is negative on the outer ones and positive in the middle branch. When I > I 1 the Jacobian stays negative along the whole critical manifold.
Reduced problem: The slow subsystem is still defined by (3.7) but now we have a different definition of C 0 . The most important change lies in the fact that the n-nullcline will, due to continuity, cross the curve E ± (I) in the (E, n)-plane for some I 2 ≤ I 1 dependent on the system parameters n B and M as we increase I. By crossing this curve the global equilibrium of the system (with n > 0) changes its stability and becomes stable. Furthermore we have already seen that the two equilibria at the E-axis collide and disappear for I = I 0 so that for I > I 0 we only have 1 equilibrium of the slow flow.
Remark. Looking at the full system we identify I = I 0 as the bifurcation parameter where the system undergoes a saddle-node bifurcation when the 2 equilibria on the E-axis collide and disappear giving rise to the curve E − (I).
For the corresponding values of I we can check again that the conditions for a generic fold point are satisfied on both curves E ± (I) everywhere except for the point I = I 1 and the singularity at . At the first one the system undergoes a cusp bifurcation where the two fold points annihilate each other. We will come back to this bifurcation later on in more detail. Last, the intersection between the n-nullcline and E ± (I) at I = I 2 satisfies the conditions of a nondegenerate fold but the slow flow is 0. We conclude that at this point we have a fold singularity.
Finally we construct the candidate orbits in the singular limit in the different parameter regimes.
• I < I 0 : If I < I 2 the orbits are equivalent to without incoming current. If I > I 2 then the fast flow will converge to one of the attracting branches of C 0 but while every orbit on the left branch still converges to the origin, contrary to the previous case, all orbits on the right branch will stay on that branch converging to the second stable equilibrium. The slow flow on the repelling branch converges like before to the saddle point.
• I 0 < I < I 1 : First every orbit follows the fast fibres to one of the attracting branches of the critical manifold. If I < I 2 the slow flow leads then to the next fold point where we can again use a fast fibre to jump to the other attracting branch forming a cycle. The flow on the repelling branch will converge away from the equilibrium to the folds following from there the cycle. If I > I 2 and assuming the n-nullcline crosses E + (I) then the flow on the left and middle branch will still converge to the minimum jumping to the right branch. There all orbits converge to the equilibrium.
The case where the n-nullcline intersects E − (I) is equivalent subject to interchange left and right and taking the maximum instead of minimum.
• I 1 < I: The entire critical manifold is attracting so every orbit flows fast to it and then converges to the unique equilibrium.
Remark. To go briefly into the biophysical implications of the above observations we note that all I 0 , I 1 , I 2 are important thresholds affecting differently the behaviour of the cell. Whenever we have a background stimulation I > I 2 , any cell which depolarises over this threshold would not be able to repolarise anymore and will therefore cease to "fire" further signals. On the other hand, a background stimulation I 0 < I < I 1 and I < I 2 results in a self-excitatory system which will "fire" regularly. Finally, when the background stimulation is higher than I 1 the cell will automatically depolarise so that any future signal is blocked.
Theorem 3.7. Whenever E − (I) = 1, candidate orbits found in the singular limit ε = 0 of equations (3.1) with I > 0 away form the bifurcation points I 0 and I 2 can be perturbed to solution curves of the full system with ε > 0.
Sketch of the Proof. Analogously to Theorem 3.4 we find that also for I > 0 away from the intersection between E = 1 and the critical manifold we can perturb every orbit as expected for ε > 0. In the case when E − (I) > 1, in particular when I > I 1 , this point lies in the left branch of C 0 . After continuing the slow manifold obtained for E < 1 over this line we can use the attracting properties of the slow manifold for E > 1 to follow that both manifolds will approach each other. Recalling that the slow manifold is not unique we can directly choose the continuation of the left part to also be our representative slow manifold for E ≥ 1. To finish the proof we need to separate the different parameter regimes when E − (I) < 1. If we first take I < I 2 we have the following cases.
• When I < I 0 the system is equivalent to the case with I = 0 and the proof of Theorem 3.5 can still be applied to derived the heteroclinic orbit between the unstable node and the saddle point.
• When I 0 > I > I 1 we have already derived a stable limit cycle with the unstable fixed point the only orbit not converging to it. In particular we know there are no further periodic orbits inside the limit cycle. This means that, defining an invariant set delimited by the cycle, we can use the Poincaré-Bendixson Theorem to show that the segment of repelling slow manifold with E > 1 will converge to the limit cycle for t → ∞ as well as that the segment with E < 1 will converge to the equilibrium for t → −∞.

Pure ODE models
Finally we look at the system with I > I 2 . By reversing time the repelling branch of the critical manifold becomes attracting and so we can use the same technique applied above and choose the continuation of the left segment of the middle branch as slow manifold. Following the analysis given by Fenichel's theorems and geometric blowup we follow that the middle branch of the slow manifold flows over the fold point diverging in backward time. In the case where I < I 0 this manifold defines a separatrix dividing the phase space into the basins of attraction of the two stable equilibria. In the case where I > I 0 this slow manifold separates the orbits converging directly to the stable equilibria and the orbits which perform a relaxation pulse over the left or right branch of C 0 and one of the fold points before converging.
The limit case with E − (I) = 1 cannot be analysed with the methods used above since the geometric blowup also requires higher regularity. By continuity we would expect that we can still perturb the candidate orbits for ε > 0 but this still has to be proven rigorously. Furthermore, when I = I 0 or I = I 2 the system has folded singularities. It is known that in small neighbourhoods around this points we can find canards and so-called canard explosions. For more details about this solutions see [13,30,31,32].
Remark. All the existence results obtained by Fenichel's Theory require ε to be "small enough". In applications we need to check for every case independently what "small enough" means specifically.
By looking at simulations we see that when I < I 2 the orbits behave as expected even for relatively large ε ≈ 10 −1 . Nevertheless, when the equilibrium changes stability for I = I 2 , the orbits oscillate around the equilibrium instead of converging through a slow manifold as expected from Fenichel's Theory even for very small ε ≈ 10 −4 . This shows that even knowing that there exists an ε for which this theory is applicable, for some values it is not the case. To understand what actually happens at this point with reasonable ε we have to look at the eigenvalues of the fast subsystem as well as of the full system.
Although the Jacobian J of the fast subsystem is strictly smaller than 0 the critical manifold stays very close to non-hyperbolicity and so the absolute value of J is very small. If we calculate the eigenvalues of the full system at the unique equilibrium we have: It holds that J < 0 and the equilibrium is away from E = 1 and E = 4 therefore we see that the parenthesis in the second term of the discriminant is always positive and of order 1 w.r.t. ε so the whole summand is in O(ε). Nevertheless, if J ∈ O(ε) then the first summand is of order ε 2 such that the eigenvalues become complex and our equilibrium is a stable spiral instead of a stable node.
Nevertheless, the equilibrium is still globally stable and every orbit will eventually converge to it.
Extended system (E, n, I). We have seen that an external current has an important impact on the Karma model. Therefore we next investigate an extended 3-dimensional systems with the additional slow equation  This representation allows for proper analysis of the fold curves. Although both models show exactly two curves of folds, their behaviour is clearly very different. Firstly we notice the extreme sensitivity of the Karma model to external currents close to I 0 which is not present for FHN. Furthermore, both fold curves in the FitzHugh-Nagumo model are parallel to each other while in the Karma model they collide and disappear like we had seen above. This collision point I = I 1 , E = 4 3 and the corresponding n given by (3.9) defines a cusp bifurcation in the Karma model. It is clearly a non-hyperbolic point and therefore it cannot be analysed using classical Fenichel theory. Nevertheless, we can do a similar analysis as for a fold point using a coordinate transformation to normal form and geometric blow-up. A detailed analysis of a cusp point using these techniques was presented by Broer et al. in [6].
We have shown above that both models exhibit relatively similar qualitative behaviour when considering I a fixed parameter. Nevertheless, the existence of a cusp singularity presents the possibility for a diverse set of behaviours if we consider for example a slowly changing external current e.g. a slow periodic input. By extending the Karma model considering a change in I we get the possibility of relaxation oscillations with a smooth return. This means that we can have oscillations whereby after a fast jump we are able to return to our starting point following only the slow dynamics. This type of behaviour is not possible in the FitzHugh-Nagumo model even after allowing changes in I.
Remark. The analysis presented in [6] assumes that the bifurcation point is not an equilibrium of the full system. In the Karma model this is in general the case, nevertheless we have the case where I 1 = I 2 when the cusp point is in fact a global equilibrium. This case has (to our knowledge) not been analyzed mathematically yet and would be an interesting future extension to the current analysis.

Travelling waves
As the next step in our analysis we want to consider also the diffusion in the models concentrating on the existence of a travelling pulse in the 1D case. Like before, the existence as well as stability of travelling waves for the FitzHugh-Nagumo equations has been studied extensively, see for example [20,22,25]. We are particularly interested in the construction of pulse solutions performed in [21]. There, Guckenheimer et al. looked at the asymmetric FHN equations with the parameter values γ = 1, a = 1 10 and D = 5. The system is very similar to (1.1) also controlled by a cubic critical manifold in the ODE case. When we add the diffusion term the system exhibits travelling pulse solutions which they proved using a numerical continuation method for the fast fibers in the co-moving frame. In the parameter space (w, c) the authors found a "V"-shaped curve of fast heteroclinic fibers connecting the left and right branches of the critical manifold. When c = 0 the system is Hamiltonian and there is a w = w * such that there is a double heteroclinic orbit. When w is smaller than w * we have a connection from the left branch to the right one while when w is bigger the connection goes in the opposite direction. A concatenation of this fibers combined with the slow flow on the critical manifold can then be perturbed analogously to the previous section for ε > 0, although the technical details become mathematically very involved.
Here we carry out a similar analysis for the Karma model starting by introducing the corresponding co-moving frame z = x + ct such that the equations are now given by We can easily transform the model into a first order system by introducing an additional variable w We now have two additional parameters with respect to the ODE model, namely c and D. The parameter c gives the velocity at which the travelling wave moves. Changing the sign of the parameter c is equivalent to inverting the direction of the wave variable z and substituting w by −w. Therefore, without loss of generality we can restrict our analysis to c > 0.
The second parameter D is the diffusion coefficient. For this parameter there are different scalings often used in the literature. Specifically in the original papers introducing the Karma model the author presents a diffusion coefficient D ∈ O(ε), introducing therefore a third scale to the system (see [28]) while a constant diffusion D ∈ O(1) was used in [29].
Below we focus on the model with D ∈ O(1) and for simplicity only the case without incoming current I = 0. In the following theorems we want to illustrate that the Karma model (2.12) can exhibit a travelling pulse solution with the resting state (0, 0) as start and end state. We sketch the geometric idea of the proof of this result. The model, after transformation to the first order system (3.17), is a (2, 1)-fast-slow system with one-dimensional critical manifold given by Reduced system. The slow flow on C 0 differs from the one in the ODE model only by a factor 1 c so we are simply scaling the flow. In particular, we have the same global equilibria as before embedded into the (E, n)-plane.

Layer problem. The fast subsystem is defined by the equations
The equilibria correspond to the points on the critical manifold for the different values of n. By choosing a different representation we have the fixed point p 0 = (0, 0) and for n M ≤ 1.0415 , 0 The Jacobian at this points is given by with eigenvalues and, when λ ± are real, corresponding eigenvectors v ± = 1 λ ± We can directly check that the equilibria p 0 and p 2 are saddles and p 1 is unstable. In addition we know that p 1 is a node when and a spiral otherwise.
Given the local structure around the critical manifold we want to find heteroclinic connections between p 0 and p 2 to later combine with the slow flow to heteroclinic candidate orbits. (ii) For n = 1 there exists a c min such that for every c ≥ c min the system has a heteroclinic connection from p 2 to p 0 .
Sketch of proof of (i). In order to prove the first statement we are going to follow the strategy in [21]. Our fist step is to compute the stable and unstable manifolds of p 0 and p 2 by taking initial conditions close to the equilibria on their tangent spaces. Next, we define the plane Σ where E = E 2 2 and calculate the intersection points q 0 and q 2 with the previously computed orbits depending on n M and c. The zeros of the function ∆(n M , c) = q 0 (n M , c) − q 2 (n M , c) (3.22) define finally the parameters which give rise to heteroclinic orbits in the fast subsystem. Once we have one such parameter pair, the complete curve in the parameter space can be found because of continuity by slowly changing n M and computing again the zeros of ∆. Figure 8 shows the computed zeros. The left branch of zeros reaching from n = 0 to n = M 15/16 corresponds to the intersection of the unstable manifold of p 0 with the stable manifold of p 2 . The right branch (see close-up) corresponds to the unstable manifold of p 2 intersecting the stable manifold of p 0 . This numerical computation could then be made rigorous, e.g., via employing rigorous numerical techniques, which are already well-established in the context of FHN [1], which concludes the proof of first part of statement (i). For the last part of the statement we observe that in the limit c = 0 the fast subsystem is Hamiltonian with the first integral given by We calculate directly that the energy level at the origin is always 0 and H(p 2 ) = 0 holds if and only if n M = 15 16 . Together with the results illustrated in Figure 8 this strongly indicates that for (n M , c) = (15/16, 0) the system has a double heteroclinic orbit. We can confirm this by computing the energy level H(E, w) = 0 as shown in Figure 9. Before we continue illustrating the geometric ideas behind the proof of the second part of Lemma 3.1, we want to make a remark regarding the previous construction.
Remark. In Figure 8 we can see the curve of heteroclinic orbits for different values of M . In particular we can see the insensitivity of the wave-front velocity with respect to the slow variable n when M 1 which is one of the important advantages mentioned in [28] of the Noble and Karma model over FitzHugh-Nagumo.
Remark. In [9] Deng proved that in the FitzHugh-Nagumo model, under certain conditions, the perturbation of a double heteroclinic orbit in the full system can result in infinitely many front and back wave solutions with an arbitrary number of oscillations. Although his results are not directly applicable in our situation as we would have to adjust the slow variable nullcline to obtain two full system equilibria on the two saddle-type branches, the existence of a double fast subsystem heteroclinic orbit in the Karma model clearly indicates already the possibility of more complex travelling waves than just single pulses.
Sketch of proof of (ii). We have seen in the previous part that the unstable manifold of p 2 and the stable manifold of p 0 connect uniquely for c = c min ≈ 0.707. Nevertheless, for c > c min we find a negatively invariant set enclosed by the E-axis, the stable manifold of p 0 and unstable manifold of p 2 and the vertical segment connecting them at E = 1 as shown in Figure 10. Since we know there are no further equilibria in this set and therefore also no limit cycle we can apply the Poincaré-Bendixson Theorem to obtain that the stable manifold of p 0 converges for t → −∞ to p 2 through the center manifold giving rise to further heteroclinic connections from p 2 to p 0 . We can see the unstable manifold of p 0 (blue) and the stable manifold of p 2 (orange) and the negatively invariant set enclosed by them (grey).
Sketch of Proof of Theorem 3.8 (continued). We can now easily construct a singular candidate orbit combining the slow and fast segments. Starting at the origin as the resting state we can jump to p 2 by a fast fibre where we follow the slow flow upwards. Since we jumped with c ≈ 1.77 we cannot jump until we reach the fold point at n = 1. Using the additional fast fibres we identified above we are able to jump back to p 0 and there follow the slow flow towards the origin.
Theorem 3.9. The homoclinic candidate orbit found in the singular limit ε = 0 of equations (3.16) can be perturbed to a homoclinic solution of the full system with ε > 0.
Idea of proof. The transition from the singular limit to the regular case can be done analogously to Section 3.1. Away from E = 1 where the system is not smooth and the non-hyperbolic fold point (2, 0, 1) we can apply Fenichel's Theory (theorems A.1 -A.3) to obtain the corresponding orbit in the regular case. Again, we can extend the orbits for E = 1 by continuity since we know that we are away from the critical manifold and finally the fold point can be analysed using geometric blow-up as introduced in Appendix A.3.
We recall that in the FitzHugh-Nagumo model a travelling wave will jump to a fast fibre directly from the normally hyperbolic part of the critical manifold.
We have now shown that in contrast to that a pulse solution for Karma model needs the jump segments generated by the fold point through the center manifold. This is a key difference between the two models. It results in a fixed position of the wave back and a slower repolarization than depolarization rate which Karma already identified as important properties for cardiomyocytes (see [28]).

Numerical simulations
In this section we simulate the full PDE systems with a focus on the Karma model. In particular, we want to interpret the numerical simulations in relation to the analysis presented above in order to understand the PDE dynamics [33] we can actually observe. For this we will use the parameter values ε = 10 −2 , D = 1, M = 4, n B = 0.5 and I = 0 except explicitly mentioned otherwise. Figure 11 shows the evolution of the system initialised with a bump function centered at x = 50. For the Karma model, we see from Figure 11 that in fact for a big enough region in x the dynamics converge to a travelling pulse as we have found analytically. Since the simulations converge to a travelling wave given an arbitrary initial profile it (most likely) follows that the travelling pulse is at least locally asymptotically stable and that it does have a substantial basin of attraction. We have not proven the local asymptotic stability analytically here but this would be an interesting point in future work as it is well-known that the FHN PDE has wide parameter ranges, where stable pulses occur and where geometric techniques allow us to prove stability [25,26].
As a comparison, Figure 12 shows a similar simulation for the FitzHugh-Nagumo model (1.1) with ε = 10 −2 , D = 1 and I = 0. At first sight we see that a big difference between Karma and FitzHugh-Nagumo is the hyperpolarization present only in the second model. Although there are heart tissues which show hyperpolarization, if we want to model e.g. ventricular cells the representation in the Karma model is notably more accurate. Furthermore we recall that the repolarization jump of the travelling wave we constructed in the previous section is ignited differently in both models, once on the fold point and once on the hyperbolic part of the manifold. Figure 13 shows that this is the case as well for the limit wave in the full PDE model. As stated before this is the reason for the slower recovery rate in the Karma equations which gives us a key difference between both models. To finish the numerical analysis we want to take a closer look at the effect of other parameters involved in the models and look first at ε. We start with the Karma model. Following the values introduced in [28,29] we have chosen for our simulations ε = 10 −2 as our basis value. In addition, to make sure that the analysis above holds and we have in fact a travelling pulse solution, we need ε to be small enough. Increasing ε shows that already for ε = 0.08 the travelling pulse dynamics seems to break down. Therefore, we will focus on smaller values of ε. By simulating the model with lower values we notice that, as expected, n becomes slower as we decrease ε so that the pulses for E as well as n elongate (see Figure 14). Further we observe in the right panel that the convergence speed towards the travelling pulse is much slower for smaller ε. Nevertheless the wave speed appears to stay unchanged for different values of ε. Since we analytically demonstrated a geometric construction for the existence of the travelling pulses taking the wave speed c as a parameter we would in fact expect changes in c of order ε with c converging to the constant value ≈ 1.77 as ε → 0. It is also intuitively clear from a biological point of view that the wave speed should depend on the properties of the medium, e.g. the diffusion D, but be quite independent of the cells recovery speed. Again we can compare this with the effects of varying ε in the FitzHugh-Nagumo model shown in Figure 15. Overall the effect of varying ε observed in both models is similar. Nevertheless, for ε = 10 −3 we find a change in the wave speed in the FitzHugh-Nagumo model while, as mentioned above, is not visible for the Karma model. We now want to consider the effects of different diffusion coefficients D again starting with the Karma model. As mentioned before, we use as basis value for the diffusion D = 1 for simplicity although the value used in [29] is 2.75. In particular, we would like to make sure that D ∈ O(1) so that the previous analysis applies. Specifically for our model with ε = 10 −2 our simulations lead to assume that D > 0.11 since otherwise the pulse seams to disappear. In Figure 16 we consider three different simulations starting with the same initial conditions for different diffusion coefficients in the range of interest. We see that in this case the wave velocity is as expected strongly affected. An increase in the diffusion rate leads to higher wave velocity. Furthermore we also see that a bigger diffusion coefficient also results in a slightly longer pulse. In the corresponding simulation of the FitzHugh-Nagumo in Figure 17 we see that the effects of different diffusion coefficients on both models are equivalent. Similarly we can look at the control parameters M and n B specific to Karma which are not fixed a priori. Using as a starting point again the values introduced by Karma in [28,29] we follow the range of interest for the parameter M that from modelling point of view varies from M = 4 up to M = 30. Even so, a higher or lower value does not qualitatively change the dynamics of the system. In Figure 18 we see that M has almost no effect on the dynamics of the slow variable n but controls the sharpness of the pulse for E. From biophysical modelling point of view this means that M controls the sensitivity of the voltage E with respect to the gating variable n. On the other hand we know that 0 < n B < 1 and, more precisely, we expect to normally encounter values lying between 0.3 and 0.8. In contrast to the previous case, if we allow n B > 1 then the unstable equilibrium changes stability and the system becomes bistable giving rise to completely different dynamics. Focusing on the range suggested by Karma we find that the parameter n B determines the position of the wave back by controlling the speed of the slow subsystem. The higher n B < 1 the slower is the slow variable and therefore the longer is the depolarisation pulse (see Figure 19). Last we can look at the effect of a small external current I in the Karma model. From the analysis of the ODE model in Section 3.1 we know that for I = I 0 ≈ 0.087 the system undergoes a saddle-node bifurcation so we cannot expect to have equivalent dynamics in the PDE case after crossing this point either. Nevertheless we want to compare the system for I < I 0 since we expect to be able to extend the analysis above in this range. In Figure 20 we see a time shot of the simulations for different values of I. At first sight we see that again the wave speed is changed where the higher the external current the faster the propagation speed of the wave. We can also see that the base line is no longer 0 but slightly higher approaching the fold point as I → I 0 as we would expect. For I = 0.08 we start being able to see that by increasing the base line we also get a weak hyperpolarization after the main pulse which we also would expect analytically due to the shape of the critical manifold. Furthermore we can shift the waves such that the wave fronts coincide and see that the wave profile is also affected by the external current (see Figure  21). Although the effect is not as noticeable as the different wave speed we see that in addition to the higher base line we also have slightly longer pulses for higher incoming current.

Discussion
We presented a systematic analysis and comparison of a polynomial version of the Karma model [28,29] with the FHN model [17,18,19] motivated by applications to modelling excitable behaviour in cardiomyocytes with regard to individual cells as well as cell populations. We started by considering their pure ODE versions. In this setting we noticed that Karma as well as FitzHugh-Nagumo present similar behaviours showing in both cases exactly three parameter regimes for the input current. When I is sufficiently small the dynamics converge to a stable resting state while in the middle range of I both models oscillate following a globally stable limit cycle. Finally, when I is high enough any orbit converges to a stable equilibrium corresponding to a depolarized state. Nevertheless, although both systems are qualitatively similar, there are also some likey important differences when applying them to model cardiomyocytes. First, in the Karma model the re-polarisation is much slower than the depolarisation because of the sharpness of the critical manifold while in the FHN model both processes are of the same order. Also, for a high external input I the dynamics of the FHN model are still controlled by an "S"-shaped critical manifold, in other words, depending on the initial conditions it is possible to undergo a depolarisation and re-polarisation before converging to the stable state. In contrast to that the Karma model does not allow any oscillation other than small fluctuations very close to the equilibrium given that in a reasonable regime for ε the fixed point is a stable spiral. Yet, the biggest difference we see occurs when considering I as a dynamic variable instead of a parameter. In the extended phase space we observe the high sensitivity to changes in I when the system is oscillatory and most importantly the cusp singularity that arises when the two folds collide. Because of this differences it would be interesting for future work to look at the models with non-constant external current I.
Next, we considered the spatially extended versions of the models and focused on travelling wave solutions in 1D without external current. This is motivated by our interest in modelling propagation of activity in populations of cardiomyocytes similar to [7,8,44]. We start by analysing the 1D PDE in the singular limit ε = 0 in order to study the existence of travelling wave solutions. Here, using similar techniques as used for FHN in [21] in addition to singular perturbation theory, we have demonstrated the existence of travelling pulses originating and converging to a fixed resting state. The first difference we have found comparing the Karma model to FHN is the insensitivity of the wave speed to different values of the slow variable. Furthermore, in contrast to FitzHugh-Nagumo, the wave back in the Karma model starts at the fold point for large parameter ranges resulting as in the ODE system in much slower repolarisation than the previous depolarisation. All the analysis in this section has been restricted to I = 0, therefore, as a future continuation of the work, it would be interesting to study if it is possible to extend the existence of travelling waves for I > 0 and especially in the range where the ODE is oscillatory.
Finally, we performed numerical simulations of the 1D PDE Karma model varying model parameters. As we would expect, the propagation velocity does not depend on the parameters controlling the reactivity of the cells but only on the parameters defining the medium, namely the diffusion coefficient D and the background current I. On the other hand, while a change in D or I also affects the shape of the pulse we have observe that the main control over the shape is given by the reaction parameters ε, M and n B . Since all these are based only on observations of different simulations it would be another interesting avenue for future work to perform an even deeper analysis of the effect of the parameters on the travelling wave solutions.
A.1 Definitions and notation with its three key theorems and finally in Section A.3 the blow-up technique is briefly introduced.

A.1 Definitions and notation
First, we introduce some important notions for the work with dynamical systems in general defined by the ODE for x : I → R n with I ⊆ R an interval and some function f ∈ C r (R n , R n ) (r ≥ 1). We denote the flow induced by the differential equation (A.1) as φ t (·).
Definition. For a set S ⊂ R n and a manifold M ⊂ R n we have: • S is called positively invariant if for all p ∈ S it holds that φ t (p) ∈ S for all t ≥ 0.
• S is called negatively invariant if for all p ∈ S it holds that φ −t (p) ∈ S for all t ≥ 0.
• M is called locally invariant under φ t (·) if for each p ∈ int(M ) there exists a time interval I p = (t 1 , t 2 ) such that 0 ∈ I p and φ t (p) ∈ M for all t ∈ I p . In other words, the flow can only leave the manifold through its boundary.
We now continue with the systems showing multiple time scales. This is the case when some of the variables are much slower than others so we can identify two separate subsystems that move at different time scales.
The setting we will be working with is a (m, n)-fast-slow system, this means we have an m-dimensional fast subsystem combined with n further variables moving at a slower time scale. The complete system is given by the (m + n)dimensional system of differential equations where f ∈ C r (R m+n+1 , R m ) and g ∈ C r (R m+n+1 , R n ). For the equation (A.2) the components of x = x(t) ∈ R m are the fast variables and those of y = y(t) ∈ R n are the slow variables of the system. The time scale separation is controlled by a small parameter ε > 0, which provides the ratio between the slow time scale t to the fast time scale τ := t/ε. Observe that (A.2) describes the system with respect to t and therefore we call it the slow system. If we write the system on the fast time scale τ = t/ε, we obtain the fast system defined as When looking at fast-slow systems we are mostly interested in the case when ε is very small and that brings up the question what happens at the limit ε → 0 for a given system, the so-called singular limit. Although the slow and fast system (A.2) and (A.3) are equivalent since we only need a reparametrization of time to transform them into each other, their corresponding singular limit is not, so we have to differentiate between the reduced system or slow subsystem which we obtain taking the singular limit of the slow system (A.2) and the layer problem or fast subsystem as the limit for the fast system (A.3). The reduced system (A.4) consists of an algebraic constraint and the differential equation for the slow variables defining the so-called reduced or slow flow. In comparison, the layer problem is given by a differential equation for the fast variables where y is a fixed parameter and the associated flow is called fast flow. The name arises as every value of y defines an independent "layer" of the system.
Definition. The algebraic constraint in (A.4) defines the critical manifold The points contained in the critical manifold correspond exactly to the equilibrium points of the fast subsystem.
Remark. The equations (A.4) define the slow flow to be naturally restricted to the manifold C 0 .
In this setting we are able to decouple the fast and the slow dynamics in the system analysing them separately. Nevertheless, to obtain a global picture of the full system we need to combine the fast and the slow trajectories and we get the following definition.

Definition.
A candidate orbit is the image of a homeomorphism γ : (a, b) → R m+n with a < b and a partition a = t 0 < t 1 < · · · < t k = b for some k ∈ N + such that • the image γ((t i−1 , t i )), i ∈ {1, . . . , k} of each subinterval is a trajectory of either the fast or the slow subsystem • the image γ ((a, b)) has an orientation that is consistent with the orientation of each trajectory γ((t i−1 , t i )), i ∈ {1, . . . , k}

A.2 Fenichel's Theory
The following statements were first introduced by Fenichel in his paper "Persistence and smoothness of invariant manifolds for flows" 1971 ( [14]) and then applied to fast-slow systems 1979 in "Geometric singular perturbation theory for ordinary differential equations" ( [15]). Fenichel's work consists of three main theorems posed in a very general setting. Since we will not need this generality, we will only present the important results already applied to fast-slow systems. We will not prove any of the statements in this section, for the proofs see [14,15,32,43].
To be able to understand the next theorem we first need some more properties of the critical manifold.
Definition. A subset S ⊂ C 0 is called normally hyperbolic if the Jacobian with respect to the fast variables D x f (x * , y * , 0) ∈ R m×m has no eigenvalue with zero real part for all (x * , y * , 0) ∈ S.
Remark. The definition shows that S ⊂ C 0 is normally hyperbolic if and only if for every (x * , y * , 0) ∈ S it holds that x * is a hyperbolic equilibrium of the fast subsystem for y = y * i.e. x * is a hyperbolic equilibrium of x = f (x, y * , 0).
Using this correspondence between fixed points of the layer problem and points of the critical manifold we can now analyse the stability properties of those equilibria and define the analogous concept for C 0 .
Definition. Let S ⊂ C 0 be a normally hyperbolic set.
• S is called attracting if for all (x * , y * , 0) ∈ S every eigenvalue of D x f (x * , y * , 0) has negative real part i.e. for all (x * , y * , 0) ∈ S the corresponding equilibrium x * of the fast subsystem is stable for y = y * .
• Similarly, S is called repelling if all eigenvalues have positive real part i.e. the fixed points are unstable.
• If S is neither attracting nor repelling, it is called of saddle type.
Finally, to measure the distance of the perturbed manifolds we will use the following metric.
Definition. The Hausdorff distance d H between to nonempty sets V, W ⊂ R k is defined by where dist(p, M ) := inf q∈M ||p − q|| gives us the distance from a point p ∈ R k to the set M ⊂ R k . In other words, the Hausdorff distance d H (V, W ) defines the maximal distance between a random point in one set to the other set.
Theorem A.1 (Fenichel's first Theorem, fast-slow version). Let S 0 be a compact normally hyperbolic submanifold of the critical manifold C 0 of (A.2) and f ∈ C r (R m+n+1 , R m ), g ∈ C r (R m+n+1 , R n ) for 1 ≤ r < ∞. Then for ε > 0 sufficiently small it holds that (F1) There exists a locally invariant manifold S ε diffeomorphic to S 0 , The flow on S ε converges to the slow flow as ε → 0, Definition. The manifold S ε obtained as conclusion of Theorem A.1 is called a slow manifold.
Remark. S ε is usually not unique. Nevertheless, in regions lying at a fixed distance from ∂S ε , all manifolds satisfying (F1)-(F4) lie at a Hausdorff distance O(e −K/ε ) from each other for some positive K ∈ O(1). For this reason, a representative of the manifolds is often called "the" slow manifold since, in most cases, it is arbitrary which representative to choose.
By Theorem A.1 we know that, starting with a fast-slow system in the singular limit, if we perturb the equations by taking ε > 0 sufficiently small the structure and behaviour of the critical manifold do not disappear. Instead, any compact subset of the manifold perturbs continuously to a slow manifold S ε . The perturbation does not only preserve the topological structure of the critical manifold but the flow on the slow manifold is also ε-close to the original slow flow. In particular, W s,u loc (S ε ) exist and furthermore S ε is normally hyperbolic and has the same stability properties with respect to the fast variables as S 0 (attracting, repelling or of saddle type).
Although the stable and unstable manifolds can only be defined for an equilibrium , with the help of Fenichel's second Theorem we are able to generalize this notion to fast-slow systems with ε small enough. The stable and unstable manifolds W s,u loc (S 0 ) that result from the union of those of the individual fixed points of the fast subsystem are not lost by the perturbation but instead we find that the topological as well as the analytical properties in new manifolds W s,u loc (S ε ) remain similar to those of the original ones.
Theorem A.3 (Fenichel's third Theorem). We start again with the same setting as in Theorem A.1. Then there exists a manifold F u (p) for each p ∈ S 0 such that is tangent to N u p at p with N u p the unstable component of the normal direction to S 0 (fast direction), (e) There exist constants C u , λ u > 0 such that if q ∈ F u (p), then for every t ≥ 0, (f ) F u (p) is C r with respect to the base point p.
The same conclusions (a)-(f ) with the obvious modifications hold for the family of manifolds F s (p), e.g. replace −t by t in (c) so that φ t (F s (p)) ⊆ F s (φ t (p)). Furthermore, the foliation persists for ε > 0 with all properties mentioned above and diffeomorphic to the foliation in the singular limit.
Definition. The manifolds F s,u (p) are called the stable/unstable fibres through p.
The families F s,u build decompositions of the stable/unstable manifolds by submanifolds characterised by initial conditions approaching each other at the fastest rate in forward/backward time. They are therefore called asymptotic rate foliation. Since this foliation persists under perturbations we can conclude that the asymptotic behaviour on the stable and unstable manifolds stays unchanged for positive ε.
These theorems are often referred to as geometric singular perturbation theory or GSPT as we do in this paper. Nevertheless, this term can also describe a wider compilation of geometric techniques for the analysis of singularly perturbed systems and consequently it is important to clarify if it refers to a specific result or to a hole branch of methods.

A.3 Geometric blow-up
While Fenichel's Theory gives us the tools to analyse the critical manifold and its stable and unstable manifold whenever it is normally hyperbolic, the method of geometric blow-up enables us to investigate the points where the hyperbolicity is lost. The most common example are fold points as we see multiple times in this paper but also i.e. other types of bifurcation points like the cusp bifurcation appearing in Section 3.1.3. This section aims at giving a basic understanding of the idea behind blow-up in general and applied to fold points in particular. It is not thought as a deep or complete study of this matter, see [32] for a more extensive introduction or the original papers [6,11,12,30] for detailed statements and proofs.
The method of geometric blow-up was first introduced in [11,12] independently of systems with multiple time scale dynamics. The most simple case is given for a planar vector fieldż = F (z) ∈ R 2 as follows Definition. Let S 1 be the unit circle with the corresponding polar coordinate transformation φ : S 1 × I → R 2 , φ(θ, r) = (r cos θ, r sin θ) for some (possibly infinite) interval I ⊂ R with 0 ∈ I and θ ∈ [0, 2φ) as parametrization of S 1 . Then the polar blow-upF of a vector field F ∈ C ∞ (R 2 ) with F (0) = 0 is defined bŷ F (θ, r) := (Dφ −1 (φ,r) • F • φ)(θ, r) (A.6) for r = 0 and by the continuous extension of (A.6) to r = 0.
Remark. The name polar blow-up clearly arises from the polar coordinate transformation used for r > 0. Since the blow-up method can also be performed with a different coordinate transformation (see directional blow-up introduced in [32]), it serves as specification of the coordinate transformation that has been applied. Nevertheless, one can proof that both methods of blowup, polar and directional, are equivalent up to a coordinate transformation.
Definition. Let F be a C ∞ vector field as above. We define the (rescaled) polar blow-up asF := 1 r kF with k such that the derivatives of F satisfy D k F = 0 and D k+1 F = 0. It is important to notice that this scaling does not change the qualitative structure of the orbits when r > 0.
To illustrate the idea behind the blow-up method we present a simple planar example that can be desingularized using a polar blow-up transformation.
Example. We consider the systeṁ It holds that F (0, 0) = 0 as well as DF (0,0) = 2z 1 − 2z 2 −2z 1 −2z 2 2z 2 − 2z 1 (0,0) = 0 0 0 0 so we see that the origin is a non-hyperbolic equilibrium. We now want to apply a blow-up transformation to construct a vector field with hyperbolic equilibria instead. The polar blow-up vector field is given bŷ F (θ, r) = 3r cos θ sin θ(sin θ − cos θ) 1 4 r 2 (cos θ + 3 cos(3θ) + sin θ − 3 sin(3θ)) so we can define the rescaled polar blow-up byF (θ, r) = 1 rF (θ, r). In this new vector field we find 6 hyperbolic equilibria on the circle S 1 × r = 0 so we are now able to use the linear structure to describe the dynamics for r small. Since the vector fields outside of r = 0 are conjugated can use this information to describe qualitatively the structure of the orbits close to the origin in the original vector field. We see that we have transformed 1 non-hyperbolic equilibrium into a ring with 6 hyperbolic ones.
We can generalize the definition above by allowing vector fields F ∈ C ∞ (R n ) in higher dimensions with the condition F (0) = 0 unchanged.
A further generalisation can be defined using instead of the classical polar coordinate transformation a generalized one as follows.
Remark. It is clear, that the definition of the rescaled blow-up can be applied not only for two-dimensional vector fields but also to higher dimensional blow-ups as well as weighted polar blow-ups. In this last case we call it the (rescaled) weighted polar blow-up.
Finally we want to briefly present the results for the analysis of a fold point and first need to introduce some notation. The general setting of a fold point (w.l.o.g. located at the origin) is a (1, 1)-fast-slow system as (A.3) satisfying the conditions f (0, 0, 0) = 0 , f x (0, 0, 0) = 0 and w.l.o.g.
We divide the critical manifold into the attracting and repelling branches S a 0 and S r 0 respectively. Analogously we define S a ε and S r ε as the attracting and repelling branches on the slow manifold, see Section A.2 for more details. Furthermore, for some ρ > 0 small we define the sections ∆ in := {(x, ρ 2 ) : x ∈ J in } and ∆ out := {(ρ, y) : y ∈ J out } where J in and J out are intervals such that the transition map Π : ∆ in → ∆ out is well define as shown in Figure 23. Now we can formulate the results of the flow near a fold point. The following theorem proofs that in fact the switch between the slow and fast flow at a fold point is conserved for ε > 0 as shown in Figure 23. As mentioned in the introduction of this section we will not prove this theorem but only give the key ideas behind it, for a full proof see [30,32].
Theorem A.4. There exists ε 0 > 0 such that for all ε ∈ (0, ε 0 ] the following holds: (1.) The manifold S a ε passes through ∆ out at (ρ, H(ε)) with H(ε) ∈ O(ε) (2.) ∆ in is mapped by Π to an interval of size O(e −C/ε ) for some C > 0 (3.) The function H(ε) has the asymptotic expansion H(ε) = c 1 ε 2/3 + c 2 ln ε + c 3 ε + O(ε 4/3 ln ε) as ε → 0 Key steps of the proof. The idea to proof Theorem A.4 is to extend the model to the 3-dimensional system with ε a dynamic variable with derivative 0 and use that the origin is non-hyperbolic. Therefore, in this setting we can in fact apply a blow-up to desingularize the fold point. The correct transformation to achieve this is a weighted polar blow-up with the generalized polar coordinate transformation ϕ(x,ȳ,ε, r) := (rx, r 2ȳ , r 3ε ) = (x, y, ε) Once we have the rescaled vector fieldF it remains to find charts to represent the blow-up in local coordinates. Then we can analyze the dynamics and finally translate the results to the original "blow-down" vector field.