Neural Excitability and Singular Bifurcations

We discuss the notion of excitability in 2D slow/fast neural models from a geometric singular perturbation theory point of view. We focus on the inherent singular nature of slow/fast neural models and define excitability via singular bifurcations. In particular, we show that type I excitability is associated with a novel singular Bogdanov–Takens/SNIC bifurcation while type II excitability is associated with a singular Andronov–Hopf bifurcation. In both cases, canards play an important role in the understanding of the unfolding of these singular bifurcation structures. We also explain the transition between the two excitability types and highlight all bifurcations involved, thus providing a complete analysis of excitability based on geometric singular perturbation theory.

Bifurcation diagrams of the canonical model (1) together with 'frequency-current' (f-I) plots: (Type I) c = 0.005; SNIC bifurcation near I = I bif = 0 where the frequency approaches zero; (Type II) c = 4; supercritical singular Andronov-Hopf bifurcation near I = I bif = 0; the subsequent canard explosion is clearly visible; note the small frequency band for the relaxation oscillation branch A first answer to the question of the neuron's computational properties was given by Hodgkin [1] in the 1940s, who identified three basic types (classes) of excitable axons distinguished by their different responses to injected steps of currents of various amplitudes. Type I (class I) axons are able to integrate the input strength of an injected current step, i.e. the corresponding frequency-current (f-I) curve is continuous (see Fig. 1). Type II (class II) axons have a discontinuous f-I curve because of their inability to maintain spiking below a certain frequency. The frequency band of a type II neuron is very limited and, hence the frequency is relatively insensitive to the strength of the injected current. It appears that type II neurons resonate with a preferred frequency input. Type III (class III) axons will only fire a single or a few action potentials at the onset of the injected current step, but are not able to fire repetitive action potentials like type I and type II neurons. Type III neurons are able to differentiate, i.e. they are able to encode the occurrence of a 'change' in the stimulus. Such phasic firing (versus tonic or repetitive firing) identifies these type III neurons as slope detectors [2]. Obviously, the f-I curve is not defined for type III neurons.
Starting in the 1980s, Rinzel and Ermentrout [3,4] pioneered a mathematical framework based on bifurcation theory that distinguishes type I and type II neural models. Recall, type I and type II neurons are able to fire trains of action potentials (tonic firing) if depolarised sufficiently strong which distinguishes them from type III neurons. This distinction points to a bifurcation in type I and type II neurons where the cell changes from an excitable to an oscillatory state. The main bifurcation parameter is given by I , the magnitude of the current step protocol, and leads to the following classical definition of excitability via bifurcation analysis under the variation of the applied current I [4]: • Type I: The stable equilibrium (resting state) disappears via a saddle-node on invariant circle (SNIC) bifurcation.
Note that there is no bifurcation occurring in type III models, i.e. the equilibrium (resting state) remains stable. For a possible biophysical interpretation of these three distinct excitability types we refer to, e.g., [5,6]. A canonical model that captures the main features of neural excitability is given by with a sufficiently smooth (here C 1 ) function with threshold parameter v th , and parameters 0 < ε 1, d, e > 0 such that d 3 < e < d holds 1 and c ∈ [0, c 1 ]. The parameter I ∈ [I 0 , I 1 ] ⊂ R is the main bifurcation parameter and is often associated with the injected current in a real neuron. In this canonical model, we are able to identify all three excitability types. Figure 1 shows bifurcation diagrams and frequency-current plots for type I and type II excitability where the parameters v th = 0.15, ε = 10 −2 , d = 2, e = 1.5 are fixed and (I, c) are varied.
Remark 1 In the canonical model (1), the nonlinear nature of G(v) is essential to guarantee relaxation type behaviour as observed in type I neurons which distinguishes this polynomial model from the classic FitzHugh-Nagumo model [7,8] which can only produce type II (and type III) behaviour. In more biophysically inspired twodimensional systems, the Morris-Lecar model [9] presents a prime example which is able to produce all three excitability types; see also [4][5][6]10].

Slow-Fast Excitable Systems
An important feature of most neural systems is that they evolve on multiple time scales; see, e.g., [11]. It is the interplay of the dynamics on different temporal scales that creates complicated rhythms. Multiple (or slow-fast) time-scale problems are usually modelled by singularly perturbed systems such as (1) where the time-scale separation of the 'fast' variable v (voltage) and the 'slow' variable w (recovery variable) is explicitly identified through the singular perturbation parameter ε 1. The interest in such slow-fast systems goes towards the presence of so-called relaxation oscillations. Along an orbit of relaxation oscillation type, parts where the velocity of the phase variables is small (the slow parts) are alternated with high velocity peaks on short time intervals (the fast parts). During the slow parts, the phase state is O(ε)close to the critical set f (w, v, I ) = 0 (because then (w , v ) = O(ε)), whereas during the fast part the phase state is at an O(1)-distance from this critical set. Tonic firing as observed in (1) is exactly of relaxation type.
We want to emphasise this inherent slow-fast time-scale structure found in many neuronal models and use geometric singular perturbation theory (GSPT) [12,13] as a mathematical framework. In this approach, we focus on a class of two-dimensional singularly perturbed models given by w = εg (w, v, ε, λ), where I ∈ [I 0 , I 1 ] ⊂ R is an external (constant) drive of the excitable system, the prime denotes the (fast) time derivative d/dt and ε 1 is a small positive parameter encoding the time-scale separation between the slow and fast variables. The parameter λ is considered in some interval [λ 0 , λ 1 ] and will serve, together with I , as a bifurcation parameter. The functions f and g are assumed to be sufficiently smooth. We stress that it is not important for f to be independent of λ and for g to be independent of I , though it does simplify the presentation. 2 By switching from the fast time scale t to the slow time scale τ = εt, system (3) where the overdot denotes the (slow) time derivative d/dτ . As ε → 0, the trajectories of (3) converge during fast segments to solutions of the one-dimensional layer (or fast) problem while during slow segments, trajectories of (4) converge to solutions oḟ which is a one-dimensional differential-algebraic problem called the reduced (or slow) problem. Note that the critical set is the set of equilibria of (5). In general, this set S defines a differentiable manifold referred to as the critical manifold, and it forms the phase space of the reduced problem (6). GSPT [12][13][14][15] uses these lower one-dimensional sub-systems (5) and (6) to predict the dynamics of the full two-dimensional system (3) or (4) for ε > 0. In Sect. 2, we provide the general setup for a class of two-dimensional systems in the context of GSPT that covers all three excitability types.
While this GSPT approach to explain relaxation type behaviour in neural systems is well known to the mathematical and computational neuroscience community, see e.g. [5,11], it is not often or consequently used to explain the underlying bifurcation structure in such singularly perturbed systems. This is the focus of Sects. 3-6 and we study singular bifurcations and their unfoldings in a normal form introduced in Sect. 3.1.
A closer look at the bifurcation diagram of type II excitability in Fig. 1 reveals that shortly after the Andronov-Hopf bifurcation the amplitude of the bifurcating limit cycles explodes dramatically under a very tiny (an exponentially small) parameter change. This is known as a canard explosion [14,15] and indicates that the singular perturbation nature of the neural model is also reflected in the bifurcation structure. Note also a similar dramatic change in frequency near this singular Andronov-Hopf bifurcation. We will review this (well-known) bifurcation phenomenon in Sect. 4.
Additional bifurcation structure is necessary to explain the transition from type II to type I excitability. This is covered in Sect. 5 where an incomplete canard explosion is identified. This lesser-known phenomenon refers to a premature termination of a canard explosion in a homoclinic bifurcation.
Similarly to the singular Andronov-Hopf bifurcation, one has to expect that the SNIC bifurcation associated with type I excitability shown in Fig. 1 must have a singular nature. We identify a singular Bogdanov-Takens/SNIC bifurcation point as the organising centre for type I excitability. Unfolding this type I singular bifurcation structure is the main focus of Sect. 6 which is based on the blow-up method, a desingularisation technique for nilpotent singularities that has been successfully implemented in geometric singular perturbation problems with loss of normal hyperbolicity [14][15][16].
Finally, we summarise our results in Sect. 7 and discuss its implications for possible numerical observations in slow-fast neural models.

The Setup for Slow-Fast Excitable Systems
We start with introducing basic assumptions on the singularly perturbed system (3), respectively, (4).

Assumption 1
For each I ∈ [I 0 , I 1 ], the critical manifold S is cubic shaped and given as a graph {w = ϕ I (v)}, i.e.
with attracting outer branches S ± a , repelling middle branch S r , and folds F ± . We also assume that the vertical fibres containing the two local folds F − , respectively, F + intersect the critical manifold one more time outside the fold points at points p, respectively, q; see Fig. 2. The two folds are assumed to be regular extremes of the graph w = ϕ I (v).
The representation of S as a graph {w = ϕ I (v)} necessarily implies that ∂f/∂v = 0 for all (w, v) ∈ S. The stability condition on the branches S ± a and S r identifies them . This refers to the stability property of S as a set of equilibria of the layer problem (5) and is expressed by The two local extremes of S denoted by (w ± , v ± ) = (ϕ I (v ± ), v ± ) are called fold points F ± ; they correspond to saddle-node bifurcation points in the layer problem (5) and we have Figure 2 can be viewed as the corresponding bifurcation diagram of the layer problem (5) with w as the main bifurcation parameter. If we want to impose the geometry shown in Fig. 2, then we need to-and will-assume that Assumption 2 For each λ ∈ [λ 0 , λ 1 ], the following holds along the w-nullcline g(w, v, 0, λ) = 0 of system (4): This assumption implies that the w-nullcline is a graph {w = ψ λ (v)}, and ψ λ is a monotonically increasing function. While mathematically not necessary, this reflects the property of a typical neural model where g = 0 is given as a graph of a sigmoidal function over v.  • there can be no more than one equilibrium on the lower attracting branch S − a ; • the existence of exactly two equilibria indicates a saddle-node bifurcation either on S r or at F − . Remark 2 Assumption 3 sets the scene for the transition from an excitable to an oscillatory state. It excludes the possibility of a transition from an oscillatory state to a(nother) steady state (known as depolarisation block in neuroscience), by restricting the parameter space from above. This is not necessary but allows us to focus on the onset of oscillations, not the termination.
Next we look at the reduced problem (6). The corresponding one-dimensional dynamics on the critical manifold S projected onto its base coordinate v is given by Generically, system (11) defines the reduced dynamics along the hyperbolic branches of the critical manifold, outside the fold points F ± where (11) is singular. Assumption 3 implies that g(ϕ I (v), v, 0, λ) = 0 on the upper branch S + a .
The first condition in (12) is equivalent to imposing the requirement that ϕ I (v) = 0, ϕ I (v) < 0; in other words, the fold F + is a regular local maximum of ϕ I . 3 Together with the second condition, this determines the direction of the reduced flow near F + and qualifies this fold point as a jump point: all orbits that come in along the upper branch S + a jump off the fold F + and follow the fast fibre towards q; see Fig. 2 (right). Note that no equilibrium on S r can approach the upper fold F + while varying parameters (see also Remark 2).
On the other hand, an equilibrium from S − a may cross or bifurcate at the lower fold F − , which is necessary to observe a bifurcation from an excitable to an oscillatory state. Assumptions 3 and 4 imply that the reduced flow on S − a is either towards an equilibrium on S − a or towards the lower fold F − (see Fig. 2 (right) and Fig. 4), another essential feature for an excitable/oscillatory system. These basic assumptions hold for many two-dimensional neuronal models including the canonical system (1). This model has a cubic-shaped critical manifold S (Assumption 1). The function G(v) is monotonically increasing (Assumption 2). It is the nonlinear nature of G imposed by Assumption 3 that allows us to explore the cases of different numbers of equilibria (restricted to S − a and S r only). If 0 < I 1 < I F + and v th < v + , where I F + is defined implicitly by ϕ I F + (v + ) = G(v + ), then Assumptions 3 and 4 are fulfilled and F + is a regular jump point for all I ∈ [I 0 , I 1 ] and c ∈ [0, c 1 ].

Singular (or Slow-Fast) Bifurcations
Our main task is to identify bifurcation structures in the singularly perturbed system (4) based on singular observations, i.e. for ε = 0. For system (1) When I < I bif , there is a unique singularity on the stable branch S − a with possible extra singularities on the middle branch while for I > I bif there is no singularity on S − a . Hence, the transition from an excitable to an oscillatory state happens near I ∼ I bif and we will distinguish different excitability types by different behaviour near the lower fold F − .

Remark 3
System (1) models type III excitability if I bif lies outside of the interval [I 0 , I 1 ], i.e. I bif > I 1 . In general, type III excitability can be characterised by the property that I bif lies outside of the interval [I 0 , I 1 ].

Normal Form Near the Singular Fold Point
In the canonical model (1) assuming it is not of type III, the fold F − is a singular fold point at I = I bif defined in (13). We can make the simple coordinate transformations w = x/d + I and v = y/d to write the system in the form where a = dI . Note that, when v > v th , an extra term ε e d (y − dv th ) 2 appears in the x equation.
While in general one has to do a bit more work, the normal form for system (3) near the singular fold F − shows similarity with the normal form (14) of the canonical model. The local shape of the vector field (3) near the fold F − is described in the following proposition.
where σ = ±1, σ c ≥ 0 and β = 0. The coefficients a, c and β can be computed explicitly in terms of (I, λ, ε). Proof Starting with system (3), and working under Assumption 1 and (10), we know from the implicit function theorem that the v-nullcline f (w, v, 0, I ) = 0 is a graph w = ϕ I (v). By the same argument, we can factor f (w, v, ε, I ) as (φ I,ε (v) − w) ·f for some strictly positive functionf (w, v, ε, I ), and withφ I,ε (v) = ϕ I (v) + O(ε). By dividing the vector field byf (i.e. by rescaling time), we obtain the (topologically) equivalent vector field 4 The reader may verify that Assumptions 1-3 remain unchanged after this manipulation. 5 Since ϕ I (v) > 0 and henceφ I,ε (v) > 0, we can solveφ I,ε (v) = 0 using the implicit function theorem to find an I -dependent point at v = v * (I, ε), ε-close to the fold of ϕ I (v). 6 Replacing v by v − v * and w by w −φ I,ε (v * ) allows one to assume in the remainder that the fold ofφ I,ε (v) is located at the origin. It suffices to see that where the coefficients A k and B 0 depend on (I, λ). Let us now write (w, v) = (Zx + εP y + ε 2 Q, v = Zy + εR) for well-chosen (P , Q, Z, R). Aided by a computer algebra program it is easy to verify that one can choose (P , Q, Z) in terms of R to make sure y = y 2 − x + O(y 3 ). Once we have obtained this, it is not so hard to see that one can choose R in such a way that we arrive at In other words, the leading-order coefficient with The first three terms can be put together with the primary part C 0 + C 1 y + C 2 x, by allowing C i to be ε-dependent. It shows that the remainder term is O(x 2 , y 3 , xy, εy 2 ). A linear rescaling of (x, y, t) allows one to further reduce to the case C 2 = ±1. The constraint C 0 C 2 ≥ 0 follows from Assumptions 2-3.

Remark 4
For system (14), the coefficients in the normal form are given by c ≥ 0 and at least when v th > 0. When v th < 0, the parameters a and c are shifted by extra terms (which we will not give details of here) that depend on v th . The extra term in the coefficient a will make a = 0 at I = I bif = ev 2 th . The coefficient c in the normal form in that case depends regularly on v th and can be used to locally distinguish the different scenarios depicted in Fig. 3. Hence, it makes sense to consider the coefficients (a, c) in the normal form as the main bifurcation parameters.
In system (15) for a < 0, there exists a stable equilibrium on the lower stable branch S − a while for a > 0, there is no stable equilibrium on the lower stable branch S − a and the fold F − is a regular jump point. Thus a transition in the dynamics must happen for a ∼ 0 near the lower fold F − and, depending on c, we classify the singular contact point F − as follows: • For a = 0 and c > 0, the fold F − is a singular (or slow-fast) Andronov-Hopf point; see Sect. 4. • For a = c = 0, the fold F − is a singular (or slow-fast) Bogdanov-Takens point; see Sect. 6.
These two local singular bifurcation points are associated with the two different neural excitability types I and II. Note that a Bogdanov-Takens bifurcation is a codimension-2 bifurcation that includes a codimension-1 Andronov-Hopf bifurcation in its unfolding. So, we also expect to find a connection between these two bifurcations as c tends to zero.
Remark 5 Although type III neurons are not associated with any bifurcation for fixed current input, these slope detectors play an important role in identifying dynamic changes and producing transient responses. We refer to [2] for details and [10] where type III neurons and excitability are discussed in the context of GSPT. From Fig. 3, we can deduce the presence of another (local) codimension-2 bifurcation, a cusp bifurcation [17] where two codimension-1 saddle-node bifurcations merge. 7 Its approximate location can be computed easily in the normal form (15), under the condition that we discard the higher order terms O(x 2 , xy, y 3  . This indicates that for fixed 0 < c < c cusp we observe three equilibrium states, two of which are located on S r (see Assumption 3). This has interesting consequences on the global bifurcation structure of our problem; see Sect. 5. As can be seen in Fig. 4 for c = 0 and I = I bif , the layer problem of a type I neuron has a saddle-node bifurcation of equilibria at the lower fold F − . This allows for the construction of a singular homoclinic orbit Γ as follows: we start at the saddle-node Then we follow the reduced (slow) flow towards the upper fold F + where we concatenate a fast fibre at F + that connects back towards the lower attracting branch S − a . Finally, we follow the reduced (slow) flow on S − a towards the lower fold F − and hence end up at the saddle-node equilibrium. This homoclinic orbit is the singular limit representation of the SNIC indicated in Fig. 1. Hence for a = c = 0, we have identified a (global) singular SNIC bifurcation together with a (local) singular Bogdanov-Takens bifurcation. The unfolding of these singular bifurcations is done in Sect. 6. Figure 5 summarises all our singular limit observations.

Type II Excitability: Singular Andronov-Hopf Bifurcation and Canard Explosion
In the case of type II excitability as shown in Fig. 3, the stable equilibrium on the lower branch S − a crosses the lower fold F − at I = I bif and moves onto the unstable middle branch S r where it becomes unstable as an equilibrium of the reduced problem (11). Hence, two eigenvalues change sign as we cross F − . The bifurcation point I bif will in general depend on λ; throughout this section, λ will be fixed and we will not further stress the dependence on λ, as the corresponding bifurcation is codimension-1 and only I is needed as a parameter for its unfolding. We introduce is a singular contact point that undergoes a singular Andronov-Hopf bifurcation with respect to the parameter I at I = I bif . More precisely, we impose Besides the singular point near F − occurring in this bifurcation, there are no other singular points on S − a .
Applying the normal form transformation outlined in the proof of Proposition 1 reveals that Assumption 5 amounts to imposing which can easily be verified on the canonical model (1).

Remark 6
System (1) satisfies Assumption 5 for a range of λ values, along a param- The nature of the singular Andronov-Hopf bifurcation can easily be seen when looking at the trace and the determinant of the Jacobian of system (3) given by Close to the fold F − , a bifurcation of equilibria defined by f = g = 0 happens for 0 < ε 1 when tr J = 0. This implies ∂f ∂v = −ε ∂g ∂w = O(ε) and, in the singular limit, this gives the fold condition ∂f ∂v = 0. From Assumption 5 we have ∂f ∂w ∂g ∂v < 0 evaluated at g = 0 and, hence, det J = O(ε) > 0. So, we are expecting a singular Andronov-Hopf bifurcation for I = I h that creates small O( √ ε) amplitude limit cycles with nonzero frequencies of order O( √ ε) [14,15,18]. Hence, the singular nature of the Andronov-Hopf bifurcation is encoded in both, amplitude and frequency. Figure 1 shows an example of a supercritical singular Andronov-Hopf bifurcation.
Note in Fig. 1 that the O( √ ε) branch of the Andronov-Hopf bifurcation suddenly changes dramatically near I = I c . This almost vertical branch marks the unfolding of canard cycles within an exponentially small parameter interval of the bifurcation parameter I . This is often referred to as a canard explosion [14,15]; it provides the necessary continuous connection between the small Andronov-Hopf limit cycles and the large relaxation cycles as shown in Fig. 1. In the singular limit, canard cycles can be identified as follows: Note that the stability switch of the equilibrium in the reduced problem (11) is due to the singular nature of system (11) at F − ; a stability switch of a single equilibrium without interacting with another equilibrium in a one-dimensional regular perturbation problem is otherwise not possible. In fact, for I = I bif there exists no equilibrium in the reduced problem (11) due to a cancellation of a simple zero. Hence, a trajectory is able to cross from S − a to S r with nonzero speed which is a hallmark of a singular canard. One can construct singular canard cycles that are formed through concatenations of slow canard segments and fast fibres as shown in Fig. 4(d). Note that these singular canard cycles have O(1) amplitude and have a frequency O(1) on the order of the slow time scale. These singular canard cycles will unfold to the above mentioned canard cycles in a canard explosion. The following summarises these observations. Theorem 1 In system (4) under Assumptions 1-5, assuming the existence of only one equilibrium and for sufficiently small ε, a singular Andronov-Hopf bifurcation and a canard explosion occur at The coefficients H 1 and K 1 and, hence, the type of singular Andronov-Hopf bifurcation (super-or subcritical), can be calculated explicitly. In fact, for the normal form (15) given in Proposition 1 and under the condition that c > 0 is sufficiently large, we find with respect to the bifurcation parameter a that a bif = 0, When K 1a > 0, the Hopf bifurcation is supercritical; when K 1a < 0, it is subcritical.
Proof We refer to [15], where the main part of the statement is shown. Here we just restrict to computing H 1a and K 1a in the normal form (15). Let us start with H 1a , given (21). A simple asymptotic analysis reveals that a singularity is located at ( , about which the linearisation of the vector field has a trace given by . The Hopf bifurcation hence occurs at H 1a = c 2 σ . Next we focus on the canard value H 1a + K 1a . Since the singular Hopf point is of generic nature, the parameter value at which canards are present are the same parameter values for which there exists a smooth asymptotic expansion x = ϕ(y) + ϕ 1 (y)ε + ϕ 2 (y)ε 2 + O(ε 3 ) representing an invariant graph. Expressing the invariance by plugging the series in the differential equations yields expressions for ϕ 1 and ϕ 2 , given a c = (H 1a + K 1a )ε + O(ε 2 ). Then imposing the requirement that ϕ 2 should not have a pole at y = 0 yields a condition on H 1a + K 1a that leads to the required result.
Remark 7 By actively using the singular nature of canards, the above calculation also presents an alternative way to find the first Lyapunov coefficient K 1a to determine the criticality of the singular Andronov-Hopf bifurcation.
In the singular limit, we have I h = I c = I bif indicating the singular nature of the bifurcation. Note that the classic definition of type II excitability refers to the slow O(ε) frequency band of the large relaxation oscillations which does not vary much (and not to the actual intermediate O(ε 1/2 ) singular Andronov-Hopf bifurcation frequency).
A closer look at the expression for the first Lyapunov coefficient K 1a in Theorem 1 indicates that there is a change of criticality at c = c bautin + O(ε 1/2 ) with which evaluates to c bautin = 2d 2 3 in the canonical model (1).
The Bautin bifurcation [17] indicates for fixed 0 < c < c bautin that we are dealing with a subcritical singular Andronov-Hopf bifurcation (under the variation of the parameter a) which is accompanied by another codimension-1 bifurcation, a saddlenode of periodic orbits (SNPO) bifurcation that has branched of the codimension-2 Bautin point. Due to the singular nature of our problem, this SNPO bifurcation is a bifurcation of canard cycles. Thus it happens exponentially close to the canard parameter value a c defined in Theorem 1.

Type I/II Excitability Transition Regime: Incomplete Canard Explosion
Recall from (17)  This lemma indicates that the cusp has no singular nature 8 with respect to the limit ε → 0. For fixed 0 < c < c cusp , we have three equilibrium states, at least two of which are located on S r (see Assumption 3). This changes the global bifurcation structure of our problem. While we still observe a singular Andronov-Hopf bifurcation with respect to the parameter a, the growth of the limit cycles is, however, bounded as one approaches a homoclinic connection towards one of these additional equilibria (which is of saddle type). The next theorem discusses this scenario, which describes a first transition from type II excitability towards the type I limiting situation. We formulate the results concerning the main system, but we will make the description according to the parameters introduced for the normal form (15) in Sect. 3.1. For the canonical model (1), the relation between a and I is trivial, while in general the relation between (I, λ) and (a, c) may be complicated, though Sect. 3.1 entails a procedure on how to compute the change of parameters.

Remark 8
We denote by c − sn > 0 the c-coordinate where one of the SN branches intersects the c-axis ('SNIC of canard type' in Fig. 5). In the normal form (15) excluding the higher order (big-oh) terms, we know that 0 < c − sn < c cusp . We assume that this is also the case with the big-oh terms included.

Theorem 2
In system (15) under Assumptions 1-5, for fixed 0 < c < c − sn < c cusp and 0 < ε 1 there exists an unstable equilibrium on the middle branch S r,ε bounded away from the lower fold F − . Furthermore, there exist functions 0 < a snpo (ε) < a (ε) < a s (ε) < a c (ε) < a h (ε) < a + sn (ε) Fig. 6 Bifurcation diagram for fixed 0 < c < c cusp . The parabola shows the two equilibra near the fold F − (the third equilibrium branch and second SN are not shown). The a-axis is not shown on scale, as the distance between a snpo and a c should be exponentially small. Keeping that in mind, an incomplete canard explosion is seen as a goes from a c to a s that all converge to zero in the singular limit ε → 0 (except a + sn ) and for which the following holds (see also Fig. 6

As a decreases from a , large-amplitude canard cycles appear that grow in amplitude until it disappears in a saddle-node bifurcation of limit cycles at a = a snpo .
Proof Before proving this theorem, we would like to mention that the distance between a h and a c is O(ε), just like in the previous singular Andronov-Hopf case, while the values a c , a s , a and a snpo are all exponentially close to each other.  (y, a). It yields At a = a + sn , the double singular point appears (in the singular limit ε = 0) at y = c 2 + O(c 2 ), i.e. on the middle branch, and for a h < a < a + sn there are two singular points: one unstable node/focus p − close to the fold, and a saddle p + that is O(c)away from the fold, both located on the middle branch.
For a > a h , the location of the singularities implies that the fold F − is a jump point and a stable relaxation cycle is present. This shows parts (1)-(3).

Remark 9
Note, there exists also a third equilibrium on the middle branch, an unstable node/focus denoted by n, bounded away from the fold F − . This third equilibrium bifurcates with p + along the second saddle-node branch (see Fig. 5).

Remark 10
The canard parameter value a c is not strictly defined, as also a "smallamplitude limit cycle" is not strictly defined. We choose a δ-neighbourhood of the fold and the moment where the canard cycles grow out of this δ-neighbourhood along the canard explosion, we define the parameter value a = a c (ε). (In other words, a c lies beyond the "birth of canards".) Incomplete canard explosion and homoclinic saddle loops. The presence of the saddle p + shows that the canard cycles cannot grow unlimitedly during the canard explosion. We give here a short overview of the proof of the presence of small canard cycles because we will use elements of the proof to show the existence of canard homoclinics.
For fixed c > 0, we rescale a = √ εA, thus studying an O( √ ε) neighbourhood of the positive c-axis in (a, c) parameter space along which we observe a singular Andronov-Hopf bifurcation; see Fig. 5. We consider, in normal form coordinates, two transverse sections: a section S = {y = 0, x 0 < x < x 1 } which is transverse to the fast flow and a section T = {y = 0, |x| = O(ε)} close to the singular fold F − . Using geometric singular perturbation theory, we find that both the forward flow and the backward flow of the vector field takes points of S to T ; indeed, the forward resp. backward fast flow takes points of S to S − a resp. S r , after which the dynamics of the slow flow governs the drift towards T . It defines two maps where the forward map F and the backward map B are known to be smooth in terms of (x S , √ ε, A); see [19]. To be more precise, let x + S be the (parameter dependent) coordinate of the intersection of S with the fast fibre separatrix of the saddle p + . Then this coordinate is the supremum of the x-coordinates for which the backward map in (26) is defined. In [19] it is then shown that F − B = 0 at ( √ ε, A) = (0, 0) (independently of x S ), and ∂ ∂A (F − B) = 0 at that point. Application of the implicit function theorem leads to the presence of a canard curve A = √ εA canard (x S , √ ε) along which the vector field has a canard periodic orbit that intersects the section S at x-coordinate x S . The canard curve A canard is smooth in ε and its value depends in an exponentially small way on x S . In other words, canard cycles are found in an exponentially small wedge along an O(ε) neighbourhood of the positive c-axis in (a, c)-parameter space.
We can further rely on the results in [20], where it is shown that the maps (26) are smooth up to and including (at its extension) the boundary x = x + S . This implies that the canard value a s = εA canard (x + S , √ ε) obtained above is actually a parameter curve along which the vector field has a homoclinic saddle loop of canard type (of 'jump-back' type). This proves part (7) of the theorem.

Remark 11
Around the unstable canard cycle (or homoclinic a bit later on), there appears a big relaxation oscillation of canard type. It lies close to the full relaxation oscillation, but travels an O(c)-distance along the middle repelling branch (see Fig. 6). The exact distance travelled along the middle branch can be computed using slowdivergence integrals and exit-entry relations; we refer to the literature [21].
By introducing an alternative sectionS between the middle and upper branch instead of S, we can treat homoclinic saddle loops of the 'jump-away' type in a completely similar way. The only thing that changes is that points ofS undergo a large-amplitude oscillation in their way to T in positive time (travelling along S + a towards the jump point F + , jumping off towards S − a ). The smoothness of the transition maps and the application of the implicit function theorem is analogous. This defines a = a < a s and proves part (9) of the theorem.
For a < a (ε), the homoclinic connection breaks into a repelling large-amplitude cycle. While a proceeds to the outside of an O(ε)-neighbourhood of 0, it encounters the relaxation-like attracting cycle that surrounded all repelling cycles. They disappear in a saddle-node bifurcation of limit cycles at a = a snpo < a . This proves part (10) of the theorem. Fig. 7 Heteroclinic connections of canard type undergo a transition from headless canard to canard with head, from the jump-back canard homoclinic to the jump-away canard homoclinic For a < a < a s in between the two homoclinic curves, the stable manifold W s connects in reverse time to an orbit that follows the large homoclinic for some time, but exiting at a time prior to the time needed to connect back to p + . Therefore, no limit cycles may appear. In fact, since we know there is an additional singularity on the middle branch S r , denoted by n, and assuming it is of node type, then the exponential gap between the two homoclinics is filled by canard curves along which canard-type heteroclinic connections appear between n and p + ; see Fig. 7. The transition from headless heteroclinic canard to heteroclinic canard with head can be seen as a natural continuation of the truncated canard explosion of the canard homoclinics. The proof of the presence of such heteroclinics is completely analogous to above. In particular, no limit cycles are present in this scenario. This proves part (8) and finishes the proof of the theorem.

Termination of Homoclinic Saddle Loops: SNICs
In the canonical form (1), the singular Andronov-Hopf curve intersects the saddlenode bifurcation curve along which the saddle p + collides with a third singularity, a node n on the middle branch S r . Expressing (1) in the local coordinates (14), this singular bifurcation point is marked in blue in Fig. 5 and has coordinates (a − sn , c − sn ) = (0, − 1 4β ). In this section, we discuss how this codimension-2 singular bifurcation point perturbs to ε > 0.
The SN-curve perturbs regularly for positive values of ε to a curve c = c sn (a, ε). Observe that the part of the AH curve between the origin and this codimension-2 bifurcation point perturbs to a wedge of canard curves, corresponding to the incomplete canard explosion discussed in Theorem 2.
As in the case of the saddle homoclinics, denote by x + S the intersection of the section S with the unstable separatrix of the singularity p + , up to the point where it becomes a saddle-node. Then the part {x < x + S } of S is the part for which the backward map B : S → T is well defined. Applying the results in [20], we know that the map B has a C k -smooth extension (for any k) to the boundary of its definition domain. It implies that we can define B(x + S , ε, a, c sn (a, ε)) and also F (x + S , ε, a, c sn (a, ε)), like in (26) but where we made the dependence of c explicit. Therefore, solving F − B = 0 using the implicit function theorem with respect to the rescaled parameter A (recall a = √ εA) allows us to prove the presence of saddle-node homoclinics of ('jumpback') canard type. As before, replacing the section S with the alternativeS, we can apply the same reasoning to prove the presence of a canard value along which there is a saddlenode homoclinic of ('jump-away') canard type. Since the canard values A canard are smooth in terms of x S , we clearly see that the homoclinic loop curves (of jumpback and jump-away type) terminate at corresponding saddle-node homoclinic at the saddle-node bifurcation curve.
The exponential wedge on the SN-curve between the two terminal points are terminal points of SNIC canard curves that correspond to heteroclinic canards as shown in Fig. 7. Visually, it is clear how in Fig. 7, a heteroclinic connection tends towards a SNIC as n and p + approach each other in a SN bifurcation. The method of proof is similar to the one exposed before.

Remark 12
In the case c − sn < c < c cusp fixed, the second SN-bifurcation is located ahead of the singular (subcritical) AH-bifurcation and we observe a complete canard explosion (including a SNPO bifurcation of canard cycles).
In the case c cusp < c < c bautin fixed, there is only the singular (subcritical) AHbifurcation and we also observe a complete canard explosion (including a SNPO bifurcation of canard cycles).

Type I Excitability: Singular Bogdanov-Takens/SNIC Bifurcation
In this section, we discuss the origin (a, c) = (0, 0) which represents a local singular Bogdanov-Takens and a global singular SNIC bifurcation point; see with G defined in (18). This condition is typically violated in 1-parameter families.
For SNIC bifurcations to appear with a saddle-node near the fold point F − , we will hence impose the following condition.

Assumption 6
For fixed (I, λ) = (I bif , λ bif ), the fold point F − = (w − , v − ) is a singular contact point that undergoes a singular Bogdanov-Takens bifurcation with respect to the parameters (I, λ) at (I, λ) = (I bif , λ bif ). More precisely, we impose (on top of Assumptions 1-4): Besides the possible singular points near F − occurring in this bifurcation, there are no other singular points on S − a .
Conditions (27) imply that the fold point F − is a local codimension-2 singular point. Conditions (28) imply that a complete unfolding of the singularity is obtained upon varying (I, λ). In fact, the conditions in (28) could be replaced by the slightly more general condition, where G v := ∂G/∂v, but we prefer to keep (28) in order to be able to identify I as the Hopf breaking parameter among the two parameters.

Remark 13
It can be seen that conditions (27) and (28) imply the following conditions on the normal form (15): which are verified on the canonical model (1).
Under these conditions it is well known that in ε-dependent rescaled coordinates, a regular Bogdanov-Takens bifurcation takes place; see [16]. As a consequence, the presence of small-amplitude homoclinics is clear in some parameter subset. Furthermore, as (singular) Andronov-Hopf bifurcations form part of the bifurcation diagram, canard-type orbits are present. Indeed, the double singularity in the slow dynamics at x = 0 may unfold in a way that the fold point becomes a canard point and an extra saddle-singularity in the slow dynamics on the middle branch may appear. In that way, an incomplete canard explosion can be observed that terminates in a canard-type saddle homoclinic ('jump-back', without 'head'). In fact, these are phenomena that appear locally near the Bogdanov-Takens fold point.
Besides the small-amplitude phenomena near the Bogdanov-Takens point, we consider orbits that are close to the singular saddle-node homoclinic loop Γ shown in Fig. 4. We expect the existence of large-amplitude saddle-node homoclinics (SNICs) and as in the previous section, we also expect large-amplitude saddle homoclinics as well as relaxation oscillations.
In order to get a hold on the parameters close to c = 0, we rescale the parameters and introduce for some large M > 0. By doing this we in fact assume that c = O(ε) and a = O(ε 2 ). After the parameter rescaling (29), we study the system The singularity at (x, y, ε) = (0, 0, 0) has been described in [16] as a slow-fast Bogdanov-Takens point. 9 In that paper, it is shown that a BT bifurcation takes place near the origin. More importantly, it is shown that the phase portraits associated with the BT bifurcation are the only phase portraits seen in a small neighbourhood of the origin. The paper does not deal with any interactions with global return mechanisms, i.e. the interaction of the (local) singular BT and the (global) singular SNIC have not been studied. We will therefore repeat part of the local analysis, with the focus on the interaction with the global return mechanism.

Blow-up of the Singular Fold
Near the fold, we study the system using blow-up [14,15]. We write where S 2 + denotes the half-sphere X 2 + Y 2 + E 2 = 1 with E ≥ 0 (also known as Poincaré or blow-up sphere). The weights are chosen in a way that the higher order (big-oh) terms in (30) have also higher order in the rescaled equation.
As is usual in geometric desingularisation, we study the flow on the half-sphere in different (coordinate) charts. Two charts are important: the chart K 1 (or the phasedirectional rescaling chart), and the chart K 2 (or the family rescaling chart). The K 1 chart is used to extend the orbits along the slow manifold (which are directed towards the fold) to a neighbourhood that is at distance O(ε) from the origin. While this Fig. 9 Blow-up of the singular fold. 'Birds-eye view' of the upper blown-up sphere X 2 + Y 2 + E 2 = 1 in (X, Y )-space. Normal hyperbolicity of the manifolds S a and S r is gained at the equator, allowing one to extend them onto the blown-up sphere near singularities p a , respectively, p r . The two additional singularities represent connection to the fast fibre at F − transition is the most technical and least obvious for the reader who is not accustomed to the blow-up method, it fortunately is that part where the study of system (30) agrees with the results in studies of a classical regular jump point. Hence, we do not present detailed computations in the chart K 1 , but focus on presenting important facts and refer to the literature [14,15] for a detailed analysis.
The phase-directional rescaling chart K 1 . Here, we explain the dynamics near the equator of the blow-up sphere S 2 + . When presenting a picture in blown-up (x, y, ε)space, where the origin is replaced by (or blown-up to) a sphere, we can position the point of view from the top of the ε-axis; looking down on the (X, Y )-plane we see the spherical surface X 2 + Y 2 + E 2 = 1 with E ≥ 0 as the interior of a circle X 2 + Y 2 = 1, and the equator as the circle with the outer slow-fast dynamics around it; see Fig. 9.
Calculations in K 1 reveal two hyperbolic saddle singularities, p s and p n , and two semi-hyperbolic singularities p a and p r along the equator. Combining information from the slow-fast dynamics near F − with information obtained in chart K 1 allows one to reconstruct the dynamics near the circle shown in Fig. 9.
Combining the global return mechanism, which defines a map Σ out → Σ in , with the information from this chart allows one to show the smoothness and exponentially contractiveness of the (ε, C, A)-family of maps Given the uniqueness of the centre separatrix issued from p a , one can prove that the image of any small section Σ n under this map limits to this centre separatrix (intersected with Σ a ) as ε → 0. To characterise the global dynamics, it is therefore important to know the dynamics of this centre separatrix. In particular, a connection of Σ a to Σ n will distinguish whether or not singular points are met, or whether a regular 'jump point-like' connection is possible. This study will be done in the family directional rescaling chart.
The family rescaling chart K 2 . Once the orbits have passed the chart K 1 , we can assume that x = O(ε 2 ) and y = O(ε). In the blow-up coordinates, this means that (X, Y, E) is bounded away from the equator {E = 0} of the sphere (since then it means ε ∼ r). It is well known that a study of that part of the sphere can be established by looking at an ε-dependent rescaling for some large R > 0. Applying this rescaling to (30), we can divide out a common factor ε, thus transforming the system into a regular perturbation family 10 (33)

Local and Global Codimension-2 Bifurcations
System (33) describes the flow in the interior of the sphere as shown in, e.g., Fig. 9, and it can be analysed by means of classic bifurcation analysis. Together with the information obtained from the global return mechanism, we are able to describe all observed local and global bifurcations in (A, C) parameter space.
Bogdanov-Takens bifurcation. Proof There are two singular points on X = Y 2 , located at Y = Y ± := 1 2 (C ± √ C 2 − 4A). The singularity at Y = Y + , denoted p + , is always a saddle. The singularity at Y = Y − , denoted p − is of focus/node type for C > 1 and undergoes a change in the sign of the trace along A h = − 1 4 + 1 2 C, which indicates an Andronov-Hopf bifurcation. Along this parameter line, p − is weakly unstable; a basic calculation shows that the first Lyapunov coefficient is positive. Hence the Andronov-Hopf bifurcation is subcritical.
The two singular points p + and p − collide along A + sn = C 2 /4 indicating a saddlenode bifurcation at p ± .
Saddle-node homoclinic bifurcation. The following proposition states some properties of system (33) for ε = 0. As mentioned before, we focus on the interaction of local dynamics with the global return mechanism (31). In particular, we want to understand the dynamics of the centre separatrix of p a .

Proposition 2
Along the saddle-node bifurcation line A + sn (C) = C 2 /4, we have the following behaviour of (33) for ε = 0 (see Fig. 11): 1. When C < 1 2 , the separatrix coming from p a connects to a centre-stable separatrix of the saddle-node singularity p ± . The unique unstable centre separatrix of p ± connects to p n . 2. When C = 1 2 , the separatrix coming from p a connects to the hyperbolic attracting separatrix of the saddle-node singularity p ± . The unique unstable centre separatrix of p ± connects to p n . 3. When C > 1 2 , the separatrix coming from p a connects directly to p n along a regular orbit. This is the jump scenario. In particular, the BT point p ± (C = 1) is not connected to the separatrix. 4. The attracting separatrix of the saddle-node point p ± and the separatrix coming from p a break regularly with respect to the parameter C.
Proof The proof uses basics from planar theory of vector fields (e.g. invariant curves, isoclines, positive invariant sets). It requires some computations, but as it concerns basic properties we have left the details for the reader.
Notice that the saddle-node singularity p ± is a point on the parabola W = 0 and thatẆ 2 ) 2 which is positive except at the SN point p ± . Using information from infinity (i.e. from chart K 1 ), we see that the separatrix from p a enters the region {W > 0} which is positively invariant. Hence the ω-limit set of the separatrix has to be the vertex of the parabola. Finally, the hyperbolic separatrix of the saddlenode singularity p ± is tangent to ∂{W = 0} which implies it lies outside the positive invariant set {W > 0}. This proves part (1).
For C = 1 2 , the singularity p ± lies on the invariant parabola from Lemma 3, which then coincides with the separatrix coming from p a . It is not hard to verify that the Fig. 11 Behaviour of (33) on the Poincaré disc for ε = 0, along the SN-curve A + sn (C) = C 2 4 . For C ≤ 1 2 , the point p a connects to the SN point p ± . At C = 1, a BT point p ± occurs, not connected to p a , however hyperbolic eigenspace of the saddle-node singularity p ± coincides with the tangent space of the parabola. This proves part (2).
It is a symbolic computation to verify that the separatrix coming from p a enters this invariant set, and hence cannot leave. 11 On the other hand, V computed at the saddle-node point p ± = ( 1 4 C 2 , 1 2 C) yields 1 4 (C − 1 2 ) > 0. We conclude that the separatrix from p a cannot reach p ± . Since p n is the only other option for a ω-limit, it shows part (3).
As for the regular breaking of the connection in part (4): we compute the stable separatrix of the saddle-node p ± and compare it with the separatrix coming from p a . For the comparison we choose an arbitrary section crossing {V = 0} and parameterise it by the levels of V . It is not hard to see that the separatrix coming from p a intersects any such section at V -values that are O(C − 1 2 ) 2 . On the other hand, a variational computation of the stable separatrix of p ± reveals that it is given by Since 1 4 e 1−4Y = 0, it explains the transversality. This finishes the proof of the theorem.
Clearly, the saddle-node curve A + sn (C) = C 2 /4 persists within a manifold A + sn (C, ε) = C 2 /4 + O(ε). The regular breaking property formulated in the proposition ensures that the results persist for ε > 0.  Behaviour of (33) on the Poincaré disc for ε = 0, along the curve A (C) = − 1 16 + C 4 . At C = 1 2 , the point p a connects to the SN point p ± . For C < 1 2 the p a connects to the node p − , for C > 1 2 , there is a sequence of two heteroclinic connections from p a to p n , via the saddle p + , which is resonant at C = 3 p + lies there. In the second case, the node p − is found to lie in {V > 0}. So, the unstable separatrix of the saddle in {V < 0} can only connect to p n . The computation of the eigenvalues is direct. Let us finally prove the regular breaking property, with a Melnikov-like approach that is adapted to planar dynamics; see [19]. DenotingẊ = f ,Ẏ = g, then . This function has a fixed sign on the separatrix from Y = −∞ up to the saddle p + at Y = C − 1 4 . We can now directly refer to [19] (Proposition 5.7) where the regular breaking is related to a Melnikov computation, where the integrand is exactly fg A − gf A (multiplied by a specific exponential that implies the convergence of the Melnikov integral). Since this function is sign-fixed, the related Melnikov integral is nonzero; see [23] for a generalisation of Melnikov theory to arbitrary dimensions.  (C, ε) and A (C, ε) are exponentially close.
Proof The presence of the homoclinic surface follows from a reasoning completely analogous to the one in the proof of Theorem 3. The change of stability is simply an eigenvalue computation: the equation ρ(C) = −1 is perturbed regularly under the ε-perturbation.
The emergence of an SNPO branch from the resonant saddle homoclinic is standard (see [24]), and based on three features: (i) the ratio of eigenvalues is perturbed regularly upon variation of a parameter (C), (ii) the separatrix connection breaks regularly upon variation of another parameter (A), and (iii) the divergence integral along the homoclinic loop is nonzero. Properties (i) and (ii) follow directly from the singular limit analysis in Proposition 3. Property (iii) follows from the slow-fast nature of the global return mechanism: the divergence integral computation is dominated by the passages along the slow branches S ± a , which are both attracting and yield a contribution of the order −K/ε, for some K > 0, while the fast parts and the parts near the folds yield an O(1) contribution. While this argument does not prove that the SNPO branch is uniformly defined up to the limit, the proof of such a result is based on combining the local Dulac map of the saddle with a return mechanism. Since all properties are uniform and since the global return mechanism is sufficiently smooth up to and including the singular limit, the method for showing SNPO branches is valid uniformly in ε.

Remark 15
There are no additional bifurcations (proof omitted).

Remark 16
The homoclinic surface A (C, ε) defined in Theorem 4 can be extended to C = 1 2 , up to and including its intersection with the SN-surface A + sn (C, ε) from Theorem 3. At the singular limit, this is seen in Fig. 10, but a proof is needed for ε > 0. In such a proof, one would need to blow up the vector field once more at the saddle-node singularity (x, y) = ( 1 16 , 1 4 ) and at the parameter value (a, c) = ( 1 16 , 1 2 ), using a family blow-up, in order to uniformly separate the saddle from the node. The technical issues involved in such a construction go beyond the scope of what we intend to expose in this paper.

Remark 17
The bifurcation curves AH, SNPO, HOM and HOM s shown in Fig. 8 and Fig. 10 are the same. To rigorously prove this, we would need to include the parameters (a, c) in the blow-up analysis, i.e. we would have to blow up the origin (x, y, ε, a, c) = (0, 0, 0, 0, 0). Again, the technicalities involved in such a construction go beyond the scope of what we intend to expose in this paper.

Discussion
Excitability is an important subject area in neuroscience and its modern treatment dates back to Alan Hodgkin's seminal work [1]. His distinction of three neural excitability classes based on injected steps of currents and observed corresponding distinct frequency-current (f-I) curves still forms the basis in understanding bifurcation mechanisms of neural excitability. FitzHugh was the first to use dynamical systems techniques for the qualitative description of action potential generation and threshold phenomena [7,8,25]. Rinzel and Ermentrout [3,4] provided then a mathematical framework based on bifurcation theory to distinguish between these excitability types: SNIC bifurcation for type I and Andronov-Hopf bifurcation for type II.
This dynamical systems approach pioneered by Rinzel and Ermentrout is also used to explain more complicated neural activity such as bursting patterns. Here, the inherent multiple time-scale structure of neural models given through inherent slow and fast cell membrane processes is actively used to explain the bursting pattern. The bifurcation structure found in the fast subsystem provides a possible key to understanding the genesis of bursting patterns; see, e.g., [5]. Interestingly enough, the spiking pattern itself within a burst is also often a result of a multiple time-scale structure. This relaxation type behaviour is typically ignored in the bursting literature, because it would mean to consider models with (at least) three time scales-fast, intermediate and slow-which has become only recently a new research focus [26][27][28].
On the other hand, the literature on neural excitability (see, e.g. [5]) clearly uses slow-fast decomposition to explain action potential generation for tonic spiking models, although the accompanying bifurcation analysis of such a system often ignores or just inconsequently uses the given slow/fast structure. This article actively explores the singular nature of (two-dimensional) neural models and identifies a novel singular bifurcation based on the slow-fast structure-the singular Bogdanov-Takens/SNIC bifurcation-that is key to understanding type I and (part of) type II excitability.
Using readily available tools and results from geometric singular perturbation theory [14][15][16][19][20][21][22] and bifurcation theory [17,24] we are able to unfold this singular bifurcation and identify important codimension-2 bifurcation points-Bautin, Bogdanov-Takens, resonant homoclinic and saddle-node homoclinic-which organise the bifurcation landscape for ε > 0 and help to explain the transitions between type I and type II excitability. For example, based on the position of the Bautin point in parameter space we identify a supercritical Andronov-Hopf bifurcation as a clear indicator of type II excitability while a subcritical Andronov-Hopf bifurcation does not necessarily guarantee a frequency band (significantly) bounded away from zero, i.e. the model could be close to type I and thus close to a homoclinic. Another important indicator for this proximity to type I is the cusp bifurcation and its corresponding saddle-node branches. Within the cusp region there are three equilibria including a saddle which is necessary to form a homoclinic loop. In this type I-II transition regime, properties of model neurons that are considered type II might show behaviour usually associated with type I and vice versa. Care has to be taken when inferring properties from such a simple excitability classification; see also [5], where many of these bifurcations and observations have been highlighted.
Our analytical bifurcation results provide important information for the computational neuroscience community. We show that a SNIC bifurcation associated with type I excitability only exists in a small parameter regime. Thus it is more likely to observe a (large) saddle homoclinic in one parameter continuation of neural models, although it might be very close to a saddle-node and, hence, be mistaken for a SNIC. Similarly, if one observes a subcritical Andronov-Hopf bifurcation close to a saddle-node, especially where the corresponding periodic orbits terminate nearby in a (small) homoclinic, then the other observed (large) homoclinic that terminates close to a saddle-node cannot be a SNIC.
A main obstacle in a numerical bifurcation analysis is not only the stiffness of the underlying problem but also the close proximity of different bifurcation branches. Our analytical bifurcation results should be seen as a helpful guide for numerical continuation. For example, the numerical bifurcation diagrams presented, e.g., for the Morris-Lecar neural model in [29], Fig. 6, or for a 2D sodium spiking model in [30], Figs. 1-2, are incomplete since the exponentially close branches HOM s , HOM , SNPO and SNIC (Fig. 8) are hard to distinguish numerically and the identification of certain codimension-2 points (Fig. 10)-Bogdanov-Takens, resonant HOM and SN-HOM -is a very difficult task.