A Novel Phase Portrait for Neuronal Excitability

Fifty years ago, FitzHugh 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 a 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 FitzHugh-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 complexification of the classical model. As a basic illustration, the new model is shown to highlight a core dynamical mechanism by which calcium channels control the two distinct firing modes of thalamocortical neurons.


Introduction
Rooted in the seminal work of Hodgkin and Huxley [1], 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 [2,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 conductance-based models has proven an indispensable complement to quantitative modeling [4][5][6]. In particular, the FitzHugh-Nagumo model [7], 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. [8,9] 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 [10] that captures the fast (almost discontinuous) behavior of spiking neurons. Such models are used to reproduce the qualitative [11,12] and quantitative [13,14] behavior of a large family of neuron types. Notably, their computational economy makes them good candidates for large-scale simulations of neuronal populations [15].
The Hodgkin-Huxley model and all reduced models derived from it [7,12] 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 [16], 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.
The inclusion of calcium currents in the HH model before its planar reduction leads to a novel phase portrait that seems to have 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 a proof of concept of the present paper, the benefits 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 [17]. As a consequence, their activation opposes 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 [1] with an additional noninactivating voltage-gated calcium current I Ca and a DC-current I pump that accounts for hyperpolarizing calcium pump currents. The inactivation of calcium channels occurs in a slower time scale than the HH dynamics [18]. It can be modeled by a slower adaptation of the calcium conductance, which does not affect the single spike generation mechanism. Figure 1 A 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. 1 A left), the presence of the calcium current is characterized by a triple electrophysiological signature (see Fig. 1 A right): N spike latency: the spike train (burst) is delayed with respect to the onset of the stimulation N plateau oscillations: the spike train oscillations occur at a more depolarized potential than the hyperpolarized state N 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 [19,20], plateau oscillations [21], ADPs [22,23]. 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 [7], 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 Fig. 1 B (left). This phase portrait and the associated reduced dynamics are well studied in the literature (see [7] for the FitzHugh paper, and [9,12] 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 below), 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. 1 A left).
In the presence of the calcium current, the phase-portrait changes drastically, as shown in Fig. 1 B (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 of _ V V 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 Vnullcline, 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 Fig. 1 A (right) in the presence of the calcium current (see also Fig. S1). In addition, a comparison of the relative position of the resting state and the spiking limit cycle in Fig. 1 B (right) explains the presence of plateau oscillations. As the stimulation is turned off the spiking limit cycle disappears in a saddle-homoclinic bifurcation (see below), 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 (Fig. 2, left). In this situation, at most one value of n solves the equation _ V V~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) (Fig. 2, right). In this situation, two distinct values of n solve the equation _ V V~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.
A bifurcation diagram with I app as the bifurcation parameter sheds more light on the transition mechanism between the resting and spiking modes (Fig. 3). We use XPPAUT [24] for this numerical analysis. We draw the bifurcation diagram without (left) and with calcium channels (right) only for small I app , corresponding to the transition from resting to limit cycle oscillations (Fig. 3) (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). Fig. 3 (left) illustrates the bifurcation diagram of the original reduced Hodgkin-Huxley model. For low values of I app , the unique fixed point is a stable focus that loses stability in a subcritical Andronov-Hopf bifurcation at I app~IAH . Beyond the bifurcation, the trajectory converges to the stable spiking limit cycle. When I app is lowered again below I SNLC , 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.  Fig. 1 B(top right). The node and the saddle coalesce in a supercritical fold bifurcation at I app~ISN , and disappear for I app wI 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 Central Ruler of Excitability is a Transcritical Bifurcation, not a Fold One
The power of mathematical analysis of the reduced planar model (1) is fully revealed by introducing two further simplifications.
N Time-scale separation: we exploit that the voltage dynamics are much faster than the recovery dynamics by assuming a  The critical current I c depends on e. In the singular limit (e~0) and for the corresponding critical current I c~Ic (0), one obtains the highly degenerate phase portrait in Fig. 4 A (center). This particular phase portrait contains a transcritical bifurcation (red circle) which is the key ruler of excitability. This is because, as illustrated in Fig. 4 B for ew0, the convergence of solutions either to the resting point (IvI c ) or to the spiking limit cycle (IwI c ) is fully determined by the stable W s and unstable W u manifolds of the saddle point. In the singular limit shown in Fig. 4 A, 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 can be rigorously established by geometric singular perturbation. The details of this analysis are available in the report [25]. 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 form where i is a re-scaled input current and with h:o:t: referring to higher order terms in v,w,e.
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 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 5 A shows the phase portrait of the classical reduced HH  Likewise, the conceptual sketch of the transcritical bifurcation will be familiar to all readers of basic textbooks in bifurcation analysis. For instance, the sketch is found in [26] 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. 5 B, the missing connection is brought to life by calcium channels. Their particular kinetics renders the transcritical bifurcation of HH model physiological in the presence of a highconductance 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 [9,12] 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 [10,11]. Mimicking Izhikevich approach, we simplify the planar dynamics into the hybrid model: 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: N the new parameter w 0 determines whether the intersection of the voltage and recovery nullclines will take place above (w 0 w0) or below (w 0 v0) 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 (Fig. S2); likewise in the hybrid model when w 0 w0. For high calcium conductance, the recovery variable nullcline intersects the lower branch of the voltage nullcline (Fig. S2); likewise in the hybrid model when w 0 v0. Fig. 6 summarizes the four different phase portraits that derive from the transcritical hybrid model for different values of I and w 0 . For w 0 w0, the model captures the classical view of the reduced HH model Fig. 6(bottom). For w 0 v0, the model reveals the novel excitability properties associated to a high calcium conductance Fig. 6(top).The reader will notice the similarity between the phase portraits of Fig. 1

Reduced Modeling of a Thalamocortical Relay Neuron
Thalamocortical (TC) relay neurons are the input to sensory cortices. These neurons exhibit two distinct firing patterns: either a continuous regular spiking (Fig. 7 A) or a plateau burst spiking (Fig. 7 B) [27][28][29][30]. The switch between the two modes is regulated by prominent T-type calcium currents that are deinactivated by hyperpolarization, thereby modulating the resting membrane potential [28][29][30][31][32]. These firing patterns have been observed during both in vitro and in vivo recordings [33][34][35][36][37], and have been shown to play an important role in thalamocortical relay [38][39][40][41]. Among others, synchronous burst firing of TC relay cells is the key component of slow-wave sleep [40,42], whereas the pathological generation of this firing pattern during wakefulness leads to absence epilepsy [43][44][45].  Because TC relay cells have an important role in physiology and pathology, they have been widely studied in the literature. In particular, their electrophysiological activity has been successfully reproduced in various conductance-based models [18,29,30,[46][47][48]. The two distinct spiking modes of TC neurons and their dependence on the activation of calcium channels make them a good candidate to test the relevance of our qualitative reduced model against quantitative conductance-based models. In order to verify this hypothesis, we perform a classical two-dimension reduction of a state-of-the-art conductance based model proposed in [46] (see Methods for the reduction details). Fig. 8 shows the phase portraits of the reduced TC model in the two different modes. When the neuron is initially depolarized, Ttype calcium channels are inactivated, and the phase portrait is similar to the traditional FitzHugh-Nagumo model: the voltage nullcline has a simple (upper) branch, which breaks into a left and a right branches when the stimulation is off (Fig. 8 A, black solid curves). The accompanying transcritical singularity appears at non-physiological values of the gating variable m T . When the stimulation is turned on, the v-nullcline upper branch rises (Fig. 8 A, light gray solid curve), the fixed point looses its stability to a limit cycle that relaxes to the resting state when the stimulation is switched off. Likewise, no plateau is observed and the repolarization phase is strictly monotonic (no ADP can be exhibited).
The situation is very different when the neuron is initially hyperpolarized, because T-type calcium channels are then deinactivated, critically affecting the phase portrait: a lower vnullcline branch is now present for physiological values of m T , and the hyperpolarized state belongs to this lower branch (Fig. 8 B, upper-left, black gray solid curves). When a depolarizing current step is applied, this branch falls below the m T -nullcline (Fig. 8 B, upper-left, light gray solid curves). In order to generate the first spike, the state travels the narrow region between the two nullclines, resulting in a pronounced latency. This latency relies on the small level of I T activation in this hyperpolarized state (the m Tnullcline is almost horizontal), and is amplified by the dynamic inactivation of T-type calcium channels, which further narrows this funnel. This observation in agreement with previous experimental and modeling data [48]. Furthermore, the relative position of the hyperpolarized state (lower branch) with respect to spiking cycle (upper branch) clearly explains the generation mechanism of plateau oscillations.
These high frequency plateau oscillations continue until the Ttype calcium channel inactivation dominates (Fig. 8 B, upperright). At the end of the burst, the system converges toward the hyperpolarized state, being attracted first by the upper part of the v-nullcline (depolarizing phase of the ADP), then by its lower part (repolarizing phase of the ADP) (Fig. 8 B, lower-left). Note the presence of a transcritical bifurcation in the phase portrait for a particular value of T-type calcium channel inactivation (Fig. 8 B, lower-left, green circle). Finally, when the stimulation is relaxed, the neuron recovers its initial hyperpolarized state and the inactivation of T-type calcium channels is released (Fig. 8 B,  lower-right).
This phase portrait analysis confirms that the transcritical singularity in indeed a key ruler of excitability in this reduced conductance based TC neuron model. As a consequence, our proposed transcritical hybrid model seems appropriate to capture the essence of its firing mechanisms.

Transcritical Hybrid Modeling of a Thalamocortical Relay Neuron
The phase portrait analysis in the previous section suggests the relevance of a reduced transcritical hybrid model to model TC relay neurons. 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 reproduce and explain the behavior of TC neurons and, in particular, the role of calcium currents. Fig. 9 compares the experimental step response of a TC neuron in vitro and the simulated step response of the transcritical hybrid model (1), both in the low and high calcium conductance modes. As discussed above, the small calcium conductance mode is obtained by choosing a positive w 0 (T-type calcium channels are inactivated), whereas the large calcium conductance mode is obtained by choosing a negative w 0 (T-type calcium channels are deinactivated), all the other parameters being identical in the two modes. An additional variable z accounts for the slow adaptation mechanisms, such as e.g. T-type calcium channels inactivation and variations of intracellular calcium concentration. The hybrid model reproduces the experimental observation: in the lowcalcium 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. Fig. 10 shows the phase-portraits of the transcritical hybrid model in the two modes. Note the great similarity with these phase portraits and the ones of the reduced conductance-based TC neuron model. When w 0 w0 (Fig. 10, left), 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. 10 A, left). When the hyperpolarized state belongs to the upper branch of the vnullcline, no plateau oscillations are possible (Fig. 10 A, left). Furthermore, the relaxation toward the hyperpolarized state is necessarily monotone (i.e. no ADP), as stressed in Fig. 10 B, left. At the generation of the first spike, z increases, which reduces the excitability of the cell (calcium accumulates in the intracellular space for instance, which activates hyperpolarizing calcium pumps), and the v-nullcline upper branch falls. As the neuron remains polarized, z decreases (the intracellular calcium is expelled), 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 (Fig. 10 B, left).
When w 0 v0 (Fig. 10, right), the hyperpolarized state belongs to the lower branch of the v-nullcline, and is more hyperpolarized than for w 0 w0, as in experiments. When a depolarizing current step is applied, this branch falls below the w-nullcline (Fig. 10 A,  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, as in the reduced TC model. These high frequency plateau oscillations (burst) continue until z is sufficiently large (T-type calcium channels are inactivated and calcium accumulates in the cytoplasm). Plateau oscillations then terminate in a (hybrid) saddle-homoclinic bifurcation (Fig. 10 A, 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 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. 10 B, right). A finer and more physiological modeling of the intracellular calcium dynamics could reproduce in vitro recordings with a higher degree of fidelity. Note that the ADP trajectories are slightly different in the reduced conductance based and the transcritical hybrid models. This minor difference is due to a lower time-scale separation in the reduced TC model for low values of V , but the generating mechanisms are similar. We discuss the impact of this difference in the next section.
In order to further verify the physiological consistence of the transcritical model, we compare its behavior with the simulated step response of another quantitative one compartment model    [29,30], a quantitative 200-compartments model [46] (simulations were run in the Neuron environment, based on the configuration files freely available at http://cns.iaf.cnrs-gif.fr/alain_demos.html) and a fold hybrid model [12] of a TC relay cell in the large conductance mode (Fig. S3). For the quantitative models, we plot the trajectory projection on the V {m T plane, where V and m T denotes the somatic membrane potential and the activation gating variable of the somatic T-type calcium current, respectively.
There is a striking similarity between the projection of the trajectories between both quantitative models 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 nullclines _ V V~0 and _ w w~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. Indeed, as highlighted above, reproducing this behavior with the standard reduced HH model necessitates a non physiological alteration of the reset rule ( Fig. S3 [12,49]). This underlines the relevance of the revisited model to capture the richness of neuronal excitability.

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. Fig. 11 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 12. In the fold model, ADPs are generated when trajectories cross the vnullcline from below (Fig. 12, see also Fig. S3). 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 (nonphysiological) spike.
Conversely, ADP generation in the transcritical hybrid model is robustly governed by the attractor S e 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 [50] ensures, as required in biologically meaningful conditions, the robustness of the ADP generation mechanism to small inputs.
The ADP apex is the most excitable part of the trajectory. A physiologically relevant external stimulation can easily generate a spike during this small period, whereas the neuron is barely excitable in the preceding and following periods. This suggests a major role for this ADP in the modulation of external inputs by neuron endogenous rhythm. This excitability can be finely tuned through variations of channel kinetics. For instance, the presence of a ''buckle'' in the ADP trajectory of the reduced conductance based model in Fig. 8 (compare to Fig S3) is an artifact that illustrates the sensitivity trajectories around the ADP apex, a mathematical illustration of the neuron excitability at this particular instant.

Calcium Channels Physiologically Unmask the Physiological Relevance of a Global View of the Reduced Hodgkin-Huxley Phase Portrait
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. It results in the presence of a new voltage nullcline branch in the phase portrait, which lies below the classical inverted N-shaped one. When the calcium channel density is high, neuronal excitability is governed by this lower branch, which accounts for the physiological signature of these currents: spike latency, plateau oscillations and afterdepolarization potential.
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. Indeed, the transcritical singularity is the core mechanism of the classical Hodgkin-Huxley model as well, and the well known fold bifurcation derives from it. 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. This complete picture and its dynamical consequences are brought to live by any slowly activating depolarizing current, calcium currents being the most representative. Note that it may be generated by a slowly activating persistent sodium current as well.

The Proposed Planar Model Differs from Earlier Planar Models that Include Calcium Channels
The proposed planar model is distinctively different from earlier reduced models that include calcium channels in that its slow variable (n or w) aggregates in the same time scale the antagonistic activation of potassium and calcium channels, thereby capturing the non monotonicity of the total ionic current (Fig. 2, right). This property is lost when the activation of calcium channels is treated as a fast variable, that is, set at steady-state in the reduction such as, for instance, in the popular Morris-Lecar model [51]. A classical reference such as [17] nevertheless suggests that the time constant of calcium and potassium activations are comparable, motivating the reduction adopted in the present paper.

The Richness of Neuronal Excitability is Captured in a Two Dimensional Transcritical Hybrid Model
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. Indeed, whereas tonic spiking of these cells is well captured by classical models based on the fold normal form, the generation of burst firing needs non physiological alterations of the phase portrait and the reset rule of these models. On the other hand, both firing patterns can be generated in the transcritical hybrid model through a change in one parameter, which reflects the proportion of calcium channels which are not inactivated.
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.

A Transcritical Hybrid Model as the Basic Unit of Large Population Studies
Detailed modeling of large scale neuronal networks is an invaluable tool in the analysis of the mechanisms underlying the collective behavior of the brain. The recent paper [15] illustrates  that reduced models such as those discussed here are simple enough to be used in large-scale population computations including tens of thousands of neurons. Because the proposed transcritical hybrid model captures a broader class of neuronal dynamics and is particularly adapted to neurons expressing a high density of calcium channels, it will be an ideal candidate for physiologically realistic studies of high-dimensional neuronal networks such as the thalamocortical circuitry or the basal ganglia motor loop.

Equation and Parameters of the Complete Model
The augmented HH model reads.
For the HH dynamics, we use the parameters of the original paper [1]. As all other HH currents, the additional calcium current obeys Ohm's law.
where g g Ca is the maximum calcium conductance, V Ca is the calcium Nernst potential, and d is the calcium activation gating variable. In our study, we consider the kinetics of L-type calcium channels (adapted from [3]): L-type calcium channels are calcium permeable channels that activate at high-threshold and are found in many cell types, such as dopaminergic or serotoninergic neurons, muscle and cardiac cells. The functions a x ,b x , x~n,m,h, can be found in the paper [1]. The value for the potassium Nernst potential V K~{ 12 is the same as in [1], while the sodium Nersnt and the leak Nernst potential are rounded to V Na~1 20 and V l~1 0:6, respectively. The values of the sodium g g Na~1 20, potassium g g K~3 6, and leak g l (maximum) conductances are the same as in [1]. The calcium Nernst potential is given by V Ca~1 50. The numerical simulations of Fig. 1(right) are obtained by taking g Ca~0 :4 and I pump~{ 17.

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 )~a m (V )=(a m (V )zb m (V )); ii) exploiting the approximate linear relation, originally proposed in [7], hzsn~c, with s,c^1. iii) exploiting the correlation between the potassium and calcium gating kinetics to approximate the calcium activation gating variable d as a static function of the potassium activation gating variable n (the simple relationship d :~n 3 provides a satisfactory fit, see Fig. S4). Applying this reduction to (1) with parameters as above, we obtain the planar system (g Ca~2 :7).
and similarly for other variables.

Reduced Model of TC Neurons
The complete TC neuron model reads (adapted from [46]) The dynamics of the ionic currents are similar to [46]. The conductance values are taken as follows (in mS): g g K,DR~1 , g g Na~6 0, g Ca,T~0 :1, g Na,leak~2 :5 10 {3 and g K,leak~7 10 {3 . The membrane capacitance C m is set to 290pF . Holding currents are {0:3 nA for bursting mode, and 0:2 nA for tonic mode. The amplitude of the current steps are 0:3 nA and 0:8 nA in bursting and tonic modes, respectively.
The reduction of the model is performed by assuming an instantaneous sodium activation m:m ? (V ) and merging h, n and m T so that they have similar time courses in the range of interest. The reduced model reads.
For the phase plane analysis, inactivation of T-type calcium channels h T is considered as a parameter, its timescale being much slower than the other variables. To ensure a similar resting potential in both models, the holding current of the reduced model is adjusted to 0:18nA in tonic mode. The reduced model finely reproduces the behavior of the complete model in our range of interest (Fig. S5).

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: _ v v~v 2 zbvw{w 2 zI{z i f v §v th , then _ w w~e(av{wzw 0 ) v/c, w/d, _ z z~{e z z, z/zzd z : with a~0:1,b~{3,c~15,d~15,e~1,v th~1 00,e z~0 :1,d z~4 0.
The low-calcium mode correspond to w 0~3 :2, the high-calcium mode to w 0~{ 4. The injected current I~{5 when the stimulation is off and I~85 when the stimulation is on. We stress two properties that are not captured by the reduced model (1): Firstly, the slow adaptation variable (i.e. z) dynamics are generally coupled with the membrane voltage also in the subthreshold phase. Secondly, depending on the type of the modeled calcium current, the calcium conductance, here reflected by w 0 (Section), might slowly change in response to voltage and intracellular calcium variations. Both properties are neglected in (1).

Numerical Simulations
Numerical simulations were run with MATLAB. (http://www. mathworks.com), apart from the 200-compartment model, which was simulated with the NEURON software environment. (http:// www.neuron.yale.edu). The numerical bifurcation analysis has been obtained with XPP environment. (http://www.math.pitt. edu/~bard/xpp/xpp.html).  [29,30], a 200compartment model of TC neurons [46] and the fold hybrid model [12] in large calcium conductance mode. In the phase-portrait of the hybrid models, vand w-nullclines are drawn as full and dashed lines, respectively, and trajectories are drawn as red oriented lines. The black full line represents the vnullcline at the onset of the stimulation. The gray full line represents the v-nullcline at the end of the burst. The phaseportrait of the conductance based models depicts the trajectory projection on the V {m T plane, where V and m T denotes the somatic membrane potential and the activation gating variable of the somatic T-type calcium current, respectively. (EPS)