A model of on/off transitions in neurons of the deep cerebellar nuclei: deciphering the underlying ionic mechanisms

The neurons of the deep cerebellar nuclei (DCNn) represent the main functional link between the cerebellar cortex and the rest of the central nervous system. Therefore, understanding the electrophysiological properties of DCNn is of fundamental importance to understand the overall functioning of the cerebellum. Experimental data suggest that DCNn can reversibly switch between two states: the firing of spikes (F state) and a stable depolarized state (SD state). We introduce a new biophysical model of the DCNn membrane electro-responsiveness to investigate how the interplay between the documented conductances identified in DCNn give rise to these states. In the model, the F state emerges as an isola of limit cycles, i.e. a closed loop of periodic solutions disconnected from the branch of SD fixed points. This bifurcation structure endows the model with the ability to reproduce the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\text{F}\to \text{SD}$\end{document}F→SD transition triggered by hyperpolarizing current pulses. The model also reproduces the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\text{F}\to \text{SD}$\end{document}F→SD transition induced by blocking Ca currents and ascribes this transition to the blocking of the high-threshold Ca current. The model suggests that intracellular current injections can trigger fully reversible \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\text{F}\leftrightarrow \text{SD}$\end{document}F↔SD transitions. Investigation of low-dimension reduced models suggests that the voltage-dependent Na current is prominent for these dynamical features. Finally, simulations of the model suggest that physiological synaptic inputs may trigger \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\text{F}\leftrightarrow \text{SD}$\end{document}F↔SD transitions. These transitions could explain the puzzling observation of positively correlated activities of connected Purkinje cells and DCNn despite the former inhibit the latter.


Introduction
The connectivity of the cerebellum places deep cerebellar nuclei neurons (DCNn) in a strategic position. Their axons are the main output of the cerebellum, projecting to the forebrain, the brain stem and the spinal cord. In turn, DCNn are the main target of the Purkinje cell (PC) axons, which are themselves the sole output system of the cerebellar cortex. As noted by Lllinás and Muhlethaler [1] the cerebella nuclear cells are therefore the main functional link between the cerebellar cortex and the rest of the central nervous system. It follows that a detailed understanding of electrophysiological properties of the DCNn is one of the fundamental prerequisites to understanding the overall functioning It is therefore tempting to conjecture that, in systems with dimension >2, the separatrix is also the global stable manifold of the uLC, that is a manifold with dimension k + 1 with k being the number of Floquet exponents <1. Such a manifold may have a complicated shape. But if solutions are bounded, the simplest manifold achieving separation of the attraction basins of the sFP and sLC would be a closed manifold (i.e. compact and without bounds). Characterizing the separatrix may have remained a theoretical issue since the coexistence of the sFP and the sLC has, to our knowledge, only be reported in squid axons bathed in non-physiological low Ca salines. However, the observation of such a coexistence in DCNn put this mathematical question into a neurophysiological context since characteristics of the separatrix obviously constrain the types of stimuli capable to trigger sFP ↔ sLC state changes.
The repertoire of membrane conductances expressed by DCNn has been the subject of several studies. These conductances comprise two calcium carried currents [14], a noninactivating high-threshold Ca current (I CaH ) and an inactivating (transient) one (I CaT ) and a calcium-dependent potassium current (I KCa , [15]). In addition, Raman et al. [4] have identified a tonic non-cationic current, I TCN . However, in spite of these characterizations, we still lack an understanding of how this set of conductances actually implements the intrinsic properties of the DCNn. In the absence of a theoretical framework, it is difficult to relate the response of the DCNn membrane voltage to its active conductances. In particular, several of the available experimental observations raise unsolved issues. For instance, transitions from the spontaneous firing state to the silent depolarized one can be triggered by inhibition of the calcium currents I CaH and I CaT [4] but it is not clear whether the inhibition must concern both currents or only one of them. Likewise, it is not known if the transitions observed in vitro using current pulse injection are also triggered by synaptic currents.
Biophysical models have proven highly useful to understand how active conductances interact to underpin active electrical signals in neurons. For instance, the acclaimed Hodgkin and Huxley model of the action potential in the squid axon has proved a very powerful approach to our understanding of the generation of action potentials [16]. For DCNn, the only detailed (compartmental) biophysical model of DCNn thus far (to our knowledge) has been proposed in a pioneer study by Steuber et al. [7]. The more recent model of Sudhakar et al. [17] actually adopts the structure and equations of the Steuber et al. 's model and only uses different values for the densities of active conductances. However, neither of these articles evidences capabilities of this model to reproduce state transitions experimentally observed by Jahnsen [2] and Raman et al. [4]. Moreover, the Steuber et al. 's model relies on several assumptions which experimental support can be questioned. Thus, it assumes that DCNn express two different Hodgkin-Huxley-type voltage-dependent Na currents, a fast one and a slowly inactivating one (i.e. persistent). Expression of a persistent Na current by DCNn was postulated by Jahnsen [2] and Lllinás and Muhlethaler [1] and its existence has since been confirmed by two experimental studies (Raman et al. [4] and Afshari et al. [18]). However, these studies suggest that the persistent Na current belongs to the class of so-called resurgent Na currents (Raman and Bean, [19]). Resurgent currents involve a specific mode of channel aperture whose faithful modeling requires a set of 12 ODEs (Raman and Bean, [19]). Moreover, the Steuber et al. 's model assumes a non-uniform distribution of active conductances between the dendrites and soma of DCNn which remained to be proved (with the exception of T Ca channels that cluster to the soma and proximal dendrites of DCNn; see McKay et al. [20]). Here we propose a new biophysical model of DCNn electrogenesis. To allow mathematical analysis of the model, we focus on an isopotential (non-compartmental) description, in which the DCNn membrane potential evolves under the influence of the 6 conductances mentioned above (I NaV , I Kdr , I CaH , I CaT , I KCa , and I TCN ). Spontaneous firing of fast spikes emerges in the model as an isola of limit cycles in a range of driving currents. Within this range, the model exhibits another coexisting stable attractor: a stable depolarized stationary-state. Owing to this peculiar bifurcation structure, the model reproduces the transitions between spontaneous firing and a silent state that have been reported experimentally [2,4]. Like in these experimental reports, the transitions in the model can be triggered by hyperpolarizing current pulses or by blocking Ca currents in the absence of hyperpolarization. Analysis of the model suggests that the latter transition is specifically due to the inhibition of the I CaH current. Finally, simulations of the model suggest that physiological synaptic inputs may trigger the state transitions.

Methods
It is currently unknown whether DCNn are electrically compact neurons nor whether their ion channels have non-uniform distributions over their membrane surface (with the exception of high threshold voltage-dependent Ca channels [14]). In the absence of evidence for a functional role of these heterogeneities, we chose to build an isopotential model of DCNn to (i) investigate the hypothesis that DCNn have two coexisting stable states of activity, (ii) identify the interactions between the active membrane ion currents that are responsible for this electric feature and (iii) investigate whether physiological inputs could trigger transitions between these states. DCNn do not form a homogeneous population from the phenotypic standpoint. Their population has been divided into two classes (see [21] for an introduction) which partially correlate with the location of their target neurons. Glutamatergic DCNn project excitatory inputs to motor centers in the brain stem, the mesencephalon, the thalamus and pre-cerebellar nuclei though the mossy fiber pathway. GABAergic DCNn achieve inhibitory synapses on neurons of the inferior olive [22] and glutamatergic DCNn [23]. A class of rare small-sized DCNn has also been identified [24]. In addition to their neurotransmitters and projection sites, these classes of neurons can also express different sets of ion channels [25]. However, it is unknown to what extent their intrinsic electrophysiological properties are heterogeneous. In the absence of this information, we propose below a generic model for the electrogenesis of DCNn.

State variables and equations of the model
Our standard model has 6 state variables. The dynamics of these variables is governed by 6 associated ordinary differential equations (ODE) that we describe below. The two main state variables of the model are the membrane potential, V (mV), and the cytoplasmic concentration of free Ca 2+ ions [Ca] (μM). Their respective ODE read and Equation (1) is standard for biophysical models of neurons (see e.g. [26]). It can be derived from the first Kirchhoff law of electromagnetism which states that electric charges can neither be created nor destroyed. The set of membrane ion currents I ion i (nAcm -2 ) appearing in Eq. (1) is detailed below. For the input current I S we distinguish tonic (I DC ) and phasic (I ϕ ) components and write I S = I DC + I ϕ . The phasic component was written as the product of an ohmic term g ϕ (V -E ϕ ) and the product of two Heaviside step functions of time H(tt start )H(t endt) = 1 for start < t < end, 0 elsewhere.
The ODE for the cytoplasmic free calcium ions concentration (Eq. (2)) is a balance equation. The source term corresponds to calcium ions entry into the cytoplasm through membrane channels underlying the I CaH and I CaT currents. The sink term corresponds to the extrusion of cytoplasmic Ca 2+ ions by cytoplasmic membrane pumps and their internalization by Ca pumps in the endoplasmic reticulum membrane, both processes being modeled with a simple linear term. The ODE actually includes a second sink term corresponding to Ca 2+ ions buffering by Ca-binding proteins. It appears as the denominator of the right-hand side of Eq. (2) owing to the hypothesis that the binding of Ca 2+ ions is very fast (see [27] for the derivation of Eq. (2) and Table 1 for parameter values).
The ion membrane conductances (for currents obeying the Nernst equation) and permeabilities (for currents obeying the Goldman constant-field equation) were modeled as the product of a maximum conductance g (μS/cm 2 ) or permeability P (cm/s) and either the product of voltage-dependent activation and inactivation variables m and h or a Cadependent variable w, i.e. g(P) = g(P)m p h q , in which p and q are integers. The remaining four states variables of the model correspond to in(activation) variables of ion currents I NaV , I Kdr and I CaT and namely are h NaV , m Kdr , m CaT and h CaT . Their dynamics obey first order differential equations of the form [16] in which x ∞ and τ x stand for the steady-state value and the exponential time constant of variable x, respectively. The available experimental data on the voltage-dependence in DCNn [4,18] were insufficient to constrain a fully detailed Hodgkin-Huxley-type model of these currents since they do not document the voltage dependence of rate functions governing state transitions of ion channels. We therefore adopted the modeling approach of Hughenard and McCormick [29] in which x ∞ and τ x are considered as independent functions of V (however, see the discussion) of the form Equation (5) is a so-called Boltzmann function in which the '-' and '+' signs hold, respectively, for the inactivation (h NaV , h CaT ) and the activation (m Kdr , m CaT ) variables. V x and k x are referred to as half-activation potentials and activation slopes of the variable x, respectively. Equation (6) is standard for reproducing the bell-shape of the activation time constant of voltage-dependent currents [29]. The parameter τ x 0 accounts for the observation that lim V →+∞ τ x > 0 for the inactivation variable of some particular currents (e.g. I Na in the squid giant axon [11]. The α τ x parameter allows one to introduce asymmetry in the shape of the time constant as observed for some currents (see e.g. [26], p. 48). We give below a detailed account of the model equation for each of the currents. I NaV : As in most neurons studied so far, the spikes of the DCNn are blocked by tetrodotoxin (TTX) [4] suggesting that the depolarizing phase of these spikes is underlain by a Hodgkin-Huxley-type of voltage-dependent Na current, I NaV . For this reason, we used the classical Hodgkin-Huxley formalism to model I NaV in DCNn. However, the time constant of m NaV in neurons is usually much smaller than that of h NaV . We therefore used the classical quasi-steady-state approximation replacing m NaV by its steady-state value m NaV∞ (see e.g. [31]). Llinás and Muhlethaler [1] initially suggested that I NaV in DCNn could have a persistent component. Afshari et al. [18] have confirmed that I NaV in DCNn has a persistent component that they attribute to the mechanism of resurgent current initially identified in cerebellar Purkinje cells [19]. Given that the resurgent component only amounts to <4% of the total I NaV in DCNn [18] and that a faithful modeling of I NaV endowed with resurgence properties requires a set of 12 ODE [19], we chose to not introduce (see, however, the Discussion) a persistent I NaV in our model and we rather used a classical Hodgkin-Huxley I NaV where m NaV∞ is given by Eq. (5) and h NaV by Eqs. (4)- (6). I Kdr : voltage-dependent K currents of the Kv3-type exhibit fast activation allowing ∼1 ms-duration Na spikes in several neuron types [32]. Lllinás and Muhlethaler [1] report spike duration (measured at half-width) ranging from 0.44 to 0.7 ms in DCNn suggesting that spike repolarization in these neurons is achieved by Kv3 channels. Raman et al. [4] have dissected out a high-threshold voltage-dependent potassium current in DCNn that is be presumably involved in the repolarization of their spikes. The activation time constant of this current (12 ms at +12 mV) is too large for Kv3 channels and rather points to the Kv2 family of K channels [33]. However, Raman et al. [4] only used 1 mM TEA whereas to block voltage-dependent K channels whereas blocking Kv2 channels require several millimolar concentration of TEA to block [34]. Moreover, no molecular study has yet, to our knowledge, demonstrated expression of Kv2 channels in DCNn. On the opposite, several studies have shown that DCNn express all four Kv3 subunits (see e.g. [35] and [32]) and the electrophysiological study of Lamont [28] suggests that Kv3 channels in DCNn are functional. DCNn can fire spikes at frequencies up to >100 Hz which seem unattainable with Kv2 serving to repolarize spikes. Owing to these findings, our model assumes that I Kdr is carried by fast-activating K channels (however, see the discussion for the involvement of putative Kv2 channels), where m Kdr is given by Eqs. (4)- (6). The parameters from this current were taken from a model of Purkinje cell dendrite [27]. I TCN : There is experimental evidence for a TTX-insensitive current in DCNn that can depolarize them beyond their spike threshold [4], this current being tonic (i.e. Vindependent) and non-selective for cations [4]. We modeled this current using Ohm's law with a constant conductance, i.e. as an effective leak current: Calcium currents I CaH and I CaT : Owing to the large gradient of Ca 2+ ions across the cytoplasmic membrane of neurons (typically [Ca] o /[Ca] i 2 × 10 5 in a neuron at rest), voltage-dependent Ca currents are not adequately modeled by the Nernst equation. This equation accurately describes ion fluxes only close to thermodynamic equilibrium (see e.g. [26]). Far from equilibrium, Ca currents exhibit a rectification in their I/V relation which is well accounted for by the constant-field equation of Goldman (see e.g. Chap. 13 in [26]). We therefore used this equation to describe Ca currents in our model. As mentioned above, two types of Ca currents have been identified in DCNn [14]: a non-inactivating high-threshold current I CaH and an inactivating current I CaT . Both currents are given by where Ca x ∈ {CaT, CaH} and α(V ) = e -2FV RT . A temperature T = 298 K was used in all simulations. The molecular identity of ion channels underlying I CaH in DCNn is currently not precisely known. Gauck et al. [14] have shown that a large fraction of I CaH in DCNn is non-inactivating and is therefore likely composed of L-, T-or R-type currents. They have also observed that a small fraction of I CaH inactivates, thereby also indicating the presence of N-type calcium channels. For the sake of simplicity and since L-type calcium channels open as quickly as N-types channels in neurons [36], we derived a simple equation for P CaH by assuming that it activates instantaneously and does not inactivate, i.e.
Setting an equation for I CaT proved more challenging. Indeed, this current can be underlain by three isoforms of calcium channels, Cav3.1, Cav3.1 and Cav3.3 which have marked different time constants for both activation and inactivation [18]. Molineux et al. [25] have demonstrated the expression of Cav3.1 in both GABAergic and non-GABAergic DCNn, whereas Cav3.3 expression is restricted to non-GABAergic DCNn. Since both types of DCNn exhibit rebound properties, these findings suggest that the expression of Cav3.1 is sufficient to generate a rebound in both cell types. McKay et al. [20] have shown that Cav3.1 channels are largely confined to the soma of the DCNn. McRory et al. [30] have found that both τ act and τ inact of Cav3.1-3 channels decay monotonically with voltage toward a voltage-independent nonzero minimum. This feature is well reproduced by the thalamic I CaT model of Destexhe et al. [37] which we adopted here for the DCNn. Its permeability equation reads where m CaT and h CaT are given by Eqs. (4)- (6). Nevertheless, we had to adapt the original model of Destexhe et al. [37] to provide a faithful account of I CaT in DCNn. Firstly, we adapted the values of the time constants τ m CaT and τ h CaT of the original model to match the properties of the Cav3.1 channels expressed by DCNn [25]. Secondly, the model by Destexhe et al. uses a discontinuous piecewise function to describe τ h CaT . This prevents guarantying existence and uniqueness of the solutions of the model according to the Cauchy-Lipschitz theorem (see e.g. [38]). We therefore replaced the original formulation of τ h CaT in Destexhe et al. by the continuous and differentiable formulation of Eq. (6). I KCa : The DCNn express calcium-dependent potassium currents of the SK (smallconductance) type [25]. Those potassium-specific channels are gated by calcium according to a Hill-type equation (see e.g. [39]) which assumes instantaneous gating of SK channels by calcium: Finally, the model includes a leakage current. The value of its conductance g Leak was set to a value allowing the model to reproduce the passive time constant of DCNn recorded in vitro (as probed by small hyperpolarizing pulses)

Boundedness of solutions
We follow the approach of Cronin [40] to prove that all physiologically significant solutions of the model are bounded. Let us first consider the dynamics of state variables h NaV , m Kdr , m CaT  Moreover, for all physiologically significant initial conditions on V we have 0 < x < 1 according to Eq. (5). It follows that dx dt > 0 if x = 0 and dx dt < 0 if x = 1. This implies that all physiological solutions remain in the Let us now turn to Eq. (1) in the particular case where P CaH = P CaT = 0 (namely the model deprived of its Ca currents) and let E max = max{E Na , E K , E Leak } and E min = min{E Na , E K , E Leak }. Recall then that currents I Na , E K and E Leak all have the form Proof of the boundedness of solutions is finalized by extending the above result to the full model (i.e. including Ca currents) as follows.

Theoretical evidence for two stable states of electric activity in DCNn
Consistent experimental data show that DCNn are pacemaker neurons in vitro [2][3][4][5][6][7]. In the absence of a depolarizing bias current (I DC = 0), our neuron model successfully repro- duces this property, with the spontaneous emission of self-sustained fast (∼1 ms duration) spikes ( Fig. 1A 1 ). We refer to this mode of activity as the firing (F) mode. The spontaneous firing frequency in the model is 28.9 Hz (versus 26 Hz in [2] and ∼20 Hz in [4]). The amplitude of spontaneous spikes is 67 mV (in the range of extreme values reported experimentally: 57 ± 5 mV in [1] and 82 mV in [4]). As shown in Fig. 1A 1 , a 300 ms hyperpolarizing pulse of -2 μAcm -2 amplitude interrupts spike firing in the model. Firing resumes at the end of the pulse with an initial phase of increased firing frequency of ∼200 ms duration ( Fig. 1A 1 and C 1 ). The concentration of free cytosolic calcium [Ca] first decreases during the pulse and quickly recovers its baseline level after the end of the pulse (Fig. 1B 1 ).
Increasing the pulse amplitude shows non-monotonous effects. A moderate increase of the amplitude of the pulse (to e.g. -2.5 μAcm -2 , Fig. 1A 2 -C 2 ) yields the same behavior, but with a larger transient firing frequency increase at pulse break. The peak frequency during this rebound is close to twice the basal frequency ( Fig. 1C 2 ), in close agreement with the experimental data of Raman et al. [4]. We also note the appearance of a rebound [Ca] increase at the pulse break, which is consistent with the observations of Raman et al. [4] (Fig. 1B 2 ). Large pulse amplitudes, however, produce a dramatic change of response (Fig. 1A 3 -C 3 ). For instance, a -3.9 μAcm -2 pulse triggers a switch of the model to a silent depolarized (SD) state in which the membrane potential settles at -38 mV (Fig. 1A 3 ). This model feature reproduces the observation of Jahnsen ([2]; see Fig. 5C) that a hyperpolarizing current pulse can switch spontaneously firing DCNn to a silent mode of activity. In the model, the SD state is stable as evidenced by the fact that, once settled in this state, the model remains in this state despite small depolarizing or hyperpolarizing pulses of current ( Fig. 1A 3 ). However, this stability is only local since larger-amplitude hyperpolarizing current pulses (Fig. 1D 1 ) readily switch the model back to its firing mode. This new behavior is also consistent with the following experimental result by Jahnsen ([2]; see Fig. 5E): when DCNn are switched to the SD state, a subsequent hyperpolarizing pulse can switch them back to the F mode. However, our model extends this finding by showing that depolarizing current should also be able to trigger this transition ( Fig. 1D 2 ). Therefore, our model reproduces the experimental observations that hyperpolarizing pulses can switch the DCNn from their spontaneous low frequency firing (F) state to a stable depolarized (SD) state. The SD state is only locally stable and injection of large hyperpolarizing or depolarizing current pulses switches the DCNn model back to its F state. In the following we provide a mathematical analysis of the model in order to dissect the contribution of the different membrane ion currents to these dynamical features.

Bifurcation analysis of the model and the silent depolarized state (SD)
Bifurcation analysis provides an explanatory scheme for the emergence of the overall properties of a dynamical system (see e.g. [41] for an introduction). In bifurcation analysis of neuron models, the steady injected current (I DC in our model) is the most important bifurcation parameter as it is used in experiments as a model of synaptic inputs to investigate the electric properties of the neuron. The power of this approach was first demonstrated by Rinzel and Ermentrout [41] who showed that the type I and II responses identified by Hodgkin [42] in crab neurons are due to different bifurcations (see Izhikevich [43] for a comprehensive exposure of the geometry of excitability of nerve cells). We therefore followed this approach and we built bifurcation diagrams of the model with I DC as the bifurcation parameter.
In the presence of hyperpolarizing bias currents, the model has a unique stable attractor, the SD state ( Fig. 2A). This state is actually stable in the global sense (i.e. with respect to any perturbation whatever its amplitude) since it is the only attractor in the phase space. Firing arises upon increasing the bias current I DC to I F 1 = -255.3 nAcm -2 , a value at which a unique neutrally stable periodic orbit appears ( Fig. 2A). Two different branches of Limit Cycles (LC) emerge from this orbit: one of them (green) corresponds to stable LC (sLC) whereas the other one corresponds to unstable LC (uLC). This type of bifurcation is usually referred to as fold bifurcation of LC. The neutrally stable LC has a frequency of 17.8 Hz at the F 1 bifurcation point (Fig. 2B). These results suggest that the DCNn are type II neurons characterized by both a nonzero frequency and a nonzero spike amplitude at the bifurcation point [38]. Upon increasing I DC , the amplitude of the sLC steeply increases while that of the uLC first decreases, exhibiting a minimum at I DC 1000 nAcm -2 (inset Fig. 2A), and increases afterwards. Likewise, the amplitude of the sLC shows a maximum at I DC 1000 nAcm -2 then decreases beyond this value. The parallel decrease of the amplitude of sLC and increase of the amplitude of uLC ultimately results in the merging of the two LC branches into another neutrally stable LC through a second fold bifurcation at I F 2 = 4256 nAcm -2 ( Fig. 2A). The theoretical maximum firing frequency of DCNn occurs at I F 2 and is 108 Hz, in the range of experimental maximum frequencies (80-120 Hz [7]). The branches of sLC and uLC thereby form a closed loop of limit cycles solutions of the model which is usually referred to as an isola of limit cycles (see e.g. Avitabile et al. [44] and Labouriau [45,46] in the case of the HH model). This isola coexists with a branch of stable FP. Figure 2A thereby suggests that DCNn could switch between silent and firing modes of activity within the range of current [I F 1 , I F 2 ] in response to proper inputs. Spontaneous firing in the model disappears upon zeroing g NaV (Fig. 2C). This result reproduces the finding that I Na underlies the rising phase of the spike in the DCNn [4]. After g NaV is zeroed, the membrane potential settles to a stable state with membrane potential V = -44.4 mV, i.e. above the spike undershoot (V U = -58.1 mV). This result is at odds with properties of all neuron types studied so far (to the best of our knowledge) of which the membrane potential converges to values below V U after blocking of the Na currents with TTX. Nevertheless, this model feature closely reproduces the experimental observation that, after blocking I NaV with TTX, the DCNn settle to a depolarized stable voltage of -42 ± 2 mV [4]. This property is well explained by the finding that I TCN provides a depolarizing current to the membrane, as illustrated in Fig. 2D which shows that the branch of FP is shifted downright when I TCN is blocked.
The report of a bistability between an isola of LC and a stable FP branch in a neuron membrane model is not unprecedented. Guttman et al. [12] have described similar dynamics in a model of the squid axon, which exhibits a stable point attractor coexisting with a stable firing mode of spikes when bathed in low [Ca] 0 solutions. This prompted us to investigate the hypothesis that [Ca] 0 may control bistability in DCNn by studying how [Ca] 0 affects the bifurcation diagram of the model (Fig. 3). With 0.5 mM external Ca 2+ (Fig. 3A), the isola is characterized by the coexistence of a stable fixed point (red), an unstable limit  Fig. 2, in particular for the envelope of the sLC (green), the uLC (blue) and the stable branch of fixed points FP (red). Around 1.3 mM extracellular Ca (panel D), the uLC collides with FP thus creating two Hopf bifurcation points (H 1 and H 2 ). Bifurcation F 2 is also altered by the appearance of two additional folds of limit cycles (F 3 and F 4 ). F 1 is also eventually modified by the appearance of another fold of limit cycles (F 5 ). Panels in (G) show the corresponding two-parameter bifurcation diagrams: G 1 is the overall diagram while G 1 and G 2 show magnifications of G 1 around [Ca] 0 = 1.3 and 2.0 mM, respectively. The curves show how the bifurcation points described in A-F change when [Ca] 0 or I DC are varied. The branch of folds F 1 and F 2 are shown in pink and the two branches corresponding to the Hopf bifurcations H 1 and H 2 are shown in cyan. The additional folds associated to F 1 (F 3 and F 4 ) are shown in purple and the additional fold that is associated to F 2 (F 5 ), is in black. Those Fold branches appear and disappear by collisions with each other at codimension-2 bifurcation points called cusps of limit cycles, show as points C. Thus, F 3 and F 4 collides at cusp C 1 , F 3 and F 1 at C 2 , and F 2 and F 5 at C 5 cycle (min-max values in blue) and a stable limit cycle (min-max values in green) within the range I DC = [0, 4] μA/cm 2 . Two folds of limit cycles (F 1 and F 2 , pink in Fig. 3G) locate the coalescence of the two limit cycles. With increasing external Ca 2+ , the unstable limit cycle gets closer to the fixed point (Fig. 3B, C), colliding it for [Ca] 0 1.26 mM (Fig. 3C). The collision with the fixed point gives birth to (Fig. 3G 1 ) two subcritical Hopf bifurcations (H 1 and H 2 , blue in G) between which the fixed point becomes unstable (black in D). In addition, a pair of limit cycles (one stable, one unstable) appears by a cusp of LC (C 1 in Fig. 3G 2 ), thus giving rise to two new folds of LC (F 3 and F 4 , purple in Fig. 3G 2 ). F 4 disappears by coalescence with H 1 and F 3 with F 1 in the cusp of limit cycles C 2 (Fig. 3G 2 ), leaving H 1 as the only bifurcation point in the zone for [Ca] 0 > 1.4 mM (Fig. 3E). At even larger [Ca] 0 (around 2 mM), a new stable LC appears close to F 2 and H 2 as a result of H 2 becoming supercritical, giving rise to F 5 (Fig. 3F, gray in Fig. 3G 3 ). F 2 and F 5 in turn disappear by coalescence at the cusp of limit cycles C 3 (Fig. 3G 3 ), leaving H 2 as the only bifurcation point in the zone for [Ca] 0 > 2.1 mM (Fig. 3G 1 ). These results therefore suggest the existence of a bistability (stable fixed point SD + stable limit cycle sLC) at physiological [Ca] 0 in the electro-responsiveness of DCNn.

Reversible F ↔ DB transitions induced by blocking I CaH : role of the separatrix between the F and SD states
Raman et al. [4] report that blocking Ca currents of DCNn switches their spontaneously firing activity to a silent depolarized mode at -37 mV. Our model reproduces this experimental result as it switches onto a silent depolarized state (SD) after blocking both of its Ca currents, I CaH and I CaT (not illustrated). The membrane voltage adopted by our model (-38 mV) is very close to the experimental value reported by Raman et al. [4]. The experimental protocol of Raman et al. prevented them to determine whether the F → SD switch results from blocking both Ca currents or from blocking only one of them because they blocked both currents with 2 mM cobalt (see [47]). We could readily address this question with the model. Blocking I CaT (i.e. setting P CaT = 0) does not induce a F → SD transition (result not shown). This finding is not surprising since I CaT is almost completely inactivated at the spikes undershoot voltage (h CaT 5.5 × 10 -4 ), which is the more negative V value spontaneously reached by the model in its F mode. On the opposite, blocking I CaH induces a F → SD transition (Fig. 4A). Following this switch, brief pulses of hyper- (Fig. 4A) or depolarizing (Fig. 4B) injected current reset firing in the model. Figure 4C explains this finding by showing that both the branch of FP and the isola of LC of the standard model are preserved after blocking I CaH , thereby allowing theoretical F ↔ SD state transitions. We now give evidence that the mechanism of the F → SD transition induced by blocking I CaH finds its origin in the modifications of the shape of the isola of LC illustrated in Fig. 4C. Recall from Sect. 1 that the SD state is only locally stable. It therefore has a basin of attraction that is the set of initial conditions from which the model converges onto the SD state. However, in the bistable regime the coexisting sLC also has its own basin of attraction. Cooley et al. [8] early conjectured that for membrane models with only two variables, this situation implicates the existence of an uLC separating the basins of attraction of the two stable states. In other words, the separatrix between the basins of attraction of the SD and the F states in two-dimensional models is the orbit of the uLC in the phase plane. In n > 2 dimension models the separatrix must be a manifold of higher dimension than the one-dimensional uLC since the separatrix partitions the phase space into two attraction basins. Guttman et al. [12] have hypothesized that the separatrix contains the uLC and it is reasonable to suppose that it is a n -1 dimension manifold. To our knowledge, neither of these two conjectures have been proved until now. Assuming that they are true (see the appendix), we propose that blocking I CaH in our six-dimensional model triggers the F → DB transition by reducing the maximum amplitude of the uLC (compare blue and yellow undershoot voltages of the uLC in Fig. 4C). This shrinks the attraction basin of the sLC at such point that a trajectory located in the basin of attraction of the sLC before blocking I CaH finds itself on the other side of the separatrix (inside the basin of attraction of the FP) when I CaH is blocked and eventually converge to the SD state. When a small tonic current (I DC = 37 nAcm -2 ) is injected, blocking I CaH is no more able to trigger the F → SD transition although the membrane voltage crosses several times the undershoot voltage of the uLC (Fig. 4D). This confirms that the separatrix between the F and SD attraction basins does not reduce to the uLC and is rather a manifold of larger dimension. We provide further evidence for this conjecture and that of Guttman et al. [12] in the Appendix. Taken together these results suggest that the coexistence of a branch of stable silent states with an isola of LC in the DCNn does not rest on the expression of Ca currents by the DCNn. I KCa cannot play a central role in these dynamical features since it is almost completely inactivated at the low [Ca] level reached by the model when the Ca currents are blocked (m KCa 10 -3 at [Ca] 50 nM). The origin of the peculiar bifurcation diagram of the standard model illustrated in Fig. 2A has therefore to be searched in the interactions between the remaining active currents in the model, namely I Na , I Kdr and I TCN . This question is addressed after the following section, in which we expose the role played by I CaT in the F → SD transitions triggered by hyperpolarizing current pulses in the standard model. Figure 5A 1 (reproducing Fig. 1A 1 above) recalls that large hyperpolarizing current pulses can trigger a F → SD transition in the model. However, Fig. 5A 2 shows that the same cur-  (12)). We first investigated whether long duration current pulses (i.e. of duration τ h CaT ) are able to trigger a F → SD transition in a variant model in which h CaT is no longer a function of time and membrane voltage but a constant value freely adjustable in the [0, 1] continuous interval. With h CaT = 0, the bifurcation diagram of this variant model is nearly identical to that of the standard model (not illustrated). Increasing h CaT to 0.025 induces a loss of stability of the FP branch with the appearance of two saddle-node bifurcations, SN 1 and SN 2 (Fig. 5Ba). The FP branch is stable below I SN 1 Fig. 5B 1 ) delimiting the attraction basins of the two stable FP states. The gap between the two saddle-node bifurcations widens when h CaT is increased (Fig. 5B 2 ). In addition, the F 1 limit of the isola of LC is shifted leftward so that I F 1 ∈ [I SN 1 , I SN 2 ]: the model becomes tristable in this range of injected currents (two stable silent FP states and a stable LC). Both modifications of the bifurcation diagram go on as h CaT is further increased up to 0.045 where the branch of uLC collides the unstable branch of FP points resulting in a homoclinic bifurcation of uLCs (Fig. 5B 3 ). Above this h CaT value, another bifurcation occurs with the homoclinic bifurcation splitting into two homoclinic bifurcations, H m1 and H m2 (see Fig. 5B 4 for h CaT = 0.05). Both bifurcations correspond to the collision of homoclinic orbits (corresponding to uLCs with a frequency equal to zero) with the unstable branch of FP. But as h CaT keeps increasing, the F 1 limit of the isola keeps shifting to the left and the H m1 bifurcation turns into a homoclinic bifurcation of sLC at saddle node (for h CaT = 0.07, Fig. 5B 5 ). In parallel with these effects on the left-hand side of the isola, increasing h CaT also shifts the F 2 bifurcation to the left (see the sequence of Fig. 5B 1-6 ). Above h CaT = 0.06, I F 2 < 0 so that the only attractor remaining in the phase space for I DC = 0 is the SD state (solid red circle in Fig. 5B [6][7] ) and the model is therefore forced to converge onto this SD state. We therefore conclude that hyperpolarizing current pulses trigger the F → DB transition because they de-inactivate I CaT which leads to the disappearance of the isola of LC. In support of this conclusion, Fig. 5C shows that the minimum amplitude of a hyperpolarizing pulse to trigger the F → DB transition decreases with the pulse duration, in agreement with the fact that having longer duration pulses allows for recruiting larger fractions of I CaT owing to the non-instantaneous activation of this current.

Mechanism of SD → F transitions triggered by current injections: crossing the separatrix
From a biological standpoint, the hypothesis of the existence of a silent depolarized state only holds if the DCNn in the SD state are one way or another capable to switch back to their F mode. Otherwise, the entire population of DCNn would eventually become trapped to the SD state. As this prediction contradicts observations of active DCNn in fully developed organisms, our model had (i) to provide evidence for physiological stimuli capable to trigger SD → F transitions and (ii) explain the mechanism of these transitions in order to sustain the hypothesis of the SD state in DCNn. Our standard model does produce such SD → F transitions when hyperpolarizing (Fig. 1D 1 ) or depolarizing pulses (Fig. 1D 2 ) are injected. We next studied the mechanism of this transition. Figure 6A shows that a hyperpolarizing pulse must have a minimum amplitude in order to trigger a SD → F transition. However, this minimum amplitude is clearly not a simple threshold value since Fig. 6B shows that if its amplitude is too large, a hyperpolarizing pulse fails to trigger the switch because it de-inactivates I CaT (Fig. 6B, the first two pulses are large enough to activate 40 to 80% of this current). On the other hand, the mechanism of the SD → F transition triggered by pulses of moderate amplitude like the third pulse in Fig. 6A cannot involve I CaT . Indeed, this current remains inactivated during and after the pulse (<0.15 % activation, Fig. 6A) so one does not expect significant changes in the bifurcation diagram due to I CaT deinactivation. In this case, the transition can only be explained by the fact that hyperpolarizing pulses of moderate amplitude bring the trajectory of the model across the separatrix between the F and SD attraction basins. After entering the F state attraction basin, the model eventually converges to this state. Therefore, to trigger a SD → F transition in the model, the amplitude of a hyperpolarizing pulse must be within an effective range: it must be large enough to drive the system on the other side of the F-SD separatrix but weak enough so as not to de-inactivate I CaT .
A fundamental consequence of the above analysis is the prediction that any perturbation driving only once the trajectory across the separatrix should be able to trigger a SD → F transition. The result illustrated in Fig. 6C confirms this prediction as it shows that not only hyperpolarizing currents but also depolarizing current pulses can trigger a SD → F transition. Altogether, these results support the hypothesis that crossing the separatrix between the basins of attraction of the F and SD state is the basic mechanism by which phasic inputs, either depolarizing or hyperpolarizing, trigger a SD → F state transition. As shown above, altering the calcium currents induces qualitative changes of the bifurcation structure (i.e. changes of the size of the basins of attraction or translations of the bifurcation points), but preserves the overall structure, i.e. the coexistence of a branch of stable silent states with an isola of limit cycles. Therefore, the origin of the peculiar bifurcation diagram of the standard model illustrated in Fig. 2A does not rest on the expression of Ca currents in the DCNn but has to be searched in the interactions between the remaining active currents in the model, namely I NaV , I Kdr and I TCN .

The {V, m Kdr , h NaV } sub-system and the role of I TCN in stabilizing the SD state
We next investigated how the interactions between I Na , I Kdr and I TCN explain the coexistence of a branch of FP with an isola of LC in the DCNn. We have shown above that neither the Ca currents nor the Ca-dependent K currents are required to produce the basic bifurcation diagram of Fig. 2A. After withdrawing these currents, the standard model reduces to a simplified model with only three variables: the membrane potential V , the activation variable of I Kdr , m Kdr and the inactivation variable of I NaV , h NaV . We therefore refer to this reduced model as to the {V , m Kdr , h NaV } sub-system. Figure 7A shows that this sub-system retains the basic dynamical features of the standard model. In particular, the FP branch is nearly identical in the two models. Moreover, the {V , m Kdr , h NaV } sub-system also exhibits an isola of LC albeit in a narrower range of tonic currents. Simulations of the reduced model reveal that pulses of injected current can trigger reversible SD ↔ F like those described above in the full model (not illustrated). The FP branch exhibits a marked sensitivity to the magnitude of the g TCN conductance. Reducing this parameter from its standard value (45 μScm -2 ) shifts to the right the entire branch of FP along the I DC axis. This result stems from the value of the I TCN reversal potential (E TCN = -34 mV) that defines this current as inward (depolarizing) for all points of the FP branch, since they all are below E TCN . Hence, reducing g TCN requires increasing the value of I DC to achieve a FP with the same voltage. As g TCN is reduced to 11.75 μScm -2 , the FP branch exhibits a point of infinite slope (red dot in Fig. 7B). This point corresponds to a cusp bifurcation. Below the cusp bifurcation, the FP branch is no longer stable for every I DC value, but exhibits two saddle-node bifurcations (SN 1 and SN 2 ) that form a hysteresis loop (green rectangle in Fig. 7B) for the example of g TCN = 0). Hence, the model has two coexisting stable silent modes of activity for values for g TCN values below the cusp bifurcation. The isola of LC still exists in the {V , m Kdr , h NaV } sub-system without I TCN (Fig. 7C). But the

Reduction of the model to 2D: evidence for a central role of I NaV in setting the DCNn electric personality
The properties of the {V , m Kdr , h Na } reduced model suggest that the interactions between I NaV and I Kdr play a crucial role in the emergence of the basic bifurcation diagram of the standard model. However, this sub-system does not allow a deep understanding of these interactions owing to their nonlinear nature. This prompted us to search for an even stronger reduction of the standard model to two dimensions, i.e. the minimum dimension for the emergence of limit cycles in ODE systems. To achieve this goal, we followed the approach of Fitzhugh [48] and searched for a relation between m Kdr and h NaV along the trajectory of a spike in the model ( [48,49] and see Rinzel [50] and Gerstner et al. [31] for developments of Fitzhugh's mathematical approach on biophysical grounds). Fitzhugh observed that, in the squid axon, the projection of a spike trajectory in the (m Kdr Oh NaV ) plane is close to a straight line. This allowed one to lump together m Kdr and h Na into a socalled recovery variable, w, exploiting their linear relation interdependence. Combined with the assumption of instantaneous equilibrium for m NaV , this process allowed for the reduction of the original four-dimensional Hodgkin and Huxley model [16] to a simpler, two-dimensional model that keeps the fundamental dynamical features of the original model. Nevertheless, the implementation of this method of reduction of dimensionality also requires an appropriate expression for the (voltage-dependent) time constant τ w for the new variable w. The solution was simple in the squid axon because the voltage dependencies of τ m Kdr and τ h NaV have similar shapes so that the voltage dependency of τ w can be set as a trade-off between that of τ m Kdr and τ h NaV (see Fig. 17, Chap. 2 in [22]). On the opposite, in our DCNn model, the magnitudes and voltage-dependent characteristics of τ m Kdr and τ h NaV are too different to be reconciled by a single, trade-off behavior model for τ w (Fig. 8A). We therefore turned to another approach and built 2d models that preserve either τ m Kdr or τ h NaV . In our model the trajectory of a spike in the {V , m Kdr , h NaV } model projected into the (m Kdr Oh NaV ) plane is not a simple curve (Fig. 8B). Like the projected curve of the Hodgkin-Huxley model, it exhibits two double points (see Fig. 6 in [49]) and we therefore searched for the best nonlinear fit of the trajectory in the (m Kdr Oh NaV ) plane. Firstly, we built a 2d model preserving the properties of I Kdr by fitting the h NaV versus m Kdr relation. We found that the power law h NaV = 0.07/m 0.65 Kdr provides a good fit of the mean trajectory of a spike in the (m Kdr Oh NaV ) plane (not illustrated). By substituting this expression for h NaV into the I NaV equation we obtained a first 2d reduced model, the {V , m Kdr } sub-system, that comprises a unique voltage-dependent time constant τ m Kdr . However, this model fails to capture the features of the {V , m Kdr , h NaV } model given that it fails to produce repetitive firing (not illustrated). We therefore turned to the alternative and fitted the m Kdr versus h NaV relation. We found that this relation is well fitted by the equation m Kdr = 0.0145/h 1.6 NaV . After substituting m Kdr by this equation into the I Kdr current equation, ones obtains a second 2d reduced model, the {V , h NaV } sub-system. Like the {V , m Kdr } sub-system, it contains a unique voltage-dependent time constant but it is τ h Na instead of τ m Kdr . Figure 8C shows that the {V , h Na } sub-system retains the qualitative dynamical features of the {V , m Kdr , h NaV } model with an isola of LC coexisting with a branch of stable FP. The global properties of the {V , h NaV } sub-system can be understood geometrically since the model is 2d. Figure 8D depicts characteristic trajectories of the model for I DC = 0. The blue closed curve is the trajectory of the uLC corresponding to I DC = 0 in the bifurcation diagram. We found that any trajectory starting from an initial condition located inside the region bounded by the blue curve converges onto the FP point (see for example the red trajectory in Fig. 8D). On the opposite any trajectory starting outside the region defined by the blue curve converges onto the corresponding sLC in the bifurcation diagram (see for instance the green trajectory in Fig. 8D). We conclude that the blue curve is the separatrix between the attraction basins of the FP and of the sLC. Moreover, we observed that the {V , h NaV } sub-system loses the characteristic features of the {V , m Kdr , h NaV } reduced model if τ m Kdr rather than τ h NaV is used in the I NaV equation: the isola of LC disappears whereas the FP branch remains unaffected. This shows that the characteristics of the voltage-dependent time of inactivation of I NaV are also fundamental for the coexistence of an isola of LC and of a branch of FP. Together these results suggest that the most important current for these dynamical features is I NaV .

The triggering of reversible F ↔ SD transitions by synaptic inputs
From a functional standpoint, it is crucial to determine whether F ↔ SD transitions in DCNn actually occur in response to synaptic inputs provided by MF and PC. Indeed, the ideal phasic current sources used in previous sections have an infinite resistance whereas synaptic currents have a nonzero conductance that may prevent state transitions to occur by shunting the neuron membrane. Moreover, synaptic currents have an inversion potential which sets limits to voltage changes that can be triggered by activating these conductances. We therefore investigated this question with a variant model in which ideal phasic currents were replaced by models of synaptic currents, i.e. current sources having a nonzero conductance and an inversion potential: I syn = g syn (V -E syn ). Experiments estimate the Nernst potential of Clions in DCNn to be E Cl = -75 mV [3,51]. Figure 9A illustrates the response of the standard model to a pair of successive pulses of synaptic current that model the inhibitory inputs of PC synapses onto the DCNn. Both pulses (of different amplitudes) interrupt the tonic firing of spikes in the model in agreement with experimental observations that PC inputs can silence DCNn [51]. After the pulse, the DCNn, however, resumes firing spikes after a ∼300 ms period of increased firing frequency. We found that such synaptic inhibitory pulses cannot trigger a F → DB transition in those conditions, whatever the duration or amplitude of the pulse (data not shown). Examination of I CaT dynamics during the pulse showed that inhibitory synaptic conductances deinactivate a smaller amount of I CaT compared to direct current injections, thus explaining the absence of F → SD transition by inhibitory synaptic conductances. These results may lead one to conclude that physiological inputs from PC are unable to trigger F → SD transitions in DCNn. However, Boehme et al. [52] have observed that large synaptic inhibition of DCNn by Purkinje cells potently elicit rebound potentials that are underlain by I CaT .
These authors suggested that these observations could result from an erroneous estimate of either the voltage dependence of I CaT or the value of the chloride Nernst potential E Cl .
They noticed that the previous estimates of the voltage characteristics of I CaT have been derived from DCNn recordings made in the soma where the density of I CaT is smaller than in the DCNn dendrites [14]. Owing to difficulties to achieve perfect space-clamp of the membrane voltage over the entire neuron surface, these authors estimated that the actual voltage-dependent parameters of I CaT could be less negative by as much as 10 mV. Figure 9B illustrates the model response to the same stimulation protocol as Fig. 9A but after a +10 mV shift of the voltage-dependence of I CaT . The first current pulse now is able to trigger a F → SD transition while the second pulse induces the opposite transition, switching back the model to its F state. Boehme et al. [52] further suggested that E Cl in the DCNn may also be wrongly estimated by whole cell recordings since the patch pipette imposes its electrolytes composition to the recorded neuron. They estimated that the physiological E Cl could be more negative, down to -85 mV. This negative shift of E Cl induces effects that are qualitatively identical to positive shifts of the I CaT voltage-dependency (Fig. 9C), the first pulse triggering a F → SD transition and the second pulse resetting the neuron to its firing state.

Discussion
The mathematical model presented in this study sketches an original electric personality for the neurons of the deep cerebellar nuclei (DCNn) which remained unsuspected until now in spite of experimental data evidencing puzzling state transitions in these neurons. The pioneer investigation of DCNn intrinsic electric properties by Jahnsen ([2] and [51]) reported that spontaneously active DCNn can be switched to a silent state and back to their firing mode by pulses of injected current (Fig. 5 in [2]). Jahnsen referred to the silent state as a 'steady depolarization and inactivation of the spikes' and our terming 'SD sate' stands as a shortcut for this description. The study of Raman et al. [4] subsequently reported strikingly similar results with patch-clamp recordings showing that depolarizing current pulses can switch DCNn up to a silent depolarized state until the membrane is actively hyperpolarized by a current injection. Notice that the calling of this state a 'depolarization block' (DB) by Raman et al. [4] actually appears misleading. Indeed, many (if not all) neurons exhibit an inactivation of their spike mechanism in response to large depolarizing currents which is called DB (see e.g. midbrain dopamine neurons [53], hippocampal CA1 pyramidal neurons [54] and layer 5 pyramidal neuron of rat's visual cortex [43]). However, the DB does not outlast the duration of its triggering depolarizing current in all cases that we aware of, showing that the DB is not a self-sustained state. In order to highlight these differences, we have chosen to term lasting depolarizations of DCNn as a 'SD state' , this acronym standing for a shortcut of properties of these lasting sustained depolarizations initially described by Jahnsen [2]. Experimental data of Jahnsen [2] and Raman et al. [4] therefore suggest that the SD state may be a genuine state of the DCNn distinct from their spontaneously firing regime (F mode). Since a previous model of DCNn intrinsic electrical properties did not reproduce the coexistence of distinct firing regimes [7], we have built a new model to investigate the hypothesis of a stable depolarizing state in the DCNn, its relation with the spontaneous firing state and interactions between membrane ion currents that give rise to the coexistence of these two states.
Our model uses the latest experimental findings on active membrane ion currents in DCNn to suggest that DCNn have two coexisting stable electrical states. One them is the F mode in which neurons fire fast spikes at low spontaneous frequencies in agreement with in vitro studies [2][3][4][5][6][7]. The other state is a silent mode of activity characterized by a depolarized membrane potential at around -38 mV. We call it the SD state (for stable depolarized state) to stress the proposition that it is locally stable and excitable according to our model. This is evidenced by the fact that our model can be switched back and forth between its SD and F modes by brief pulses of de-or hyperpolarizing currents (see Fig. 1). According to our results (Fig. 2), these bistability properties of the DCNn may reflect the fact that these neurons lay at rest between two-fold bifurcations delimiting an isola of limit cycles which coexists with a branch of stable fixed points. The FP for I DC = 0 would correspond to the plateau depolarization reported by Jahnsen [2] and by Raman et al. [4]. To our knowledge, such dynamical properties have been previously described only in the squid axon bathed in low [Ca] 0 saline [12].
Our model is a single-compartment, isopotential simplification of a cell that actually exhibits a more complex morphology with dendrites and soma. This simplification is in part justified by previous modeling studies of the passive properties of these cells that suggest that these cells are moderately electrotonically compact [7]. A first extension of our model consists in lumping the whole dendritic tree into a unique effective isopotential compartment and connecting it to our initial soma model with a coupling conductance. Preliminary investigation of the resulting two-compartment system suggests that the capacitive and passive resistance loads imposed by the dendritic compartment to the soma are unlikely to challenge the basic mechanism that we propose to explain the transitions from the firing to the SD state in DCNn. In fact, the documented enrichment of DCNn dendrites in CaT current (14) could even boost the facility with which hyperpolarizing inputs can trigger these transitions.

Experimental testing of the model predictions
Our analysis of the mechanisms underlying F ↔ SD transitions with the standard version of the model and variant versions mimicking pharmacological block of the different membrane ion currents suggests three experiments to unambiguously test the hypothesis of F-SD states bistability in the DCNn. These experiments could be achieved with standard intracellular recording of DCNn in cerebellum slices. Experiment I-evidence for a branch of SD states: our model predicts that recorded DCNn should reach a silent state of activity upon injecting them with sufficiently large amounts of steady depolarizing currents. Notice that this behavior is expected for any neuron type whose rising phase of the spike is underlain by a voltage-dependent Na current exhibiting inactivation. Unlike other neurons however, our model predicts that a DCNn should remain in this silent state when the amplitude of injected currents is decreased back at a rate sufficiently slow to prevent crossing of the separatrix between the F and SD states. Then, starting from hyperpolarized membrane potentials, our model predicts that the recorded DCNn should remain silent when the magnitude of the hyperpolarizing current is slowly reduced down to zero and even upon injecting large depolarizing currents. Experiment II-evidence for an isola of limit cycles: this second experiment investigates the firing mode of the DCNn with the help of standard (firing frequency vs I DC ) curves. Injecting depolarizing steady currents of growing magnitude into spontaneously active DCNn should increase their firing frequency up to a maximum current value (I F2 in Fig. 1) beyond which the DCNn should switch to a silent state. Conversely, in a spontaneously firing DCNn, increasing the magnitude of the tonic hyperpolarizing current above I F1 should abruptly switch the neuron to its silent mode. The observation of this result would lead one to categorize the DCNn as type II neurons according to the Hodgkin classification [42]. The latter part of this experiment should especially be achieved in the presence of broad-spectrum inhibitors of postsynaptic channels to remove spontaneous synaptic activity in slices. Membrane potential fluctuations due to ongoing synaptic activity would hinder determining if the DCNn minimal firing frequency is actually nonzero. Experiment III-changing the external Ca 2+ concentration: this experiment tests the model prediction that the DCNn can be made to adopt the classical bifurcation scenario for type II neurons, by turning the left-hand side fold bifurcation of the isola into a subcritical Hopf bifurcation and the right one into a supercritical Hopf bifurcation (see [43]) by increasing [Ca] 0 . These changes of bifurcation would not be accompanied by qualitative modification of the f/I relationship. But they could be identified by a continuous decrease of the spike amplitude down to zero upon increasing the driving current up to the right-hand side bifurcation. Experiment IV-extent of I Na inactivation in the SD state: intracellular recordings could (i) check that I Na is not fully inactivated in the SD state and (ii) investigate whether the full blockage of I Na with TTX prevents DCNn state transitions by removing the SD state.

Ionic currents and dynamics of the model
The analysis of the standard model and its lower dimension variants suggest that the voltage-dependent properties of Na currents are of paramount importance in the state transitions displayed by the DCNn in vitro. In order to justify this major conclusion, we come back to the biophysical foundations of our model, namely the set of ion currents extracted from the experimental literature on the DCNn to build our model and their mathematical formulation. To avoid introducing unjustified hypotheses, we assumed that voltage-dependent currents obey the classical Hodgkin-Huxley formalism which has proven capable to reproduce the qualitative dynamics of most neuron types studied so far with mere parameter adjustments. For I Kdr this hypothesis is consistent with the finding by Raman et al. [4] of a classical delayed rectifier K current in DCNn. However, the large time constant of activation reported in this study (12 ms at +12 mV) suggests that this voltagedependent current is carried by (slow) Kv2 channels (see e.g. [29]) whereas other studies give no evidence for the expression of Kv2 channels by DCNn. These studies rather show that DCNn express (fast) Kv3 channels and that these channels are functional (see Methods). Our model was accordingly designed with a voltage-dependent time constant for I Kdr corresponding to fast Kv channels. Nevertheless, we investigated the possible involvement of Kv2 channels in the DCNn electro-responsiveness with a variant model comprising a slowly activating K current, I KdrS , in addition to I Kdr . I KdrS was given the same voltagedependence as I Kdr but its time constant was multiplied by 11 to achieve a value of 12 ms at +12 mV. Adding I KdrS to the model with the same conductance as I Kdr did not change the overall structure of the model's bifurcation diagram (compare panels A and B in Fig. 10). However, it increased the width of the isola of LCs and reduced the spontaneous firing frequency from 28.9 Hz to 25.5 Hz. This frequency could even be reduced to 20 Hz (the spontaneous firing frequency of DCNn in vitro [2][3][4][5][6][7]) by increasing the g KdrS value to 3.3 times that of g Kdr (not illustrated). These results suggested that previous studies of Kv channels in DCNn may have missed the Kv2 channels implied by the results of Raman et al. [4] and that currents through Kv2 channels underlie the low spontaneous frequency firing of DCNn. To address this hypothesis, we fixed the value of g KdrS and decreased that of g Kdr to determine to what extent slow K V 3 channels are mandatory to explain the observed firing dynamics of DCNn. Decreasing g Kdr diminishes the width of the isola of LC down to g Kdr = 3 × 10 3 μScm -2 , where the branch of SD states loses its stability (Fig. 10C check). Nevertheless, overall stability of the FP branch was restored by increasing the conductance of I TCN to 3 × 10 3 μScm -2 , suggesting that the functional role of I TCN is to stabilize the SD state of DCNn. We tested this proposal by reducing g Kdr down to the complete withdrawal of the fast K V current in the model. The range of tonic currents over which the FP branch is unstable in the hybrid Kdrs-Kdr variant model (Fig. 10B) was even more enlarged after switching the I Kdr model to solve Kv2-like dynamics (Fig. 10E). Nevertheless, increasing g TCN proved again capable to resolve this issue as shown by panel F in Fig. 10 displaying the bifurcation diagram of the model endowed with g TCN to 215 μScm -2 . Beyong strengthening our proposal on the functional role of ITCN, these results demonstrate that the basic bifurcation scenario of the standard model can be achieved by a model comprising only slow Kv channels.
Adequacy of the classical HH formalism to model voltage-dependent Na currents in DCNn proved more challenging owing to two reasons: (i) the integer power to which raise variable m NaV and (ii) the evidence for a resurgent Na current in DCNn (Ref ) not   [55]. B 1 . I NaV replaced by I NaRB with conductance g NaRB = 5 × 10 3 (same as that of I NaV in the standard model). B 2 . g NaRB = 3.5 × 10 3 . B 3 . g NaRB = 3 × 10 3 . B 4 . Same as B 3 with g TCN = 75 (standard value = 45). B 5 . With g Kdr = 1.25 × 10 4 (standard value = 4.5 × 103). B 6 . g NaRB = 3 × 10 3 , g Kdr = 1.25 × 10 4 and g TCN = 9.5 × 10 1 introduced in our model. Raman et al. [4] fitted their data for the steady-state activation of Na currents with a Boltzmann function of the form given by Eq. (5) for m NaV∞ (see their Fig. 3C). The standard HH formalism for this current uses the same function but raised to a power of three (see e.g. [26]), i.e. m 3 NaV∞ . Assuming that the error between the data points and the model prediction are independently and identically distributed according to a normal law, we computed the Bayesian information criterion (BIC) for both models [56]. We found that BIC m ∞ -BIC m 3 ∞ > 16, which means that the standard m 3 NaV∞ formalism yields a better description of the experimental data than a m Na∞ formalism. The classical HH formalism hence appears to be capable to describe adequately the transient component of I NaV currents in DCNn. Notice that the two studies on I NaV in DCNn that we are aware of ( [4] and [18]) provide divergent results regarding the inactivation parameters of this current. Our model adopts as standard parameters the values reported in [18] in physiological salines rather than values in [4] as the latter were obtained in low external Na conditions. Nevertheless, we found that the basic bifurcation diagram and F ↔ SD state transitions are robust to changes in the inactivation parameters of I NaV which include the values reported in [4], thereby showing that these properties do not result from an unrealistically precise setting of the I NaV parameters. We then investigated the possible role of the resurgent component of I NaV reported in DCNn by Afshari et al. [18]. This current, which was initially identified in cerebellar Purkinje cells [55], corresponds to a transient block of NaV channels in their open state that is not accounted for by the classical HH formalism. For this reason, we investigated the properties of a variant model in which I NaV was modeled with the state transition scheme of Raman and Bean [19], which produces a resurgent current. Figure 10B summarizes the results of this investigation. When endowed with I NaRB instead of the standard I NaV (same conductance as I NaV ), the model can no longer produce F ↔ SD transitions due to the loss of overall stability of the FP branch (Fig. 10B 1 ). The exchange of Na current models also induces the appearance of saddle-node (SN) bifurcation in this branch. Rather than forming an isola of LC as in the standard model, the branch of sLC ends into a homoclinic bifurcation at regular saddle at the left of the bifurcation diagram and into a Hopf bifurcation (H 1 ) at the right of the diagram. As I NaRB adds a resurgent component to the transient Na current modeled with I NaV , we reasoned that loss of stability of the FP may stem, at least partly, from the Na current having become excessively large. In agreement with this reasoning, panels B 2 and B 3 in Fig. 10 show that a progressive decrease of g NaRB shifts the H 1 point to the left of the bifurcation diagram. It also removes the SN bifurcation point to leave a second Hopf bifurcation point (H 2 ) which shifts to the right of the bifurcation diagram as g NaRB is decreased furthermore. However, decreasing g NaRB alone proved unable to restore the standard bifurcation diagram as the branch of LCs was lost when the FP resumed overall stability (not illustrated). As our study suggests that the I TCN current contributes to the stability of the FP branch, we then increased g TCN in addition to decreasing g NaRB . Figure 10B 4 shows this additional parameter change brought the H 1 and H 2 closer to each other and allowed the model to recover two short branches of uLC at ends of the branch of sLC. However, further increases in g TCN proved unable for the variant model to recover the isola of LC as they exchanged the stability of the left uLC branch (not illustrated). Given that I Kdr exerts a counteracting effect to I Na , we finally examined the effects of increasing g Kdr . Figure 10B 5 shows that increasing g Kdr has similar effects to that of increasing g TCN . However, increasing g Kdr alone also proved unable for the model to recover the isola of LC (not illustrated). Nevertheless, panel B 6 shows that com-strongly suggest to study computational models of the cerebellum network endowed with the basic DCNn state transition features disclosed by our study.

Appendix: On the separatrix between the basins of attraction of the SD and F states
We show in the Methods section that solutions of our model are bounded. Moreover, our model has no strange (fractal dimension) attractor in addition to the F and SD attractors. It follows that any trajectory of the model, whatever its initial conditions, must eventually converge on either the F or the SD attractor. This implies the existence of a closed manifold (i.e. compact and without bounds) separating the phase space. This leads us to propose that the separatrix delimiting the basins of attraction of the F and SD attractors is a (n -1)dimension closed manifold. We give here information supporting this hypothesis.
A virtue of the reduced {V , m Kdr , h Na } model is that it possesses a three-dimensional phase space allowing one to visualize its dynamics, contrary to the standard six-dimensional system. We used this property to study in more detail the separatrix between the basins of attraction of the SD and F states. We were able to sample the separatrix in the {V , m Kdr , h NaV } model by imposing a fixed value to two its three state variables model and then finding by dichotomy values of the remaining variable that are closest to the separatrix. We illustrate the results in Fig. 11A for three different orientations in the {V , m Kdr , h NaV } phase space. The black mesh locates the points of the separatrix sampled for different fixed values of either m Kdr or h NaV . This mesh clearly forms a closed surface showing that the separatrix in the {V , m Kdr , h NaV } reduced model is two-dimensional. The small parallelogram in panels A represents a flat surface locally tangent to the separatrix at a point in the phase space with coordinates (V = -30.1, m Kdr = 0.3, h NaV = 0.06). The inner side of this surface is depicted in light gray and its outer side in dark gray. Figure 11B illustrates an example of the fate of all trajectories starting from initial conditions located inside the volume delimited by the separatrix: all these trajectories eventually converge to the SD state (red dot). On the opposite, all trajectories of the model starting from a point located outside of this volume converge onto the sLC (green curve) as illustrated by Fig. 11C. These results confirm that the surface delineated by the black mesh in Fig. 11A [1][2][3] is the separatrix between the attraction basins of the F and SD states. Numerical investigation of the uLC in the {V , m Kdr , h NaV } reduced model revealed that this limit cycle has one of its Floquet exponents equal to 1 (as expected from Floquet theory, (see e.g. [13])), a second one being <1 and the remaining one being >1. Hence, the local unstable manifold of the uLC is of dimension 1 whereas its local stable manifold is of dimension 2. Since the global stable and unstable manifolds of the uLC are locally tangent, respectively, to the local stable and unstable manifolds in a neighborhood of the uLC and according to the above results showing that the separatrix is a closed surface, we are led to propose that the SD/F separatrix is actually the global stable manifold of the uLC. The finding that this manifold is of dimension 2 agrees with the hypothesis that the separatrix is of dimension n -1 in the {V , m Kdr , h NaV } reduced model. Moreover, numerical analysis of the uLC in the standard six-dimensional model revealed that one of its Floquet exponents is equal to 1, four of them are <1 and the last one >1 so that the sLC stable manifold is of dimension 5. According to the conclusion derived from the {V , m Kdr , h NaV } reduced model, it follows that the F/SD separatrix is likely of dimension 5 in the standard model, thereby providing support for the n -1 dimension hypothesis on the separatrix. Gray box: flat surface patch tangent to the separatrix; light and dark grays, respectively, indicates inside and outside of the separatrix surface. B. The model converges to its SD state from a state point of the phase space located inside the phase space volume bounded by the separatrix. C. Convergence of the model to the sLC from initial conditions located outside of this volume. All data depicted in the figure were computed with I DC = 0