The geometry of rest-spike bistability

TL;DR

This paper introduces an augmented Morris-Lecar model that captures bistability between resting and spiking states, providing a physiologically interpretable and mathematically tractable framework for studying complex neuronal dynamics.

Contribution

The paper presents a modified Morris-Lecar model with an added slow inward current to accurately model bistability, bridging physiological relevance and mathematical simplicity.

Findings

The augmented model exhibits bistability between resting and spiking states.

It effectively captures slow spiking and bursting phenomena.

The model offers a balance of physiological interpretability and analytical tractability.

Abstract

Morris-Lecar model is arguably the simplest dynamical model that retains both the slow-fast geometry of excitable phase portraits and the physiological interpretation of a conductance-based model. We augment this model with one slow inward current to capture the additional property of bistability between a resting state and a spiking limit cycle for a range of input current. The resulting dynamical system is a core structure for many dynamical phenomena such as slow spiking and bursting. We show how the proposed model combines physiological interpretation and mathematical tractability and we discuss the benefits of the proposed approach with respect to alternative models in the literature.

Equations

$$\epsilon \dot{v}$$

$$\dot{n}$$

$$\tau \dot{p}$$

$$i_{ion}$$

$$=c(v)+n(v+1)+p(v-1)$$

$$S_x=\frac{g_x}{2}\left(\tanh \left(\frac{v-a_x}{b_x}\right)+1\right)$$

$$0$$

$$\dot{n}$$

$$\tau \dot{p}$$

$$v^{\prime }$$

$$n^{\prime }$$

$$p^{\prime }$$

$$i_{ion}(v,n,p)=i$$

$$\frac{\partial i_{ion}}{\partial v}=\frac{dc}{dv}(v)+n+p=0$$

$$\frac{\partial i_{ion}}{\partial v}\dot{v}$$

$$\tau \dot{p}$$

$$\dot{v}$$

$$\tau \dot{p}$$

$$\frac{\partial i_{ion}}{\partial v}=0\frac{\partial i_{ion}}{\partial n}\dot{n}+\frac{\partial i_{ion}}{\partial p}\dot{p}=0$$

$$\frac{\partial i_{ion}}{\partial n}\dot{n}+\frac{\partial i_{ion}}{\partial p}\dot{p}\ne 0$$

$$\frac{{\partial }^2{i}_{ion}}{\partial v^2}\ne 0$$

$$i_s(v):=i_{ion}(v,S_n(v),S_p(v))$$

$${\Pi }_f\circ {\Pi }_l(I_l)\subset I_h{\Pi }_f\circ {\Pi }_h(I_h)\subset I_l$$

$$\Pi ={\Pi }_f\circ {\Pi }_h\circ {\Pi }_f\circ {\Pi }_l:I_l\to I_l$$

$$i_s(v)=i_{ion}(v,S_n(v),S_p(v))$$

$$\epsilon \dot{v}$$

$$\dot{s}$$

$$i=i_{ion}(v,n(s),p(s))$$

$$n=n(s)p=p(s)$$

$$i_{ion}$$

$$\frac{\partial i_{ion}}{\partial v}$$

$$\frac{\partial }{\partial s}(i_{ion}(v,n(s),p(s)))$$

$$i=i_{ion}(v,n,p)\frac{\partial i_{ion}}{\partial v}(v,n,p)=0$$

$$\Pi \left(\begin{array}{c}v\\ p\end{array}\right)=\left(\begin{array}{c}R(v,p,i,\epsilon )\\ G(v,p,i,\epsilon )\end{array}\right)$$

$$R(v,p,i,0)=0$$

$$G=G_0(p)+\mathcal{O}(\epsilon \ln (\epsilon ))$$

Full text

1 Introduction

Conductance-based models are by now well established as a fundamental modeling framework to connect the physiology and the dynamics of excitable cells. Ever since the seminal work of Hodgkin and Huxley [14], there has been a continuing effort in the literature to develop models that combine mathematical tractability and physiological interpretation. An interesting example is the two-dimensional model published by Morris and Lecar in 1981 [19]. Like Hodgkin-Huxley model, it captures the essential physiology of excitability: a spike results from the fast activation of an inward current followed by the slow activation of an outward current. The former provides positive feedback in the fast time-scale whereas the latter provides negative feedback in the slow time-scale. Because it is only two-dimensional, the model is also amenable to phase-portrait analysis without any reduction. Its geometry is similar to the one of FitzHugh-Nagumo model [6], the first mathematical model proposed to understand the core dynamics of the Hodgkin-Huxley model. In that sense, Morris-Lecar model combines the physiological interpretation of Hodgkin-Huxley model and the mathematical tractability of Fitzugh-Nagumo circuit.

In the present paper, we aim at capturing in a similar way the essence of rest-spike bistability, that is, the coexistence of a stable spiking attractor and a stable fixed point in a slow-fast model. The importance of this phenomenon is well acknowledged in the neurodynamics literature due to its role as a building block of neuronal patterns such as bursting [4, 16]. We obtain rest-spike bistability by adding one extra current in Morris-Lecar model: an inward current with slow activation. The resulting model combines the three following features:

1. i)

For a range of input currents, the model is rest-spike bistable, that is, a stable equilibrium coexists with a stable limit cycle. The geometry of the two attractors is robust to the time-scale separation in the sense that it persist in the limit of infinite time-scale separation. 2. ii)

The model has the direct physiological interpretation of dynamics and attractors being shaped by three distinct currents (fast positive feedback (e.g. sodium activation), slow negative feedback (e.g. potassium activation), and slow positive feedback (e.g. calcium activation)). 3. iii)

The model is amenable to a mathematical analysis by geometric singular perturbation theory.

We are not aware of other single-cell models in the literature combining those three features. Mathematical models of rest-spike bistability often lack the first feature above. For instance, a homoclinic bifurcation in the Morris-Lecar model only exists for a specific time-scale separation (see e.g. Table 3.1 in section 3.2 of [4]). Limitations of such models with respect to the geometry of the attractors and the robustness of bursting are discussed in [7, 9]. We are only aware of two published models in which rest-spike bistability persists in the singular limit of infinite time-scale separation. The first one is the model proposed by Hindmarsh and Rose in 1982 [12] as a mathematical model aimed at capturing low-frequency spiking. The second one is the transcritical model proposed in 2011 as a two-dimensional reduction of a physiological model combining the currents of the Hodgkin-Huxley model with a slow inward (calcium) current [2]. Both models are planar and lack the second feature, that is, they can only be regarded as a mathematical reduction of a physiological conductance-based model.

This paper aims to contribute to the idea that balancing positive and negative feedback in the slow time scale is a key mechanism to generate rest-spike bistability. This viewpoint is at the core of the planar model in [7] and its importance from a physiological viewpoint is highlighted by [9]. Here we complement those works by studying how this mechanism can be naturally implemented in a physiological context: using two distinct slow currents, one providing negative feedback to restore the membrane potential, the other providing positive feedback to obtain two attractors separated by the stable manifold of a saddle.

The remainder of the paper is organized as follows. Section 2 presents the model and recalls the notions of geometric singular perturbation theory needed for its analysis. In section 3 we study numerically the dynamics on the critical manifold, highlighting its persistence properties. Section 4 builds on this picture to derive conditions for multistability and monostability, we focus on the singular case and mention what hypotheses guarantee persistence. In section 5 we discuss some variations of the same geometric picture, while in section 6 we relate the Hindmarsh-Rose and the transcritical models to the one we are studying. We draw some conclusions in section 7. Two appendices report additional details.

2 A Model of Rest-Spike Bistability

We consider a three-dimensional slow-fast conductance-based model defined by

[TABLE]

where $\varepsilon$ is a small parameter. The total ionic current is the sum of a leak current and three voltage-gated currents:

[TABLE]

The parameters $-1<v_{l}<+1$ that appear in the equation can be thought of as reversal potentials. In the absence of an external current $i$, the voltage range $[-1,1]$ is positively invariant. The currents $S_{m}(v)(v-1)$ and $p(v-1)$ are then negative (inward currents) whereas the current $n(v+1)$ is positive (outward current). The variables $n$ and $p$ are gating (positive) variables that model the slow activation of the inward current $p(v-1)$ and of the outward current $n(v+1)$. The inward current $S_{m}(v)(v-1)$ has instantaneous activation, a standard simplification for currents that activate in the fast time-scale. The functions $S_{x}(v)$ correspond to activation functions that we assume of the form

[TABLE]

Here the multiplicative factor $g_{x}$ corresponds to the maximal conductance associated to the current $x$. We find it convenient to write the equations in this form, rather than including maximal conductances in the voltage equation, because this allows us to change maximal conductances of slow currents without modifying the critical manifold of the system. A consequence of this is that the dynamics of $p$ and $n$ lies between zero and the corresponding maximal conductance, i.e. $p\in[0,g_{p}]$ and $n\in[0,g_{n}]$.

The key property of the model is the presence of the slow inward current $p(v-1)$. In the absence of this current, the model is two-dimensional and has a phase portrait similar to the classical FitzHugh-Nagumo model. With this additional slow inward current, both continuous spiking and rest coexist for the same value of applied current, as shown by the simulation in fig. 1 (see appendix B for numerical values of the parameters).

We note that similar phenomena can be obtained with a current of the type $p(v+1)$ where $p$ inactivates, i.e. decreases as $v$ increases. Physiologically this corresponds to an outward current that inactivates slowly, rather than an inward current that activates slowly. Both types of currents model a source of positive feedback in the slow time-scale [10]. A classical example of slowly inactivating outward current is the $A$-type potassium current [3].

We use geometric singular perturbation theory [5] to study the slow-fast system eq. 1 as $\varepsilon$ tends to zero. The singular limit of this model is the differential-algebraic system

[TABLE]

which we call slow dynamics or reduced system. After rescaling time, the same limit leads to the layer dynamics

[TABLE]

where ′ refers to differentiation with respect to the fast time $\tau=t/\varepsilon$.

The reduced system eq. 4 is constrained to the critical manifold $\mathcal{C}_{0}$, defined by

[TABLE]

and corresponding to fixed points of the layer dynamics eq. 5. Normally-hyperbolic compact subsets of $\mathcal{C}_{0}$ persist as invariant manifolds of eq. 1 for $\varepsilon$ small enough. This manifolds are not necessarily unique, but we assume one family of perturbation $\mathcal{C}_{\varepsilon}$ has been fixed and call them slow manifolds.

Perturbations of subsets of $\mathcal{C}_{0}$ maintain their type of stability with corresponding (local) stable and unstable manifolds. These admit invariant foliations, with each point on the critical manifold acting as base for a fiber. Invariance of the foliation can be interpreted as points on each fiber “shadowing” the corresponding base point, in forward time for the stable manifold and backward for the unstable. Points on $C_{\varepsilon}$ follow a dynamics that is a regular perturbation of the reduced system eq. 4; in the following we refer to this perturbation as slow dynamics.

A point $x$ on the critical manifold is normally hyperbolic if it is a hyperbolic fixed point of the layer dynamics eq. 5. If this is the case, as $\varepsilon\to 0$ the fibers based at $x$ tend to its stable and unstable manifolds in the layer dynamics eq. 5. For eq. 1 the layer dynamics is one dimensional, so that hyperbolic fixed points are either attractive or repulsive, with their invariant manifolds corresponding to lines with $n$ and $p$ constant.

We consider parameter ranges for which the critical manifold can be divided in three normally-hyperbolic branches. These are separated by two lines of folds that we call $F_{l}$ and $F_{h}$, and verify

[TABLE]

The two lines of folds are connected by an unstable branch $M$. The other branches, $S_{l}$ and $S_{h}$, are both stable. Figure 2 shows the typical shape of the critical manifold for fixed $i$.

Around any point away from the lines of folds, the critical manifold admits a parametrization in the slow variables $n$ and $p$. However, this local parametrization cannot be made global due to the presence of folds. Following [23, 24, 25], we use $v$ and $p$ to obtain a parametrization valid in the interval $v\in(-1,\,1)$. This is achieved by solving eq. 6 for $n(v,p,i)$. The corresponding projection is shown in fig. 2.

The reduced dynamics in these coordinates is obtained differentiating eq. 6:

[TABLE]

The first equation becomes singular on the lines of folds $\frac{\partial i_{ion}}{\partial v}=0$. Multiplication by $\frac{\partial i_{ion}}{\partial v}$ recovers a regular differential equation

[TABLE]

The two systems eq. 8 and eq. 9 share the same trajectories with different time parametrizations. Moreover, in eq. 9 time is reversed on the unstable branch $\frac{\partial i_{ion}}{\partial v}<0$ and new fixed points can appear on the lines of folds. These verify

[TABLE]

They are called folded singularities [23].

Away from the lines of folds the two systems eq. 8 and eq. 9 are largely equivalent, but important differences occur in the neighborhood of $F_{l}$ and $F_{h}$. Moreover, near these lines the perturbed dynamics is no longer constrained by normal hyperbolicity, in particular it cannot be obtained as a regular perturbation of the reduced system eq. 8. Different phenomena are possible. The least degenerate situation occurs when the desingularized vector field is never zero along these lines:

[TABLE]

Under this assumption the desingularized vector field eq. 9 can point either to the unstable branch or the stable one. Assuming the additional nondegeneracy condition

[TABLE]

the first case corresponds to jump points, at which the reduced system eq. 8 admits two solutions backwards in time but none in forward time. For $\varepsilon>0$ a stable branch of $\mathcal{C}_{\varepsilon}$ near these points can be continued using the flow [24]. Doing so shows that trajectories on the slow manifold pass the folds and reach a fiber contained in the stable manifold of the other stable branch of $C_{\varepsilon}$, with the flow contracting the direction transverse to the manifold.

Condition eq. 11 corresponds to the vector field being transverse to the critical manifold, a condition which is violated at folded singularities. These are fixed points of the desingularized system eq. 9, but not necessarily fixed points of the reduced dynamics eq. 8. As a consequence, they can be reached in finite time. Depending on the type of fixed point they can correspond to the singular limit of canard trajectories, i.e. intersections between stable and unstable branches of the slow manifold [23]. Generically, the desingularized flow changes direction at these points. Hence, a folded singularity delimits the set of jump points on a line of folds [25].

3 Reduced Dynamics

We will now study the reduced system eq. 8, often with the aid of its desingularized version eq. 9. Fixed points can be parametrized by $v$ through the steady-state $i$-$v$ curve

[TABLE]

This is shown in fig. 6 and is an S-shaped curve, with two folds separating three families of fixed points $\mathcal{X}_{l}$, $\mathcal{X}_{m}$ and $\mathcal{X}_{h}$; $\mathcal{X}_{l}$ corresponds to low voltages, $\mathcal{X}_{m}$ to intermediate voltages and $\mathcal{X}_{h}$ to high voltages. For fixed $i$, we denote points in each family with corresponding lower-case letters $x_{l}$, $x_{m}$ and $x_{h}$. In addition to these three fixed points, the desingularized dynamics eq. 9 has a folded singularity $x_{f}\in F_{l}$. For parameter values reported in appendix B, and $i$ in the range of interest in this section, this point is a focus and does not lead to canard trajectories [23]; it only delimits jump points on $F_{l}$.

Figure 3 shows the typical phase portrait of the reduced system eq. 8. The fixed point $x_{l}$ is a stable node, while $x_{m}$ and $x_{h}$ are both saddle points. Their stable and unstable manifolds do not extend beyond $F_{l}$ and $F_{h}$ due to loss of existence and uniqueness along these lines. In particular, unstable manifolds terminate at jump points.

For $\varepsilon>0$ hyperbolic fixed points persist in the slow dynamics with their stable and unstable manifolds [5, 22]. In the perturbed system eq. 1 these fixed points are still hyperbolic. In particular, saddle points remain saddle, with their invariant manifolds being obtained as a combination of trajectories in the slow dynamics and fast fibers. The unstable manifold of $x_{m}$ is completely contained in the slow manifold. Its stable manifold, instead, is two dimensional; it includes the stable manifold in $\mathcal{C}_{\varepsilon}$ as well as all fast fibers based on that curve. In the singular limit this surface tends to the stable manifold of $x_{m}$ in the reduced system eq. 8 and all nearby segment with constant $n$ and $p$ that intersect it. Similarly, $x_{h}$ perturbs to a saddle with a one-dimensional stable manifold and a two-dimensional unstable manifold.

Adding a trivial equation for $i$ to eq. 1, the same is true for the family of fixed points $x_{m}(i)$. At least for $i$ in a small interval, this family persists together with its two-dimensional unstable manifold and three-dimensional stable one. Sections of these manifolds for fixed $i$ coincide with the invariant manifolds of the corresponding fixed point.

When $i$ varies on larger domains, the outlined phase portrait can undergo two distinct qualitative changes: increasing $i$ leads to a fold of the $i$-$v$ curve at which $x_{l}$ and $x_{m}$ merge in a saddle-node bifurcation, leaving only one fixed point $x_{h}\in M$. Likewise, decreasing $i$, $x_{m}$ and $x_{h}$ reach a similar fate, leaving $x_{l}\in S_{l}$ as the only fixed point.

To obtain this second bifurcation it is necessary that one of the two fixed points crosses the line $F_{l}$ and changes branch111This assumes that the bifurcation does not happen exactly on $F_{l}$.. In our case $x_{m}$ crosses $F_{l}$. This passage corresponds to an exchange of stability with the folded singularity through a folded saddle-node [17]. Beyond this crossing, the folded singularity is a saddle, while $x_{m}$ is a node of the reduced system. In a similar fashion increasing the applied current leads to $x_{h}$ crossing $F_{h}$, which happens once $x_{h}$ is the only fixed point left. After this crossing, $x_{h}$ is a stable fixed point on an attractive branch.

Finally, varying $i$ can lead to changes in the type of folded singularity. As already mentioned $x_{m}\in F_{l}$ corresponds to a folded saddle-node, thus varying $i$ and moving $x_{m}$ between branches leads to different types of folded singularity: it is a saddle when $x_{m}\in M$ and a node when $x_{m}\in S_{l}$. Both situations lead to canard trajectories [23]. Moreover, since $x_{f}$ is a focus in the phase portrait described above, it has to change to a node before becoming a folded saddle-node.

4 Rest-spike bistability

Returning to the phase portrait in fig. 3, we now analyze the global return mechanism that leads to rest-spike bistability.

In the singular limit $\varepsilon=0$, trajectories on stable branches of the critical manifold $\mathcal{C}_{0}$ stay on it until they reach a line of fold in correspondence of a jump point. Once one of these points is reached the singular trajectory is continued along a fast fibre with constant $n$ and $p$, reaching the opposite branch as shown in fig. 2. The points at which these singular trajectories arrive correspond to the projections of $F_{l}$ and $F_{h}$ along fast fibers. We call these projections $P_{l}\subset S_{l}$ and $P_{h}\subset S_{h}$.

Based on this property, we can analyze the singular system referring only to the $v$-$p$ plane and the reduced dynamics: when a trajectory reaches a jump point it is transported to the corresponding projection keeping $p$ fixed, as shown in fig. 2 for a limit cycle.

Rest-spike bistability follows from how the stable and unstable manifold of $x_{m}$ constrain trajectories. The role of the stable manifold is simple, it separates initial conditions on $S_{l}$ that reach a jump point on $F_{l}$ from those that remain on the critical manifold and tend to $x_{l}$. The unstable manifold, instead, determines if the system is multistable. This is the case if the unstable manifold stays away from $x_{l}$. Otherwise almost all trajectories converge to $x_{l}$. We treat these two situations separately in the next sections.

4.1 Bistability

In the following we denote by $x_{1}$ the intersection of the unstable manifold of $x_{m}$ with $F_{l}$, and by $x_{-1}$ the intersection of the stable manifold of $x_{m}$ with $P_{l}$. Following the singular flow from $x_{1}$ leads to $x_{2}\in P_{h}$, then to $x_{3}\in F_{h}$ and back to $P_{l}$ at $x_{4}$ (see fig. 4). We recall that given the dynamics eq. 1 we can assume that $p$ lies in the interval $[0,g_{p}]$, where $g_{p}$ is the maximal conductance appearing in eq. 3.

Assume that the trajectory starting at $x_{4}$ reaches a jump point on $F_{l}$ ($x_{5}$), as shown in fig. 4. Consider the segment $I_{l}\subset P_{l}$ between $x_{4}$ and $p=g_{p}$. The reduced dynamics maps this segment to $F_{l}$ in finite time, defining a map $\Pi_{l}:I_{l}\to F_{l}$. Clearly the same map can be defined using the desingularized reduced system eq. 9, thus as long as this vector field is transverse to $F_{l}$ at all points in $\Pi_{l}(I_{l})$ the map is smooth. We note that this is equivalent to $\Pi_{l}(I_{l})$ not containing folded singularities. Similarly, on $S_{h}$ we define the segment $I_{h}\subset P_{h}$ between $x_{2}$ and $p=g_{p}$, and a corresponding map $\Pi_{h}:I_{h}\to F_{h}$. We denote the projection along fast fibers by $\Pi_{f}$ (from $F_{l}$ to $P_{h}$ and from $F_{h}$ to $P_{l}$). Since the dynamics is bounded by the line $p=g_{p}$, by construction we have

[TABLE]

which allows us to define the singular Poincaré map

[TABLE]

This construction shows that the stable manifold of $x_{m}$ divides the state space into two invariant sets. One is the basin of attraction of $x_{l}$, while the other one has dynamics characterized by the Poincaré map eq. 15. Since this is a smooth map of an interval into itself it admits at least one fixed point, which corresponds to a singular relaxation oscillation. As shown in [24], if this fixed point is hyperbolic, under the additional hypothesis that the singular trajectory intersects $P_{l}$ and $P_{h}$ transversally, it perturbs to a hyperbolic limit cycle for $\varepsilon>0$. In fact the Poincaré map eq. 15 is (up to conjugacy) a global version of the one used in that reference.

We remark that this construction only guarantees multistability. Further analysis of the map eq. 15 is required to obtain a more accurate picture. While this is beyond the scope of this work, numerical simulations confirm that this map has a unique attracting fixed point.

4.2 Monostability

Constructing the Poincaré map eq. 15 requires that $x_{4}$ falls inside the interval defined by $x_{-1}$ and $p=g_{p}$ on $P_{l}$. The situation in which this assumption fails is illustrated in fig. 5. In this case most trajectories on $S_{l}$ and $S_{h}$ are attracted by the stable fixed point $x_{l}$, the only exception being the stable manifold of $x_{m}$.

To see this, we start from $x_{-1}$ and consider its anti-image through $\Pi_{f}$ on $F_{h}$. Continuing to follow the singular flow “backwards”, as shown in fig. 5, leads back to $P_{l}$ at a point that we call $x_{-2}$. Any compact segment in $P_{l}$ that lies between these points is mapped by the singular flow strictly inside the segment delimited by $x_{4}$ and $x_{-1}$. Since any point strictly inside this second segment converges to $x_{l}$, the same conclusion extends to all points in the original segment.

The same argument shows that points in the portion of $S_{l}$ delimited by the trajectories starting at $x_{-1}$ and $x_{-2}$ tend to $x_{l}$. The only exceptions are these boundary trajectories that reach $x_{m}$ and belong to its stable manifold. As long as the stable manifold of $x_{m}$ is unbounded in the $p$ coordinate, the same argument can be iterated on all $S_{l}$ and adapted to $S_{h}$, leading to the conclusion that almost all points on $S_{l}$ and $S_{h}$ are in the basin of attraction of $x_{l}$. This situation persists for small enough $\varepsilon>0$, and since most points are attracted to stable branches of the slow manifold we obtain that for almost all initial conditions the perturbed dynamics converges to $x_{l}$.

4.3 Homoclinic trajectory and bifurcation diagram

Transitions between monostability and bistability in system eq. 1 are controlled by the applied current $i$. The phase portraits in figs. 4 and 5 suggest the presence of a homoclinic trajectory, which can be obtained by decreasing the applied current from the bistable case. In the singular limit this trajectory corresponds to the condition $x_{4}=x_{-1}$ and delimits the boundary of bistability. We denote by $i_{H}$ the value of current at which this happens. While we cannot expect this homoclinic trajectory to persist for $\varepsilon>0$ with $i$ fixed, it is natural to ask whether for $\varepsilon>0$, fixed and small, we can find an $i_{H}(\varepsilon)$, close to $i_{H}$, at which a homoclinic trajectory exists. There is a natural transversality condition that guarantees this property. The family of fixed points $x_{m}(i)$ admits a three-dimensional stable manifold and a two-dimensional unstable one. Their intersection is a homoclinic trajectory. In the singular limit, following the unstable manifold of $x_{m}(i)$ leads back to $S_{l}$ after two jumps. Extending $\mathcal{C}_{0}$ to include $i$, $x_{m}(i)$ is a (normally hyperbolic) invariant set in it, with two-dimensional invariant manifolds. The continuation of the unstable one using the singular flow, after two jumps intersects the stable manifold in the plane $i=i_{H}$. If this intersection is transverse then it persists for small $\varepsilon$ and $i$ close to $i_{H}$. We show this in appendix A adapting the arguments used in [24] to prove existence of relaxation oscillations.

To conclude this section, fig. 6 shows the bifurcation diagram of the whole system eq. 1 computed with AUTO-07p [1] for parameter values reported in appendix B. The numerics confirms the presence of a family of limit cycle (red curves) and its coexistence with a family of fixed points (blue curve). The family of periodic solutions terminates in a homoclinic trajectory for low values of $i$ (the numerical continuation was stopped at period $T=10^{4}$).

5 A common geometric picture

The bifurcation diagram illustrated in the previous section is understandably only one among many possible scenarios compatible with the three-dimensional geometry of fig. 2. While a detailed study of all possible cases is beyond the scope of this work, we wish to highlight how different types of bistability could share the same geometric structure. To do this we use ideas and techniques from [20]. As in section 3 we identify fixed points with the $i$-$v$ curve

[TABLE]

and divide them in three families $\mathcal{X}_{l}$, $\mathcal{X}_{m}$ and $\mathcal{X}_{h}$, separated by two folds. As noted in section 3 there is a value of current $i_{c}$, between the two folds, at which $x_{m}$ crosses $F_{l}$ to enter the unstable branch $M$. The scenario studied in section 4 assumes $i_{c}<i_{H}$ since the homoclinic bifurcation occurs when $x_{m}\in S_{l}$.

As a first variation we consider what happens when the bistable range extends to current values for which $x_{m}\in M$. The bifurcation $x_{m}\in F_{l}$ corresponds to a folded saddle-node. Beyond this bifurcation $x_{m}\in M$ is a node of the reduced dynamics while $x_{f}$ is a saddle. In this case the analysis is easily adapted from section 4. One must simply substitute the stable manifold of $x_{m}$ with the one of $x_{f}$, and use $\Pi_{f}(x_{f})$ in place of $x_{2}=\Pi_{f}(x_{1})$. Figure 7 shows the corresponding geometric construction. A classical example where this scenario occurs is the Hodgkin-Huxley model with the reversal potential of potassium increased. This situation of bistability has been studied in the early work [21]. Its planar reduction leads to the transcritical model [10]. Also in this case the boundary of bistability is a singular homoclinic trajectory. This trajectory, however, has to go through the folded singularity $x_{f}$ to reach $x_{m}$ on the unstable branch $M$.

Both cases discussed so far assume that $x_{m}$ and $x_{h}$ collide in a fold on $M$. Yet another scenario corresponds to this fold occurring on $S_{l}$, after $x_{h}$ crosses $F_{l}$. Also this crossing leads to a folded saddle-node, after which $x_{h}\in S_{l}$ can perturb to a stable fixed point. Local analysis around folded saddle-node shows the possibility of Hopf bifurcations [17], which are indeed found numerically. After this the system presents two stable fixed points. The relevant part of the reduced dynamics in this case is shown in fig. 7: the stable manifold of $x_{m}$ acts as separatrix between the basins of attraction of the two stable fixed points, while the one of $x_{f}$ (a folded saddle) separates initial conditions that reach a jump point on $F_{l}$ from those that remain on the critical manifold.

The examples above suggest that many possible variants for transitions between monostability and bistability are possible. We also note that many of the geometric constructions used in [7, 9, 20] have an analog in our setting, allowing, for example, non-plateau oscillations, contrary to the case showed in fig. 1. This flexibility is interesting in the perspective of connecting the present approach to the classification of bursting types according to the transitions that occur from rest to spike and vice-versa (see e.g. [15]).

6 Connections with phase portrait analysis

We close this paper by clarifying the connection between the proposed three-dimensional model and two published slow-fast phase portraits of rest-spike bistability.

The first phase portrait goes back to the seminal work of Hindmarsh and Rose [12, 13]. In one of the earliest attempts to model slow spiking and bursting, Hindmarsh and Rose proposed to modify FitzHugh-Nagumo model with a recovery variable that has a nonmonotonic activation function. Geometrically, this situation corresponds to a degenerate case of the planar pictures described in section 3 and section 4, in which all essential elements are contained on a line. As a result, the main elements of the three dimensional dynamics can be captured by constraining it to a plane, resulting in a simplified two-dimensional model of rest-spike bistability. This is characterized by the classical N-shaped critical manifold, as shown in fig. 8. The price paid for this simplification is that the flexibility of the two-dimensional slow dynamics described in section 5 is lost. For instance, bistability is only possible if $x_{l}$ lies out of the stripe delimited by $P_{l}$ and $F_{l}$, ruling out patterns in which the voltage of the resting state is between maximum and minimum of the spike. We note that the nonmonotonicity of the activation function in Hindmarsh-Rose model has the natural interpretation of summarising in one variable the distinct roles of an inward and an outward slow current.

The second rest-spike bistable phase portrait is the transcritical model of [7]. This model was obtained as a two-dimensional reduction of a conductance-based model that adds a slow calcium current to the Hodgkin-Huxley model [2]. The analysis of [7] rests on the presence of a transcritical bifurcation of the critical manifold. This bifurcation also directly relates to the mixed role of the slow variable as a source of both positive and negative feedback in the slow time-scale. A main motivation of the present paper was to understand the geometric picture generated by this motif in conductance-base models, where these two roles are often played by distinct variables.

To connect the transcritical bifurcation of the planar model [7] to the three-dimensional geometry of the present paper we consider how this planar reduction can be obtained. Referring to our model eq. 1 for simplicity, a planar reduction is typically obtained imposing an algebraic constraint between $n$ and $p$, which can be interpreted as a path $n(s)$, $p(s)$ [21]. After obtaining a dynamic equation for $s$ from a combination of $\dot{n}$ and $\dot{p}$, the system becomes

[TABLE]

which is a slow-fast planar model. Its critical manifold is given by

[TABLE]

It corresponds to the intersection of the critical manifold of the larger system with the surface

[TABLE]

A transcritical bifurcation is obtained when

[TABLE]

Geometrically this corresponds to a point at which the surface eq. 19 is tangent to the line of folds of the critical manifold, as shown in fig. 8. Similar geometric constructions lead to the presence of a transcritical bifurcation when reducing the Hodgkin-Huxley model with increased potassium reversal potential, as well as when reducing the same model augmented with a calcium current, as done in [7].

An equivalent interpretation of how the transcritical bifurcation arises is that the path $(n(s),\,p(s))$ defining the surface eq. 19 is tangent to the line of folds

[TABLE]

projected onto the $n$-$p$ plane. This is the simplest example of how singularities in the sense of [11] can be generated from elementary catastrophes, the core idea in the path formulation of [11, Ch.3 §12]. This is particularly interesting in view of [8], where singularity theory is used to obtain a global description of the critical manifolds of slow-fast planar systems relevant to neuronal dynamics. Two singularities play a prominent role: hysteresis, in connection with spiking, and winged cusp, for rest-spike bistability. Both these singularities can be realized as paths in the unfolding of the cusp catastrophe [11]. Interestingly, this bifurcation is often found in the fast subsystem of neuronal models (an early example being [26]), and it is typically related to the appearance and disappearance of bistability. For example, decreasing the sodium conductance in the Hodgkin-Huxley model leads to the appearance of this bifurcation, and the same is achieved by reducing $g_{m}$ in eq. 1. The presence of this type of bifurcation in these models suggests that those singularities can arise from model reduction similarly to what happens in the transcritical case.

7 Conclusions

We studied a simplified slow-fast model of neuronal activity that exhibits rest-spike bistability. The simplest physiological models of excitability include a fast-activating inward current and a slowly-activating outward current. Our model adds a slowly-activating inward current to this basic motif. We think of this model as a core structure for the generation of multistability in more general and realistic conductance-based models. We speculate that similar results are possible using a slowly inactivating outward current, which would have the same functional role of a slow positive feedback.

Through geometric singular perturbation theory we could analyze the geometry of this three-dimensional model. This geometry is rather simple, with the slow dynamics taking place on a classical N-shaped critical manifold. The saddle point on the critical manifold is a key feature of the proposed model. Its stable manifold acts as separatrix, while its unstable manifold determines whether multiple attractors are present. Moreover, a same geometric picture captures different types of bistability, suggesting a common framework to study different phenomena important to neuronal dynamics.

This is by no means the first study of a slow-fast systems with one fast and two slow variables, nor the first single-cell model of bistability. The value of this model is in that it explains how bistability can arise in a physiologically relevant context using a mechanism that is generic but not widely acknowledged. Our hope is that it contributes to the view that a combination of positive and negative feedback in the slow time-scale is a core element in the generation of neuronal patterns.

Appendix A Existence of a homoclinic trajectory

In this section we show that under the assumption of transversality, the intersection of stable and unstable manifolds that leads to a singular homoclinic trajectory persists for $\varepsilon>0$. We do this using the setting of [24] and in particular their results on maps defined by the flow of eq. 1. We recast these results in the notation of section 4 and refer the reader to the original work for details.

As in section 4, $i_{H}$ is the value of $i$ at which a singular homoclinic trajectory exists. We consider the reduced dynamics for this value of $i$, and fix a point $x_{u}$ on the unstable manifold of $x_{m}$ between $x_{m}$ and $F_{l}$. Similarly, we fix a point $x_{s}$ on the stable manifold between $P_{l}$ and $x_{m}$.

After a local change of coordinates we can find two neighborhoods of these points, $N_{s}$ and $N_{u}$, such that the critical manifold $\mathcal{C}_{0}$ corresponds to the plane $v=0$. The intersections of these neighborhoods with the planes $n=n_{s}$ and $n=n_{u}$ determine two surfaces $\Sigma_{s}$ and $\Sigma_{u}$. Rotating $n$ and $p$ if necessary, we can assume that $\Sigma_{u}\cap\mathcal{C}_{0}$ intersects the unstable manifold of $x_{m}$ transversally and only at $x_{u}$, and similarly for $\Sigma_{s}\cap\mathcal{C}_{0}$. For fixed $\delta$ we let $N_{\delta}=(i_{H}-\delta,\,i_{H}+\delta)$ and consider $\Sigma_{s}\times N_{\delta}$ and $\Sigma_{u}\times N_{\delta}$. If $\delta$ is small enough, stable and unstable manifolds of $x_{m}(i)$ intersect transversally these extended neighborhoods (in the critical manifold extended to include $i$). In the following we assume that $N_{s}$, $N_{u}$, $N_{\delta}$ are shrunk whenever necessary.

In section 3, we have characterized the stable manifold of $x_{m}$ for small $\varepsilon>0$, this is composed of a line on $\mathcal{C}_{\varepsilon}$ and the fibers based on it. In the limit $\varepsilon\to 0$, the singular stable manifold intersect $\Sigma_{s}$ transversally along one of these fibers. Thus if $\varepsilon$ and $\delta$ are small enough the same will be true for the stable manifold of $x_{m}(i)$ for fixed $i$ and $\varepsilon$. Moreover, since at $\varepsilon=0$, $i=i_{H}$ this intersection is a line of constant $p$, we can find a parametrization of it that has the form $p=p_{s}(v,i,\varepsilon)$. Similarly, the intersection of $\Sigma_{u}$ with the unstable manifold of $x_{m}(i,\varepsilon)$ defines two functions $v_{u}(i,\varepsilon)$ and $p_{u}(i,\varepsilon)$.

Notice that in this section we use $v$ and $p$ to parametrize the two slices $\Sigma_{s}$ and $\Sigma_{u}$, so that $v$ preserves its nature of fast variable. This differs from the use of $v$ and $p$ to parametrize the critical manifold as done in section 3 and section 4.

We can now use the same construction of [24] to obtain a map $\Pi:\Sigma_{u}\to\Sigma_{s}$ corresponding to the action of the flow. This has the form

[TABLE]

$R$ is exponentially small in $\varepsilon$ ($\lvert R\rvert+\lVert\nabla R\rVert<\exp(-c/\varepsilon)$) and in particular verifies

[TABLE]

$G$ has the form

[TABLE]

where $G_{0}:\Sigma_{u}\cap\mathcal{C}_{0}\to\Sigma_{s}\cap\mathcal{C}_{0}$ is the map defined by the singular flow. Smooth dependence on $i$ follows from standard results.

The only difference between this map and the Poincaré map defined in [24] is that we consider two different sections $\Sigma_{s}$ and $\Sigma_{u}$ rather than one.

Applying this map to $(v_{u},\,p_{u})$ we obtain the intersection of the unstable manifold of $x_{m}$ with $\Sigma_{s}$

[TABLE]

In this setting an intersection of stable and unstable manifolds corresponds to a solutions of

[TABLE]

where $v_{u}=v_{u}(i,\varepsilon)$ and $p_{u}=p_{u}(i,\varepsilon)$. Thus, we can define

[TABLE]

and a homoclinic trajectory corresponds to $P_{u}-P_{s}=0$.

At $\varepsilon=0$

[TABLE]

and the existence of the singular homoclinic trajectory at $i=i_{H}$ means that

[TABLE]

Assuming that

[TABLE]

an application of the implicit function theorem222We need a relaxed version of the implicit function theorem since the dependence on $\varepsilon$ is not smooth but only continuous, this is proven in [18]. guarantees the existence of a continuous functions $i_{H}(\varepsilon)$ such that $P_{u}(i_{H}(\varepsilon),\varepsilon)=P_{s}(i_{H}(\varepsilon),\varepsilon)$.

At $\varepsilon=0$, $P_{u}(i,0)$ is the intersection of the singular unstable manifold (after two jumps) with $\Sigma_{s}\cap\mathcal{C}_{0}$, while $P_{s}(i,0)$ corresponds to the intersection of the stable manifold of the reduced flow and $\Sigma_{s}\cap\mathcal{C}_{0}$. Condition eq. 30 corresponds to transversality of the intersection between the manifolds $p=P_{s}(i,0)$ and $p=P_{u}(i,0)$ in the extended neighborhood $\Sigma_{s}\times N_{\delta}$. Since the invariant manifolds of $x_{m}(i)$ can be obtained applying the singular flow to these two sections, we see that condition eq. 30 is equivalent to transversality of the intersection between the invariant manifolds of $x_{m}(i)$ on the critical manifold (where the unstable manifold has been continued past two jumps using the singular flow).

Appendix B Parameters

The analysis in sections 3 and 4 uses the following numerical values for the parameters of eq. 1

[TABLE]

Parameters related to $m$, $n$ and leak are approximated from those in [4] for Morris-Lecar model in the case Hopf, after normalizing $v$ between -1 and 1. Other parameters have been chosen to obtain the desired bistability.