A Novel Phase Portrait to Understand Neuronal Excitability

Fifty years ago, Fitzugh introduced a phase portrait that became famous for a twofold reason: it captured in a physiological way the qualitative behavior of Hodgkin-Huxley model and it revealed the power of simple dynamical models to unfold complex firing patterns. To date, in spite of the enormous progresses in qualitative and quantitative neural modeling, this phase portrait has remained the core picture of neuronal excitability. Yet, a major difference between the neurophysiology of 1961 and of 2011 is the recognition of the prominent role of calcium channels in firing mechanisms. We show that including this extra current in Hodgkin-Huxley dynamics leads to a revision of Fitzugh-Nagumo phase portrait that affects in a fundamental way the reduced modeling of neural excitability. The revisited model considerably enlarges the modeling power of the original one. In particular, it captures essential electrophysiological signatures that otherwise require non-physiological alteration or considerable complexication of the classical model. As a basic illustration, the new model is shown to highlight a core dynamical mechanism by which the calcium conductance controls the two distinct firing modes of thalamocortical neurons.


Introduction
Rooted in the seminal work of Hodgkin and Huxley [12], conductance-based models have become a central paradigm to describe the electrical behavior of neurons. These models combine a number of advantages, including physiological interpretability (parameters have a precise experimental meaning) and modularity (additional ionic currents and/or spatial effects are easily incorporated using the interconnection laws of electrical circuits [10,3]). Not surprisingly, the gain in quantitative description is achieved at the expense of mathematical complexity. The dimension of detailed quantitative models makes them mathematically intractable for analysis and numerically intractable for the simulation of large neuronal populations. For this reason, reduced modeling of conductancebased models has proven an indispensable complement to quantitative modeling. In particular, the FitzHugh-Nagumo model [9], a two-dimensional reduction of Hodgkin-Huxley model, has played an essential role over the years to explain the mechanisms of neuronal excitability (see e.g. [22,8] for an excellent introduction and further references). More recently, Izhikevich has enriched the value of reduced-models by providing the Fitzugh-Nagumo model with a reset mechanism [13] that captures the fast (almost discontinuous) behavior of spiking neurons. Such models are used to reproduce the qualitative [5,4] and quantitative [19,21] behavior of a large family of neuron types. Notably, their computational economy makes them good candidates for large-scale simulations of neuronal populations [6].
The Hodgkin-Huxley model and all reduced models derived from it [9,4] focus on sodium and potassium currents, as the main players in the generation of action potentials: sodium is a fast depolarizing current, while potassium is slower and hyperpolarizing. Initally motivated by reduced modeling of dopaminergic neurons in which calcium currents are essential to the firing mechanisms [7], the present paper mimicks the classical reduction of the Hodgkin-Huxley model augmented with an additional calcium current. The calcium current is a distinct player because it is depolarizing, as the sodium current, but acts on the slower timescale of the potassium current.
To our surprise, the inclusion of calcium currents in the HH model before its planar reduction leads to a novel phase portrait that has been disregarded to date. Mimicking earlier classical work, we perform a normal form reduction of the global HH reduced planar model. The mathematical normal form reduction is fundamentally different in the classical and new phase portrait because it involves a different bifurcation. The classical fold bifurcation is replaced by a transcritical bifurcation.
The results of these mathematical analysis lead to a novel simple model that further enriches the modeling power of the popular hybrid model of Izhikevich. A single parameter in the new model controls the neuron calcium conductance. In low calcium conductance mode, the model captures the standard behavior of earlier models. But in high calcium conductance mode, the same model captures the electrophysiological signature of neurons with a high density of calcium channels, in agreement with many experimental observations. For this reason, the novel reduced model is particularly relevant to understand the firing mechanisms of neurons that switch from a low calcium-conductance mode to a high calcium-conductance mode. Because thalamocortical (TC) neurons provide a prominent example of such neurons, they are chosen as the main illustration of the present paper, the benefits of of which should extend to a much broader class of neurons.

Planar reduction of Hodgkin-Huxley model revisited in the light of calcium channels
Calcium channels participate in the spiking pattern by providing, together with sodium channels, a source of depolarizing currents. In contrast to sodium channels whose gating kinetics are fast, calcium channels activate on a slower time-scale, similar to that of potassium channels [11]. As a consequence, they oppose the hyperpolarizing effect of potassium current activation, resulting in bidirectional modulation capabilities of the post-spike refractory period. We model this important physiological feature by considering the HH model [12] with an additional voltage-gated (noninactivating) calcium current I Ca and a DC-current I pump that accounts for hyperpolarizing calcium pump currents. The resulting model is similar to HH model when the conductance of the calcium currents is low, but becomes strikingly different when the calcium conductance is high. Note that the inactivation of calcium channels is not included in the HH dynamics because it is known to take place in a much slower time scale [25]. The inactivation will typically be accounted for by a slower adaptation of the calcium conductance, see Section 5. Figure 1a illustrates the spiking behavior induced by the action of an external square current I app in the two different modes. As compared to the original HH model (Fig. 1a left), the presence of the calcium current is characterized by a triple electrophysiological signature (see Fig. 1a right): • spike latency: the spike train (burst) is delayed with respect to the onset of the stimulation • plateau oscillations: the spike train oscillations occur at a more depolarized voltage than the hyperpolarized state • after-depolarization potential (ADP): the burst terminates with a small depolarization This electrophysiological signature is typical of neurons with sufficiently strong calcium currents. See for instance: spike latency [20,18], plateau oscillations [2], ADPs [1,5]. However, the mechanisms by which these behaviors occur have never been analyzed using reduced planar models to date. Following the standard reduction of HH model [9], we concentrate on the voltage variable V (that accounts for the membrane potential) and on a recovery variable n (that accounts for the overall gating of the ion channels) as key variables governing excitability (see methods). The phase-portrait of the reduced HH model is shown in Figure 1b (left). This phase portrait and the associated reduced dynamics are well studied in the literature (see [9] for the FitzHugh paper, and [8,4] for a recent discussion and more references). We recall them for comparison purposes  only. The resting state is a stable focus, which lies near the minimum of the familiar N-shaped V -nullcline. When the stimulation is turned on (spiking mode), this fixed point loses stability in a subcritical Andronov-Hopf bifurcation (see the numerical bifurcation analysis of Supplementary Section S.2), and the trajectory rapidly converges to the periodic spiking limit cycle attractor. As the stimulation is turned off (resting mode), the resting state recovers its global attractivity via a saddle-node of limit cycles (the unstable one being born in the subcritical Hopf bifurcation), and the burst terminates with small subthreshold oscillations (cf. Fig. 1a left).
In the presence of the calcium current, the phase-portrait changes drastically, as shown in Figure  1b (right). In the resting mode, the hyperpolarized state is a stable node lying on the far left of the phase-plane. The V -nullcline exhibits a "hourglass" shape. Its left branch is attractive and guides the relaxation toward the resting state after a single spike generation. The sign ofV changes from positive to negative approximately at the funnel of the hourglass, corresponding to the ADP apex. The right branch is repulsive and its two intersections with the n-nullcline are a saddle and an unstable focus.
When the stimulation is turned on, the V -nullcline breaks down in an upper and a lower branch. The upper branch exhibits the familiar N-shape and contains an unstable focus surrounded by a stable limit cycle, very much as in the reduced Hodgkin-Huxley model. In contrast, the lower branch of the V -nullcline, which is not physiological without the calcium currents, comes into play. While converging toward the spiking limit cycle attractor from the initial resting state, the trajectory must travel between the two nullclines where the vector field has smaller amplitude. As a consequence, the first spike is fired with a latency with respect to the onset of the stimulation, as observed in Figures 1a (right) in the presence of the calcium current (see also Supplementary Fig. S1). In addition, a comparison of the relative position of the resting state and the spiking limit cycle in Figure 1b (right) explains the presence of plateau oscillations. As the stimulation is turned off the spiking limit cycle disappears in a saddle-homoclinic bifurcation (see Supplementary Sections S.2 and S.3), and the resting state recovers its attractivity.
The presence of the lower branch of the V -nullcline has a physiological interpretation. In the reduced HH model, the gating variable n accounts for the activation of potassium channels and the inactivation of sodium channels. Their synergy results in a total ionic current that is monotonically increasing with n for a fixed value of V ( Supplementary Fig. S2 left). In this situation, at most one value of n solves the equationV = 0 and there can be only one branch for the voltage nullcline. In contrast, when calcium channels are present, the reduced gating variable must capture two antagonistic effects. As a result, the total ionic current is decreasing for low n (the gating variable is excitatory), and increasing for large n (the gating variable recovers its inhibitory nature) (Supplementary Fig. S2 right). In this situation, two distinct values of n solve the equationV = 0, which explains physiologically the second branch of the V -nullcline. To summarize, the lower branch of the voltage nullcline accounts for the existence of an excitatory effect of n, which is brought by calcium channel activation.

The central ruler of excitability is a transcritical bifurcation, not a fold one
The power of mathematical analysis of the reduced planar model (2) is fully revealed by introducing two further simplifications.
• Time-scale separation: we exploit that the voltage dynamics are much faster than the recovery dynamics by assuming a small ratioṅ = O(ε)V (the approximation holds away from the voltage nullcline) and by studying the singular limit ε = 0.
• Transcritical singularity: by comparing the shape of the voltage nullcline in Fig. 1(b) (I = 0) and Fig. 1(d) (I = 12), one deduces from a continuity argument that a critical value 0 < I c < 12 exists at which the two branches of the voltage nullcline intersect. The critical current I c depends on . In the singular limit ( = 0) and for the corresponding critical current I c = I c (0), one obtains the highly degenerate phase portrait in Figure 2a (center). This particular phase portrait contains a transcritical bifurcation (red circle) which is the key ruler of excitability. This is because, as illustrated in Figure 2b for > 0, the convergence of solutions either to the resting point (I < I c ) or to the spiking limit cycle (I > I c ) is fully determined by the stable W s and unstable W u manifolds of the saddle point. In the singular limit shown in Figure 2a, these hyberbolic objects degenerate to a critical manifold that coincides with the voltage nullcline near the transcritical bifurcation. It is in that sense that the X-shape of the voltage nullcline completely organizes the excitability, i.e. the transition from resting state to limit cycle.
The persistence of the manifold W s and W u away from the singular limit is proven by geometric singular perturbation (Supplementary Section S.3). The same analysis also establishes a normal form behavior in the neighborhood of the transcritical bifurcation: in a system of local coordinates centered at the bifurcation, the voltage dynamics take the simple forṁ where i is a re-scaled input current and with h.o.t. referring to higher order terms in v, w, ε.
It should be emphasized that it is the same perturbation analysis that leads to the classical view of the Hodgkin-Huxley reduced dynamics: the transition from Figure 1b It is of interest to realize that the addition of the calcium current in the HH model unmasks a global view of its phase portrait that has been disregarded to date for its lack of physiological relevance. Figure 3  bifurcation will be familiar to all readers of basic textbooks in bifurcation analysis. For instance, the sketch is found in [12] as a prototypical example of non-generic bifurcation. It is symptomatic that this particular example is described at length but not connected to any concrete model in a texbook that puts much emphasis on the relevance of bifurcation analysis in neurodynamics applications. As shown in Fig 3 (bottom), the missing connection is brought to life by calcium channels. Their particular kinetics renders the transcritical bifurcation of HH model physiological in the the presence of a high-conductance calcium current.

Transcritical hybrid modeling of neurons
The singular limit of planar reduced models reveals that the excitability properties of spiking neurons are essentially determined by a local normal form of bifurcation of the resting equilibrium. This property is at the core of mathematical analysis of neuronal excitability (see [8,4] and the rich literature therein). In recent work, Izhikevich showed that, for computational purposes, the combination of the local normal form dynamics with a hybrid reset mechanism, mimicking the fast (almost discontinuous) spike down-stroke, is able to reproduce the behavior of a large family of neurons with a high degree of fidelity [13,5]. Mimicking Izhikevich approach, we simplify the planar dynamics into the hybrid model:v The proposed transcritical hybrid model is highly reminiscent of the hybrid model of Izhikevich, but it consideraly enlarges its modeling power by including two features of importance: • the transcritical normal formv = v 2 − w 2 + I replaces the fold normal formv = v 2 − w + I, in accordance with the normal form analysis of Section 3. • the new parameter w 0 determines whether the intersection of the voltage and recovery nullclines will take place above (w 0 > 0) or below (w 0 < 0) the transcritical singularity.
The parameter w 0 is a direct image of the calcium conductance: for small calcium conductances, the recovery variable nullcline only intersects the upper branch of the voltage nullcline (Supplementary

Hybrid modeling of a thalamocortical relay neuron
Thalamocortical (TC) relay neurons are the input to sensory cortices. As illustrated in Figure 5, these neurons exhibit two distinct firing patterns: either a continuous regular spiking Fig. 5a, or a plateau burst spiking Fig. 5b. The switch between the two modes is regulated by prominent Ttype calcium currents that are deinactivated by hyperpolarization, thereby modulating the resting membrane potential [17]. Thalamocortical neurons are widely studied in the literature, and their two spiking modes make them a prototypical example to illustrate the relevance of the reduced model. We emphasize that our objective is not a fine tuned quantitative modeling of the TC neuron firing pattern. Rather, we attempt to provide a qualitative picture of how the proposed simple hybrid dynamics permits to explain the behavior of TC neurons and, in particular, the role of calcium currents. modes. As discussed in Section 4, the small calcium conductance mode is obtained by choosing a positive w 0 , whereas the large calcium conductance mode is obtained by choosing a negative w 0 , all the other parameters being identical in the two modes. The hybrid model reproduces the experimental observation: in the low-calcium mode, it responds with a slow regular train of action potentials; in the high-calcium mode, it responds with a long spike latency, plateau oscillations, and an ADP. Phase portrait analysis permits to explain the observed behavior (Supplementary Section SS.5). Among other, it shows that the switch from the upper branch of the v-nullcline (low calcium conductance mode) to its lower branch (high calcium conductance mode) is sufficient to reproduce and explain the strongly different step-responses of the TC neuron. In order to verify the physiological consistence of the transcritical model, we further compare its behavior with the simulated step response of a quantitative 200-compartments model of a TC relay cell 1 [6] in the large conductance mode (Fig.6). For the quantitative model, we plot the trajectory projection on the V − d plane, where V and d denotes the somatic membrane potential and the activation gating variable of the somatic T-type calcium current, respectively.
As shown on the figure, there is a striking similarity between the projection of the high dimension trajectory and the phase portrait of the second-order transcritical model. In both cases, the ADP is generated during a decrease of the activation variable, and plateau oscillations are exhibited far from the resting state. Moreover, the spike latency is a robust property of the transcritical model because the trajectory must visit the neighborhood of both the nullclinesV = 0 andẇ = 0 before converging to the spiking limit cycle. It should be stressed that there are no comparable ways to reproduce this behavior in a fold hybrid model. As shown by Izhikevich [4], reproducing this behavior with the standard reduced HH model necessitates a non physiological alteration of the reset rule (see also Supplementary Sections S.6 and S.7). This underlines the importance of the revisited model to capture the richness of neuron excitability. The phase-portrait of the compartmental model depicts the trajectory projection on the V − d plane, where V and d denotes the somatic membrane potential and the activation gating variable of the somatic T-type calcium current, respectively. In both cases, the ADP is generated during a decrease of the activation variable, and plateau oscillations are exhibited far from the resting state.

Conclusion
The inclusion of calcium channels in Hodgkin-Huxley model has a dramatic impact on its mathematical reduction: the firing mechanisms are governed by the local normal form of a transcritical rather than fold bifurcation. Interestingly, it is not the phase portrait of the reduced HH model that is affected by calcium, but only the subregion of the plane where it is physiologically relevant. As a consequence, the classical FitzHugh Nagumo phase portrait is a particular (because localized) view of the more complete picture studied in the present paper. Although this enlarged phase portrait is the source of rich and diverse forms of excitability, its essence is captured in a simple and physiologically grounded hybrid model. The illustration of its modeling power on the thalamocortical neuron excitability shows the impact of revisiting the classical view. This illustration is just the top of the iceberg because the same principle will apply to many important families of neurons that are thoroughly studied and that have so far largely resisted reduced modeling. The proposed model will impact the understanding of excitability of e.g. dopaminergic, serotonergic, and subthalamic nucleus neurons, whose various firing patterns have a direct and critical impact in physiology and diseases, such as Parkinson's disease and depression. Because of its simplicity and computational efficiency, it is also an ideal candidate for physiologically realistic studies of high-dimensional neuronal networks.

Equation and parameters of the complete model
The augmented HH model reads For the HH dynamics, we use the parameters of the original paper [12]. As all other HH currents, the additional calcium current obeys Ohm's law whereḡ Ca is the maximum calcium conductance, V Ca is the calcium Nernst potential, and d is the calcium activation gating variable. Such kinetics fit the majority of calcium channel subtypes [4]. Exploiting the similarities between the potassium and calcium gating kinetics, we further assume that the behavior of the calcium activation gating variable d is well approximated as a static function of the potassium activation gating variable n. The simple choice d := n, and a := 3 suffices for our paper. The exponent accounts for the delayed activation of calcium currents, described often empirically by second to sixth powers (see [11] page 112, and references therein). These parameters do not reflect any precise physiological calcium current. We choose them as a prototypical example. The functions α x , β x , x = n, m, h, can be found in the paper [12]. The value for the potassium Nernst potential V K = −12 is the same as in [12], while the sodium Nersnt and the leak Nernst potential are rounded to V N a = 120 and V l = 10.6, respectively. The values of the sodiumḡ N a = 120, potassiumḡ K = 36, and leak g l (maximum) conductances are the same as in [12]. The calcium Nernst potential is given by V Ca = 150. The numerical simulations of Figure  1

Planar reduction and phase portrait analysis
We follow the standard reduction of the original HH model to a two dimensional system by: i) assuming an instantaneous sodium activation, m ≡ m ∞ (V ), where m ∞ (V ) = α m (V )/(α m (V ) + β m (V )); ii) exploiting the approximate linear relation, originally proposed in [9], h + sn = c, with s, c 1. Applying the same reduction to (2) with parameters as above, we obtain the planar system V -"nullcline" refers to the set {(V, n) ∈ R 2 :V = 0} and similarly for other variables.

Hybrid modeling of TC neurons
For modeling convenience, we add a fitting parameter b ∈ R in the sub-threshold (continuous) voltage dynamics: The extra parameter does not affect the nature of the bifurcation, but tunes the slope of the v-nullcline branches. In order to account for the effect of intracellular calcium variations and for the associated activation of calcium pump currents, we also add a slow adaptation variable z as follows:

S.2 Bifurcation analysis
A bifurcation diagram of (5) with I app as the bifurcation parameter sheds more light on the transition mechanism between the resting and spiking modes (Figure 1). We use XPPAUT [1] for this numerical analysis. In both the I Ca -on and off modes, we draw the bifurcation diagram only for small I app , corresponding to the transition from resting to limit cycle oscillations (for larger I app , the stable limit cycle disappears in a supercritical Andronov-Hopf bifurcation in both cases, which leads to a stable depolarized, i.e. high-voltage, state).  Beyond the bifurcation, the trajectory converges to the stable spiking limit cycle. When I app is lowered again below I SN LC , the spiking limit cycle disappears in a saddle-node of limit cycles, the unstable one (not drawn) emanating from the subcritical Andronov-Hopf bifurcation, and the trajectory relaxes back to rest. Figure S4 (right) illustrates the bifurcation diagram of the reduced Hodgkin-Huxley model in the presence of calcium channels. For I app < I SN , a stable node (lower branch), a saddle (central branch), and an unstable focus (upper branch) are present, as in Figure 1b(top right). The node and the saddle coalesce in a supercritical fold bifurcation at I app = I SN , and disappear for I app > I SN , letting the trajectory converge toward the stable limit cycle. The spike latency observed in the I Ca -on configuration unmasks the ghost of this bifurcation. The stable limit cycle disappears in a saddle homoclinic bifurcation as I app falls below I SH , which lets the trajectory relax back to the hyperpolarized state.
The homoclinic bifurcation exhibited by the Hodgkin-Huxley model with calcium channels is a key mathematical difference with respect to the standard HH model. The next section unfolds this bifurcation by exploiting the time-scale separation between the fast voltage V and the slow recovery variable n.

S.3 Transcritical and saddle-homoclinic bifurcations of in the presence of calcium channels
By exploiting the sharp time-scale separation between the membrane voltage and the recovery variable dynamics, in this section we assume thatṅ = O( )V , where > 0, and study the singularly perturbed limit = 0. We highlight in particular the existence of a transversal V -nullcline selfintersection in the passage from Figure 1b(top right) to Figure 1b(bottom right) that is associated to a transcritical singularity (see also Fig. 2a center). We then rely on some recent results [9] to qualitatively study the phase portrait of (2) away from the singular limit (i.e. > 0) and for different values of the input current.

S.3.1 Transcritical bifurcation in the presence of calcium channels
The presence of the V -nullcline self intersection is confirmed by a closer inspection of the phase portrait near the nullcline break-up, as illustrate in Figure S5. Let us discuss this transition in more details. The existence of the self intersection follows directly by ∂V ∂Iapp = 1 and the implicit function theorem [11]. We claim that it is also transversal. Let I tc ∈ (2.4, 2.46) be the input current for which the V -nullcline has the self-intersection. Let f V (V, n, I app ) be the voltage dynamics of the reduced model (2) in the presence of calcium current, that is whereV is defined in (2) and all the other parameters are unchanged. Then, we want to show that there exists (V tc , n tc ) ∈ R 2 , such that ∂f V ∂n (V tc , n tc , I tc ) = 0 (S5c) describing a non-degenerate (transversal) self-intersection of the V -nullcline at (V tc , n tc ) (see e.g. [9] and Section 5.5.2, Th. 5.7 in [12]). The self-intersection satisfies the additional constraint that both intersecting branches are not parallel to the V axis, as implied by (S5e). Condition (S5a) is the V -nullcline equation. Noticing that ∂f V ∂n and ∂f V ∂V do not depend on I app , conditions (S5b) and (S5c) follows by the fact that, as I app < I tc varies, as in Figure S5(a), the right and left extrema of, respectively, the left and right branches of the V -nullcline lie, by definition, on the line ∂f V ∂n = 0, and, similarly, as I app > I tc varies, as in Figure S5(b), the minimum and the maximum of, respectively, the upper and lower branches of the V -nullcline lie, by definition, on the line ∂f V ∂V = 0. Since at the intersection the four extrema coincide, conditions (S5b) and (S5c) follow. We stress that (S5b) and (S5c) define the point (V tc , n tc ) as the intersection of two lines (cf. Figure S5) in a unique way. Conditions (S5d) and (S5e) are generic and can be easily verified numerically.
In the singular limitṅ = 0, the self intersection described by conditions (S5) corresponds to a transcritical bifurcation (see e.g. Section Section 3.2 in [13]) of the voltage dynamics, with n as the bifurcation parameter. That is, as sketched in Figure 2a(center), for I app = I tc , the two intersecting branches exchange their stability at n = n tc . More precisely, as n varies above and below the intersection point, i.e. n = n tc , the voltage dynamics has two fixed points, the left one is stable, whereas the right is unstable. The stability of the fixed points follows from the fact that, on the left of the self-intersection, ∂f V ∂V < 0, while ∂f V ∂V > 0 on the right (cf. Figure S5). The two fixed points collide in a transcritical singularity at the self-intersection.

S.3.1.1 A normal form lemma
In this technical section we compute a normal form of (2) associated to the self-intersection described by conditions (S5).
Give > 0, let whereṅ is defined in (2), and g 0 := g(V tc , n tc ), where V tc and n tc satisfy the defining conditions (S5b)-(S5b). The form (S6) is just a rescaling (through ) of the recovery variable dynamics that highlights the timescale separation between V and n. The following lemma is an application of Lemma 2.1 in [9] to the reduced HH model (2) with calcium current.
Lemma 1 Suppose that g 0 < 0. Then there exists an affine change of coordinates and a rescaling of and I that transforms the dynamical systeṁ where f V and g are respectively defined in (S4) and (S6), into the dynamical systeṁ  w, )).
Noticing that f V is affine in the input current, the extra term α(I − I tc ) must be added tov in the case I = I tc . The result follows by defining the rescaled input currentĨ := λ + α(I − I tc ).
Note that for the dynamical system (S9) with I = 0, the self intersection of the V -nullcline discussed in Section S.3.1 becomes evident, with the two intersecting branches given by w = ±v.

S.3.2 Singularly perturbed saddle-homoclinic bifurcation in the presence
of calcium channels Figure S6 provides a close look of the saddle-homoclinic bifurcation in the reduced HH model (2) with calcium current. In order to rigorously prove the existence of this bifurcation, we rely on the normal form introduced in Lemma 1 and exploit the timescale separation between v and w through geometrical singular perturbations theory (GSPT). To this aim we must briefly recall some basic of GSPT, relying on (S9) as an explicit example. The interested reader will find in [7] an excellent introduction to the topic, and in [10,9,8] some recent extensions on which we rely for the forthcoming analysis. The time rescaling τ := t transforms (S9) into the equivalent system which describes the dynamics (S9) in the slow timescale τ . In the limit = 0, commonly referred to as the singular limit, one obtains from (S9) and (S10) two new dynamical systems: the reduced dynamics evolves in the slow timescale τ , while the layer dynamicṡ evolves in the fast timescale t. Figure S7(a) depicts the fast-slow dynamics associated to (S11)-(S12). The main idea behind GSPT is to combine the analysis of the reduced and layer dynamics to derive conclusion on the behavior of the nominal system, i.e. with > 0. The reduced dynamics (S11) is a dynamical system on the set S 0 := {(v, w) ∈ R 2 : v 2 − w 2 + I + h 1 (v, w, ) = 0}, usually called the critical manifold. The points in S 0 are indeed critical points of the layer dynamics (S12). More precisely, portions of S 0 on which ∂V ∂V is non-vanishing are normally hyperbolic invariant manifolds of equilibria of the layer dynamics, whose stability is determined by the sign of ∂V ∂V . Conversely, points in S 0 where ∂V ∂V = 0 constitute degenerate equilibria. In particular, the layer dynamics (S12) exhibits, for I = 0, two degenerate equilibria. As depicted in Figure S7(a), they are given by the self-intersection of the V -nullcline, corresponding to a transcritical singularity (Section S.3.1), and by the fold singularity at the maximum of the upper branch of the V -nullcline.
The basic result of GSPT, due to Fennichel [2], is that, for sufficiently small, non-degenerate portions of S 0 persist as nearby normally hyperbolic locally invariant manifolds S of (S8). More precisely, the slow manifold S lies in a neighborhood of S 0 of radius O( ). The dynamics on S is a small perturbation of the reduced dynamics (S11). We point out that S may not be unique, but is determined only up to O(e −c/ ), for some c > 0. That is, two different choices of S are exponentially close (in ) one to the other. Since the presented results are independent of the particular S considered, we let this choice be arbitrary. The trajectories of the layer dynamics perturb to a stable and an unstable invariant foliations with basis S .
The analysis near degenerate points is more delicate. Only recently some works have treated this problem in its full generality for different types of singularities [10,9,8]. Figure S7 (b),(c),(d) sketch the extension of the attractive slow manifold S a after the fold point, and the three possible ways in which S a and the repelling slow manifold S r can continue after the trascritical singularity, depending on the injected current.
The result depicted in Figure S7 relies on the following theorem, adapted from [9]. Figure S7, where ρ < 0 and |ρ| is sufficiently small, and v − , v + are such that ∆ ∩ S − r = ∅. For a given > 0, let q a, := ∆ ∩ S a and q r, := ∆ ∩ S r be the intersections, whenever they exist, of respectively the attractive and repelling invariant submanifolds S a and S r with the section ∆. The following theorem reformulates in a compact way the discussion contained in Remark 2.2 and Section 3 of [9] 2 for systems with inputs of the form (S9). Theorem 1 implies the existence of the saddle-homoclinic bifurcation in the reduced Hodgkin-Huxley model (2) with calcium current. As stressed by Figure S7(b,c,d), the slow attractive S a (resp. repelling S r ) manifold coincides with the unstable W u (resp. stable W s ) manifold of the saddle point, as it can be proved via qualitative arguments. Thus, for I < I c ( √ ), the unstable manifold W u continues after the transcritical singularity on the left of W s , toward the stable node. See Figure S6(a,b) and Figure S7(b). For I = I c ( √ ), W u extends after the transcritical point to W s , forming the saddle-homoclinic trajectory, as depicted in Figure S6(c) and Figure S7(c). For I > I c ( √ ), the unstable manifold of the saddle W u continues after the transcritical singularity on the right of W s , and spirals toward an exponentially stable limit cycle, whose existence can be proved with similar GSPT arguments (see for instance [10]). This situation is the one depicted in Figure S6(d) and Figure S7(d).

S.4 Hybrid singularly perturbed saddle-homoclinic bifurcation
The saddle-homoclinic bifurcation analysis provided in Section S.3.2 for the reduced calcium-gated HH model naturally extends to the hybrid dynamics (1) with w 0 < 0. Indeed, by construction, if w 0 < 0, this model can be transformed in the normal form (S9) derived in Lemma 1, and Theorem 1 applies directly.
Similarly to the derivation in Section S.3.2, we can associate to (1) two new dynamical systems, describing its singular dynamics in the fast timescale t and in the slow timescale τ := t, respectively. More precisely, the hybrid reduced dynamics evolves in the fast time-scale t. Figure S8(a) depicts the associated slow-fast hybrid dynamics for I = 0, a ∈ (0, 1), and w 0 < 0. The analysis of the non-singular limit follows the same line as the analysis developed in Section S.3.2 for the continuous time case. The only difference is that the return mechanism, provided in the continuous time case by the right attractive branch of the critical manifold S 0 , is now replaced by the hybrid reset. The result is summarized in Figure S8   coincides with the image of the unstable manifold of the saddle W u through the hybrid reset map. As in the continuous time case, the stable manifold of the saddle W s can be shown to coincide with the slow repelling manifold S r . Let I c ( √ ) be defined as in Theorem 1. For I < I c ( √ ), the unstable manifold W u is brought back on the left of W s , and the resting state is globally stable, as in Figure  S8(b). At I = I c ( √ ), the hybrid reset connects W u and W s , corresponding to a hybrid homoclinic bifurcation. The associated phase-portrait is shown in Figure S8(c). Finally, for I > I c ( √ ), the unstable manifold W u is brought by the reset mechanism on the right of W s and directly into the newborn hybrid limit cycle attractor, as in Figure S8(d).
We emphasize that the existence of the hybrid homoclinic bifurcation and the associated critical value I c are independent of the reset point (c, d), provided that the trajectory starting from (c, d) be attracted toward the left upper branch of the critical manifold, as in Figure S8. As discussed in Section S.3.2, under this condition, two different choices of the reset point are associated to two slow manifolds S a that are O(e −c/ ) near, for some c > 0. By the result in Theorem 1, the two values of the critical current I c for which the slow attractive manifold extends to the slow repelling manifold are again O(e −c/ ) near. This also ensures that the value for which the hybrid homoclinic bifurcation happens is independent of the reset point, modulo variations that are O(e −c/ ). The same robustness properties are not shared by fold hybrid models (see Section S.7 below).
S.5 Phase-plane analysis of the transcritical hybrid TC neuron model  (3) in the low and high calcium conductance modes, as well as the corresponding phase portraits. The analysis of the phase-portraits gives insights on the mechanisms by which TC cells exhibit tonic firing or bursting when submitted to a similar step input, according to the initial resting potential. When w 0 > 0 (Fig. S9, left), which corresponds to a small calcium conductance (T-type calcium channels are inactivated), the hyperpolarized state belongs to the upper branch of the v-nullcline.
Application of a depolarizing current step lifts the voltage nullcline above the resting state, thus generating a transient non-delayed action potential (marked with a * in Fig. S9b, left). Note that, when the hyperpolarized state belongs to the upper branch of the v-nullcline, no plateau oscillations are possible (Fig. S9b, left). Furthermore, the relaxation toward the hyperpolarized state is necessarily monotone (i.e. no ADP), as stressed in Fig. S9c, left.
At the generation of the first spike, a certain amount of calcium ions enter the cell due to the presence of high voltage activated calcium currents. The subsequent activation of outward calcium pump currents 3 transiently hyperpolarizes the cell. As the intracellular calcium is expelled (Fig. 9c, left), calcium pump currents slowly deactivate and the cell slowly depolarizes (interspike period), until the spiking threshold is reached and a new action potential is fired.
When w 0 < 0 (Fig. S9, right), which corresponds to a large calcium conductance (T-type calcium channels are available), the hyperpolarized state belongs to the lower branch of the vnullcline, and is more hyperpolarized than for w 0 > 0, as in experiments. When a depolarizing current step is applied, this branch falls below the w-nullcline (Fig. S9b, right). In order to generate the first spike, the state travels in the narrow region between the two nullclines, resulting in a pronounced latency. Furthermore, the relative position of the hyperpolarized state (lower branch) with respect to (hybrid) spiking cycle (upper branch) clearly explains the generation mechanism of plateau oscillations. These high frequency plateau oscillations (burst) continue until the intracellular calcium, which enters at each action potential, and the associated pumps current (i.e. z) are sufficiently large. Plateau oscillations then terminate in a (hybrid) saddle-homoclinic bifurcation (Fig. S9c, right). At the end of the burst, the system converges toward the hyperpolarized state following the left branch of the v-nullcline, thus generating a marked ADP at the passage near the nullcline's funnel. The subsequent slow phase is mainly ruled by the variations of the intracellular calcium. With the adopted simple dynamics it consists in a slow depolarization that follows the decrease of the intracellular calcium (Fig. S9c, right). A finer and more physiological modeling of the intracellular calcium dynamics could reproduce in vitro recordings with a higher degree of fidelity.

S.6 Modeling TC neurons with a fold hybrid model
This section reproduces the modeling of TC neuron with the fold hybrid model. Such a model was recently proposed by Izhikevich and used in a large-scale model of the mammalian thalamic system [6]. The model succeeds in reproducing the firing pattern of TC neurons by modifying the reset mechanism. Fig. S10 shows that the voltage time-course reproduces with fidelity the three hallmarks of large calcium conductances, that is (A) spike latency, (B) plateau oscillations, and (C) ADP. However, the phase portrait illustrates that the modified v-reset law has no obvious physiological interpretation. As a consequence, plateau oscillations are not endogenously generated in this model. Furthermore, as opposed to both the quantitative model and the proposed hybrid model, the ADP is generated following an increase of the recovery variable, when the reset point crosses the stable manifold of the saddle.

S.7 Robust ADP generation
For a neuron model to be biologically relevant, it should be robust to exogenous disturbances (small synaptic inputs, thermal noise, etc.). The firing pattern, in particular, should remain unchanged. Figure S11 compares the perturbation robustness of three TC neuron models to small current impulses. It suggests that the fold hybrid model is less robust than the transcritical model, because a tiny pulse is sufficient to generate an extra action potential at the ADP apex.
The difference in robustness is explained by the different ADP generation mechanisms, as illustrated in Figure S12. In the fold model, ADPs are generated when trajectories cross the v-nullcline from below (cf. Figures S10 and S12). The absence of any robust attractor in the ADP generation region makes the ADP height and shape heavily dependent on the exact reset point. Moreover, when small current pulses are applied, the ADP generation is disrupted, and the model fires an extra (non-physiological) spike.
Conversely, ADP generation in the transcritical hybrid model is robustly governed by the attractor S a that steers the trajectories through the ADP apex and toward the resting point. That is the reason why the ADP height and shape barely depend on chosen reset point. At the same time, the persistence to small perturbations of this invariant manifold [3] ensures, as required in biologically meaningful conditions, the robustness of the ADP generation mechanism to small inputs.    Figure S12: Comparison of the ADP generation mechanisms in the fold (left) and in the transcritical hybrid models (right). The stable manifold of the saddle (×) is depicted in green. In the Izhikevich model ADPs are generated by sliding near the stable manifold of the saddle and crossing the v-nullcline from below. In our proposed hybrid model, ADP are robustly generated along the attractor S a .