Excitable neurons, firing threshold manifolds and canards

Abstract

We investigate firing threshold manifolds in a mathematical model of an excitable neuron. The model analyzed investigates the phenomenon of post-inhibitory rebound spiking due to propofol anesthesia and is adapted from McCarthy et al. (SIAM J. Appl. Dyn. Syst. 11(4):1674-1697, 2012). Propofol modulates the decay time-scale of an inhibitory GABAa synaptic current. Interestingly, this system gives rise to rebound spiking within a specific range of propofol doses. Using techniques from geometric singular perturbation theory, we identify geometric structures, known as canards of folded saddle-type, which form the firing threshold manifolds. We find that the position and orientation of the canard separatrix is propofol dependent. Thus, the speeds of relevant slow synaptic processes are encoded within this geometric structure. We show that this behavior cannot be understood using a static, inhibitory current step protocol, which can provide a single threshold for rebound spiking but cannot explain the observed cessation of spiking for higher propofol doses. We then compare the analyses of dynamic and static synaptic inhibition, showing how the firing threshold manifolds of each relate, and why a current step approach is unable to fully capture the behavior of this model.

Fig. 1

The inhibitory gating variable, s and the inhibitory current, $I_{\mathrm{syn}}$. The dimensionless modified propofol model is simulated with inhibition modeled at $t=100$. Here ${\tau }_s=10$. a A time-trace of s. The dynamics of the gating variable s are decoupled from the rest of the system. The gating variable decays exponentially with decay constant ${\tau }_s$ where the value of this decay constant is linked to the amount of propofol present. b A time-trace of $I_{\mathrm{syn}}$. Recall the synaptic current is given by $I_{\mathrm{syn}}=g_{i}s(V-E_i)$. At the onset of dynamic behavior the synaptic current exhibits a rapid spike followed by an approximately exponential decay. Note the small kink at approximately $t=180$ due to the spike (action potential) in v

Fig. 2

Rebound spiking with dynamic inhibition. The dimensionless modified propofol model is simulated for four different values of ${\tau }_s$, 7, 8, 21, and 22. a${\tau }_s=7$, b${\tau }_s=8$, c${\tau }_s=21$, d${\tau }_s=22$. Voltage time-traces are shown within each panel, and the corresponding time-traces of s are shown below. Note that there exists a range of ${\tau }_s$ values for which rebound spiking occurs, specifically for ${\tau }_s\in [8,21]$

Fig. 3

Rebound spiking with step inhibition. A classical current step protocol is applied to the dimensionless modified propofol model. Here, the value of $I_{\mathrm{syn}}$ is held constant at $3.5\times {10}^{-4}$, roughly mimicking the slowly decaying portion of the dynamic inhibition current. The protocol is simulated for four different durations, 10, 20, 50, and 200 ms. a Duration = 10 ms, b Duration = 20 ms, c Duration = 50 ms, d Duration = 200 ms. Voltage time-traces are shown within each panel, and the corresponding time-traces of $-I_{\mathrm{syn}}$ are shown below. As opposed to dynamic inhibition, here we observe no cessation of spiking for prolonged inhibitory input. On the contrary, we observe an isolated couplet of spikes for prolonged inhibition, and as the duration of the inhibitory step is further increased, we observe a triplet of spikes

Fig. 4

Critical manifold in $(v,w,s)$-space. The critical manifold, $S_0$, is the set of stationary equilibria within the layer problem, (7). This surface defines the interface between fast and slow dynamics: the layer problem describes rapid movement toward or away from $S_0$, the reduced problem describes slower dynamics on the manifold itself. Within a physiologically relevant domain, $S_0$ has a roughly cubic structure, consists of upper and lower fold curves, $F^{+}$ and $F^{-}$ respectively, and three sheets of stationary equilibria. There exists a single lower attracting sheet, $S_a^{-}$, and upper and lower repelling sheets, $S_r^{+}$ and $S_r^{-}$, respectively. On the uppermost surface lies a line of Andronov–Hopf bifurcations from which a family of limit cycles emanate. The line of Andronov–Hopf bifurcations lies out of range while the family of limit cycles terminate within (not shown). Note the shape and stability properties of the critical manifold are independent of ${\tau }_s$

Fig. 5

Folded saddle geometry and associated trajectories. The geometry of a generic folded saddle. The folded saddle (purple) is denoted ${\mathbf{p}}_{\mathrm{fs}}$, while the fold curve (black dashed) is denoted F. a Folded saddle geometry according to the singular reduced problem. The folded singularity resembles an ordinary saddle equilibrium with stable and unstable invariant manifolds (red). The trajectories (blue) follow these invariant manifolds moving away from the stable manifold and toward the unstable. b Within the desingularized reduced problem the dynamics on the repelling surface, $S^r$, are reversed due to the rescaling of time (desingularization). Trajectories may pass through the folded saddle with non-zero velocity traveling either of the invariant manifolds. These trajectories correspond to singular canards; the stable invariant manifold to the true canard and the unstable invariant manifold to the faux canard. c Folded saddle geometry in 3D. The true canard acts as a separatrix on the attracting surface, $S^a$. If a trajectory lands within the region enclosed by the true and faux canards, then it is bound away from the fold curve. However, if the trajectory lands within the region enclosed by the fold curve and true canard it travels toward the fold curve. Here, the trajectory “jumps off” due to a blow up in finite time of the desingularized reduced problem, where subsequent dynamics are dictated by the layer problem. The region within which trajectories necessarily ‘jump off’ the critical manifold is indicated (purple shaded)

Fig. 6

The desingularised reduced problem. The desingularized reduced problem, (15), near the lower fold is projected onto $(v,s)$-space. Here, ${\tau }_s=15$, although the basic structure shown here is common to all values of ${\tau }_s$ analyzed. The lower fold curve, $F^{-}$ (gray dashed), is indicated. The folded saddle, ${\mathbf{p}}_{\mathrm{fs}}$ (purple), lies on the fold curve and gives rise to two invariant manifolds: a stable invariant manifold, $W^s({\mathbf{p}}_{\mathrm{fs}})$ (red), and an unstable invariant manifold, $W^u({\mathbf{p}}_{\mathrm{fs}})$ (blue). The dynamics in $(v,s)$-space above $F^{-}$ are reversed, and so the stable and unstable manifolds have reversed stability properties above $F^{-}$. Arrows indicate motion along the invariant manifolds. Each manifold terminates at a stable node equilibrium, within the reduced flow. The stable manifold terminates at ${\mathbf{eq}}_2$ on $S_a^{-}$ (orange) and the unstable manifold terminates at ${\mathbf{eq}}_3$ on $S_r^{-}$ (green). Note, if a singular trajectory lands onto the shaded region of $S_0$, it eventually undergoes a rebound spike. Inset: A magnification of the desingularized reduced problem near the folded saddle

Fig. 7

Singular global trajectories in $(v,w,s)$-space. The desingularized reduced problem is projected onto the critical manifold near the lower fold (gray dashed) for ${\tau }_s=15$. Singular solution trajectories from the layer problem, (7), and the reduced problem, (15), are concatenated to produce singular global trajectories (black). a From the initial condition (IC), the layer problem dictates that the singular trajectory falls onto $S_a^{-}$. Since the singular trajectory base point lies within the region bound between the canard separatrix and $F^{-}$, the reduced problem dictates that the trajectory evolves toward the fold curve. At the fold, due to a singular blow up of the reduced problem, the trajectory undergoes fast oscillations within the layer problem. This trajectory corresponds to a successful post-inhibitory rebound spike. b The corresponding system projected onto $(v,s)$-space. c The three-dimensional system from a different angle. Note the initial approach of the trajectory onto the critical manifold. d The corresponding system projected onto $(v,w)$-space. Note, this view provides a clear delineation between the singular canard and singular trajectories while the others do not. Hereafter, this projection is used when comparing canards and their respective trajectories

Fig. 8

Singular limit spiking activity transitions. Singular solution trajectories from the layer problem, (7), and the reduced problem, (15), are concatenated to produce singular global trajectories (black). The singular limit predicts a range of ${\tau }_s$ values for which rebound spiking occurs; ${\tau }_s\in [5,24]$. The layer problem dictates that the trajectory has a base point on $S_a^{-}$ independent of ${\tau }_s$. Once on the manifold, the reduced problem dictates that the trajectory remains to one side of the canard separatrix. a The singular trajectory for ${\tau }_s=24$ in $(v,w)$-space. Since this trajectory lies to the right of the separatrix, it evolves in time toward the fold curve, $F^{-}$ (gray dashed), at which point the layer problem describes the onset of oscillatory behavior. This singular prediction corresponds to a successful rebound spike. b The singular trajectory for ${\tau }_s=25$ in $(v,w)$-space. This trajectory lies to the left of the separatrix and evolves in time toward ${\mathbf{eq}}_3$ (green). This singular prediction corresponds to an unsuccessful rebound spike. An animation of this figure under variation of ${\tau }_s$ is given within Additional file 1

Fig. 9

Singular limit prediction of spiking. Within the singular limit, the position of the base point determines whether a singular trajectory results in a rebound spike. We plot the s-value of the trajectory initial condition (black dashed), which coincides with the base point on the critical manifold. Note that the initial condition, and thus the base point is independent of ${\tau }_s$. For integer values of ${\tau }_s$, the s-value at which the corresponding canard separatrix crosses the v-value of the base point (red/blue points) is plotted. Thus, if an initial condition lies at a higher s-value than the corresponding canard separatrix (red), the trajectory falls onto the region of $S_a^{-}$, for which a trajectory rebound spikes. Otherwise, if the initial condition lies at a lower s-value than the corresponding canard separatrix (blue), the trajectory does not go on to spike. A linear interpolation shows the roughly parabolic shape, which gives rise to two rebound spiking transitions, and thus a single range of ${\tau }_s$ for which there exists rebound spiking. An animation of this figure under variation of ${\tau }_s$ is given within Additional file 1

Fig. 10

Nonsingular canards and global trajectories. The true nonsingular canard (blue) for the three-dimensional modified propofol model is located for ${\tau }_s=15$. This canard is then continued in ${\tau }_s$ using a suitable boundary condition problem. We thus obtain the true canard for any value of ${\tau }_s$. This true canard is compared to the corresponding trajectory (black) of the three-dimensional modified propofol model. Here, we see that the true nonsingular canard correctly bounds the trajectories into spiking and quiescent regimes. The folded singularity (purple) and ${\mathbf{eq}}_3$ (green) are indicated. a The nonsingular trajectory and associated canard for ${\tau }_s=8$ in $(v,w)$-space. b The nonsingular trajectory and associated canard for ${\tau }_s=9$ in $(v,w)$-space

Fig. 11

Singular limit analysis of step inhibition. A classical current step protocol is applied to the dimensionless modified propofol model. A geometric singular perturbation analysis is then performed on this system. The synaptic current is held at $I_{\mathrm{syn}}=3.5\times {10}^{-4}$ to simulate inhibition. We compare the critical manifold for each value of $I_{\mathrm{syn}}$. On each manifold, the lower branch (magenta) is linearly stable, the upper branch (gray) linearly unstable. The stable node equilibria (green), ${\mathbf{eq}}_3$ and ${\overline{\mathbf{eq}}}_3$, and saddle equilibrium (orange), ${\mathbf{eq}}_3$ are indicated. Initially, the system starts at rest on the node equilibrium of the $I_{\mathrm{syn}}=0$ critical manifold. At the onset of inhibition, the critical manifold is shifted in the direction of negative w. In the singular limit, the layer problem dictates that the trajectory falls to the lower branch of the critical manifold. Once on the manifold, the reduced problem dictates that the trajectory slowly approaches ${\overline{\mathbf{eq}}}_3$. Once the inhibitory current is removed, the manifold shifts back to its original position, at which point the layer problem determines that the trajectory shoots upward in the direction of positive v. In the singular limit, the threshold manifold (blue dashed) is the concatenation of the middle branch of the critical manifold with the fast fibre through the lower fold of the critical manifold. If the singular limit trajectory passes this manifold, the singular limit predicts a spike event

Fig. 12

Non-singular analysis of step inhibition. A classical current step protocol is applied to the dimensionless modified propofol model for 40 ms, and 250 ms. We compare the critical manifolds (v-nullclines) for $I_{\mathrm{syn}}=0$ and $I_{\mathrm{syn}}=3.5\times {10}^{-4}$. On each manifold, the lower branch (magenta) is linearly stable, the upper branch (gray) linearly unstable. The corresponding stable node equilibria (green), ${\mathbf{eq}}_3$ and ${\overline{\mathbf{eq}}}_3$, and saddle equilibrium (orange), ${\mathbf{eq}}_2$, are indicated. Initially, the system starts at rest at the node equilibrium of the $I_{\mathrm{syn}}=0$ critical manifold. At the onset of inhibition, the manifold is shifted in the direction of negative w. The trajectory falls to the lower branch of the shifted manifold, slowly approaching ${\overline{\mathbf{eq}}}_3$. Once the inhibitory current is removed, the critical manifold shifts back to its original position, at which point the trajectory rapidly shoots in the direction of positive v. The threshold manifold for trajectory spiking is indicated (orange dashed). a If the trajectory passes this separatrix, when released from inhibition, the system spikes. The trajectory is reset with a net shift in the direction of positive w. Upon resetting, the trajectory lies to the left of the separatrix and falls within the basin of attraction of the stable node equilibrium. b As the length of inhibition is increased, the trajectory evolves closer toward the shifted equilibrium position. Here, the trajectory undergoes three spikes before being reset to the left of the separatrix and finally coming to rest. Note the system used here to calculate the threshold manifold makes use of the reduction $\tilde{x}={\tilde{x}}_{\mathrm{\infty }}(v)$, for $x=m,h,n$, in order to allow calculation along the repelling middle branch. Locally (for sub and perithreshold regimes) this reduced system well approximates the modified propofol model, despite being unable to repolarize after a spiking event

Fig. 13

Firing threshold manifolds: dynamic vs. static inhibition. Firing threshold manifolds of dynamic inhibition and static inhibition are here plotted together. Recall the firing threshold manifolds for dynamic inhibition are the nonsingular canards, here shown for odd values of ${\tau }_s$ between 1–15 (solid purple). The firing threshold manifold for static inhibition is the stable manifold of ${\mathbf{eq}}_2$ (orange dashed) shown earlier in Fig. 12. For decreasing values of ${\tau }_s$ (darker shades of purple), the nonsingular canards more closely resemble the static inhibition threshold manifold. Note each system (dynamic or static inhibition) used to calculate a threshold manifold required the reduction $\tilde{x}={\tilde{x}}_{\mathrm{\infty }}(v)$, for $x=m,h,n$, to allow calculation along repelling manifolds

Fig. 14

Singular and nonsingular trajectories of the original and modified propofol model. Singular global trajectories (black, solid) are here overlaid with their corresponding nonsingular propofol model trajectories (black, dashed). Here, we compare the original and modified propofol model solution trajectories and singular limit predictions for ${\tau }_s=20\text{ ms}$ (spiking occurs in both models for this value). a Within the original propofol model, solution trajectories are kept in close proximity to the lower fold curve. Near the fold, normal hyperbolicity is lost and the time-scale separation between fast and slow variables breaks down. Hence, the singular trajectory does not accurately describe the fast approach onto the critical manifold. This results in a substantial deviation of solution trajectories from the singular limit prediction of subsequent slow dynamics; note, for example, the discrepancy between the s-values at which the non-singular and singular trajectories spike. b The modified propofol model solution trajectories are no longer constrained near the fold curve, and thus the time-scale separation is preserved. This results in a more accurate singular limit prediction of nonsingular trajectories

Fig. 15

Gating variable steady-state and time-scale functions. a The four gating variable steady-state functions, $m_{\mathrm{\infty }}(V)$ (red), $h_{\mathrm{\infty }}(V)$ (orange), $n_{\mathrm{\infty }}(V)$ (green), and $w_{\mathrm{\infty }}(V)$ (blue). These functions give the opening probability of their respective ion channels for a given value of the membrane voltage, V. b The four gating variable time-scale functions, ${\tau }_m(V)$ (red), ${\tau }_h(V)$ (orange), ${\tau }_n(V)$ (green), and ${\tau }_w(V)$ (blue). These functions describe the rough time-scales on which each gating variable evolves for a given value of the membrane voltage. c A closer look at panel b reveals that ${\tau }_m(V)$, ${\tau }_h(V)$ and ${\tau }_n(V)$ evolve on much smaller time-scales than ${\tau }_w(V)$

Fig. 16

Variable time-scales. The time derivatives of each variable are plotted over a spike event; v (red), m (orange), h (green), n (blue), w (purple), and s (gray). We find that there is a clear order of magnitude difference between the time derivatives of the variables v, m, h, n, and w, s. Inset: A magnification of the derivative time-traces; here the relative magnitudes of $\frac{dw}{dt}$ and $\frac{ds}{dt}$ can be seen