Skip to main content

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.

1 Introduction

Excitable neurons are typically at rest, but can fire action potentials in response to certain forms of stimulation. Many textbooks refer to excitability as an all-or-none response, i.e. a “subthreshold” synaptic input leads to a small graded postsynaptic potential while a “superthreshold” input evokes an action potential. One then seeks to find an action potential threshold, i.e., a particular voltage value that demarcates the all-or-none response. However, this is a misconception. Using geometric analysis of neural models, it was FitzHugh [10] who first noticed that a firing threshold, if it exists, is never a number but a manifold.

Another characteristic feature of neurons is the existence of processes that evolve on multiple time-scales. The interaction of ionic currents acting on different time-scales is responsible for the creation of action potentials in neurons. This time-scale feature leads to mathematical models of neurons often referred to as singular perturbation problems. These models are particularly amenable to analysis using geometric singular perturbation theory (GSPT) [9, 14] with the specific aim of giving predictions of model dynamics based on singular limit observations. Threshold phenomena are closely related to folded critical manifolds of such singularly perturbed neural problems. In the famous 2-dimensional singularly perturbed FitzHugh–Nagumo model, it is the repelling middle branch of the folded (cubic-shaped) critical manifold that forms the firing threshold manifold in the singular limit. In the full system, this firing threshold manifold perturbs to a nearby repelling slow manifold.

In higher dimensional models with more than one slow variable, the geometry of such singularly perturbed problems becomes quite intricate and folded singularities play a prominent role. Folded singularities lie on a fold of the critical manifold where stable and unstable branches of this (higher dimensional) manifold meet. Canards [2, 3, 8, 15, 21, 23, 24] are trajectories associated with folded singularities and connect the stable and unstable branches of the critical manifold. These special solutions have been identified as important objects in explaining complex oscillatory patterns known as mixed-mode oscillations (MMOs) [4]. There now exists a substantial amount of literature on applications of canard theory, and we refer the interested reader to detailed tables on relevant literature provided in a review on MMOs [5].

Canards form boundaries of different dynamic behavior—they are separatrices by nature. More importantly, they encode slow time-scales, i.e., temporal information is reflected in the geometry of a canard. In the two-dimensional case, Izhikevich [13] clearly highlights the fact that canard trajectories and repelling slow manifolds provide the best approximation to the firing threshold manifold, hence giving a mathematical “structure” to the famous No Man’s Land by R. FitzHugh [11]. More recently, Desroches et al. [6] discuss the relationship between canards and excitability thresholds in planar slow-fast systems by identifying inflection lines of the flow.

In the present study, we focus on a higher dimensional neural model and the specific role of folded saddle canards as firing threshold manifolds. We observe that varying one of the slow time-scales changes the boundaries of different dynamic behaviors. The slow time-scale is thus encoded in the position of the saddle canard separatrix, and so, remarkably, a change in a slow time-scale can be “seen” in a geometric object. Folded saddle canards have been shown to form firing threshold manifolds in general Morris–Lecar/FitzHugh–Nagumo type neural models with dynamic current input [25] and even in a climate model of the so-called “compost-bomb instability” by Wieczorek et al. [26]. In both, under variation of a particular parameter, the position of a folded saddle canard varies explaining the excitability properties of the model.

This geometric observation becomes important when we try to understand changes in neural dynamics. Neurons constantly sense their environment, e.g., temperature, acidity or glucose, and they encode fluctuations of these parameters by changing their ion channel dynamics. The famous Hodgkin–Huxley model of the squid giant axon [12] incorporates temperature changes through a Q 10 -temperature factor that increases the speed of the ion channel gates with increasing temperature. If such environmental changes are encoded in different speeds of slow ion channels, then we expect this to be reflected in a change of the firing threshold manifold of the neuron. Thus, identifying the cause of different neural dynamics through the specific position of a canard in a singularly perturbed system becomes a valuable diagnostic tool to understand this phenomenon, and it is the focus of this work.

2 Propofol and Rebound Spiking

As a case study, we investigate the role of the general anesthetic propofol on rebound spiking in the central nervous system [17]. Many general anesthetics, including propofol, prolong the duration of GABAergic inhibitory postsynaptic currents (IPSCs), and this action contributes to the behavioral properties of these drugs [1]. Mathematically speaking, propofol changes the slow time-scale of the deactivation of the GABAa receptor channel. Paradoxically, low doses of propofol causes excitation rather than sedation. This behavior can already be observed in an isolated single cell model that receives GABAergic IPSCs [18]. We adapt this propofol neuron model formulated in [18] slightly in order to more clearly emphasize the role of folded saddle canards in the observed dynamic behavior. We note that only minimal adjustments have been made in order to preserve the qualitative behavior of the model, namely the observation of post-inhibitory rebound spiking for a window of GABAa synaptic time-scale values. The modification consists of two parameter changes, the details of which are given below. The essential difference is that this modification shifts the resting membrane potential to a lower, more hyperpolarized, voltage value, allowing a more uniform separation of time-scales over a range of GABAa time-scale values. This modification enables us to make full use of the machinery of GSPT. More details of the relation between the two models can be found in the last section of the paper.

The modified propofol model is given as a six-dimensional conductance based model with Hodgkin–Huxley type dynamics:

c m d V d t = I Na I K I L I M I syn + I app , d x d t = ( x ( V ) x ) τ x ( V ) , x = m , h , n , w , d s d t = s / τ s ,
(1)

where the ionic currents are defined as

I Na = g Na m 3 h ( V E Na ) , I K = g k n 4 ( V E k ) , I L = g L ( V E L ) , I M = g m w ( V E k ) , I syn = g i s ( V E i ) .
(2)

The membrane currents consist of a fast sodium current, I Na , a potassium current, I K , and a leak current, I L , collectively referred to as the spiking currents, and a slow muscarinic potassium M-current, I M . This model neuron receives GABAa IPSCs, here modeled by the current I syn . The applied current, I app , models tonic external drive.

As is typical in Hodgkin–Huxley type conductance based models, gating activation and deactivation variables model changes in ion channel conductances. Each take values restricted within the interval [0,1] reflecting opening probabilities. Here, m, n, and w are gating activation variables, whereas h is a gating deactivation variable. Similarly, the variable s is a gating variable for synaptic GABAa receptor channels. The value of s corresponds to the inhibitory input the neuron receives. In order to model inhibitory postsynaptic input in the absence of a network of neurons, the variable s is set to s 0 =0.714 at the instantaneous activation of the GABAa receptor channel. This s value is used as it most closely reflects the dynamic change in synaptic input within a network of coupled neurons. The deactivation of this channel is then simply modeled by (slow) exponential decay with time constant τ s ; see Fig. 1 for a comparison of the synaptic gating variable s and the synaptic inhibitory current I syn . Since propofol alters the deactivation time-scale of GABAa receptor channels, we consider τ s as our main study parameter.

Fig. 1
figure 1

The inhibitory gating variable, s and the inhibitory current, I syn . The dimensionless modified propofol model is simulated with inhibition modeled at t=100. Here τ 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 τ s where the value of this decay constant is linked to the amount of propofol present. b A time-trace of I syn . Recall the synaptic current is given by I 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

The explicit forms of the steady state functions, x (V), and the time-scale functions, τ x (V), are given in the Appendix (see Fig. 15 for plots). The model constants are given in Table 1. Note the two modifications from the original propofol model in [18]: g i =4 mS (originally g i =0.04 mS) and the w-nullcline is shifted by 3 mV in the direction of negative v.

Table 1 Propofol neuron and network model system parameters

2.1 Dynamic Inhibition, but not the Current Step Protocol, Leads to Cessation of Spiking with Increased Inhibition

Classically, an applied current step protocol has been used to understand rebound spiking [13, 19]. However, as was pointed out in [18], we find here that it is necessary to consider the dynamic interplay between the inhibitory synaptic current and the slow potassium M-current in order to explain the observed behavior of spiking within the modified propofol model (compare Figs. 2 and 3). Within the modified propofol model, there exists a specific range of values of τ s for which the system exhibits spiking behavior. Outside this range, the system fails to spike (Fig. 2). This matches well with the experimental observation that only low doses of propofol are able to create the state of paradoxical excitation [20].

Fig. 2
figure 2

Rebound spiking with dynamic inhibition. The dimensionless modified propofol model is simulated for four different values of τ s , 7, 8, 21, and 22. a τ s =7, b τ s =8, c τ s =21, d τ 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 τ s values for which rebound spiking occurs, specifically for τ s [8,21]

Fig. 3
figure 3

Rebound spiking with step inhibition. A classical current step protocol is applied to the dimensionless modified propofol model. Here, the value of I syn is held constant at 3.5× 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 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

Using a traditional approach, the spiking behavior is studied with a step protocol. Here, an applied current is switched on, kept at a constant level, then removed. By holding I syn constant, the dynamics in s are lost, thus rendering the step protocol system five-dimensional as opposed to the six-dimensional dynamic inhibition system (1). This step protocol is usually able to reproduce spiking patterns, while simplifying input dynamics, and thus give insight into the associated spiking behavior. However, applying a step protocol to the present model, we find a single transition from inactivity to isolated spiking as τ s is increased. As τ s is further increased a transition from single spiking to a doublet, and then to a triplet of spikes is observed (Fig. 3). A maximum of three spikes is observed despite further increases in duration of inhibition. The step protocol is unable to reproduce a cessation of spiking for increased synaptic inhibition. It thus becomes apparent that there is a necessary dynamic mechanism required to yield a specific range of spiking under variation of synaptic inhibition time-scale.

In the present study, our aim is to identify firing threshold manifolds of the propofol model, both with dynamic and with static inhibition, thus explaining the spiking behavior under variation of propofol dosage and duration of current step inhibition, respectively. Key to this aim is the use of GSPT for which a detailed analysis of time-scales is necessary. This time-scale analysis is presented in Sect. 3. Identifying multiple time-scales in system (1) implies a splitting of solution trajectories of (1) into segments of fast and slow dynamics. These fast and slow dynamics are captured by lower dimensional subsystems, termed the layer and reduced problems, respectively. GSPT uses these lower dimensional subsystems, studied in Sect. 4, to provide insight into the geometric structures which govern the behavior of the model, and thus to predict the dynamics of the full (higher dimensional) system (1). Results are given in Sect. 5 for the case of dynamic inhibition. In particular, we identify a singular canard of folded saddle type as the separatrix that forms the firing threshold manifold. In Sect. 6, a similar analysis is carried out for the case of static inhibition; again a singular limit prediction of the spiking threshold manifold is identified. In both the dynamic and static inhibition cases, a numerical confirmation of singular limit predictions is made by calculating the true firing threshold manifolds. This emphasizes the predictive power of GSPT. We also point out that the firing threshold manifold in the case of static inhibition is in fact the structure which the firing threshold manifold in the case of dynamic inhibition approaches in the limit as τ s 0, i.e., in the transition from slow to fast synaptic GABAa inhibition decay rates. In Sect. 7, we discuss the original and modified propofol models and their respective analyses. Finally, in Sect. 8, we make some concluding remarks about canards and excitability in neural models.

3 Time-Scales and Dimensional Analysis

By observation of the time traces in Figs. 2 and 3, it can be argued that there exists a multiple time-scales structure within the typical solution trajectories of the modified propofol model (1). In order to investigate this behavior further, a dimensional analysis of the system is performed so as to roughly determine the time-scales on which each variable evolves.

3.1 The Modified Propofol Model Has Three Distinct Time-Scales

As part of the nondimensionalization procedure, we introduce the three scaling constants, k v , k t , and g max . The membrane voltage is scaled according to V=v k v , where k v is the typical range over which the modified propofol model trajectories evolve in the V-direction. Time is scaled as t=τ k t , where k t has units of ms and is determined below. The ion channel conductances are scaled according to g x = g ¯ x g max , where g max =max{ g Na , g k , g L , g m , g i }. This leads to a dimensionless form of the modified propofol model (1),

c m k t g max d v d τ = I ¯ Na I ¯ K I ¯ L I ¯ M I ¯ syn + I ¯ app , τ x k t d x d τ = ( x ˜ ( v ) x ) τ ˜ x ( v ) , x = m , h , n , w , τ s k t d s d τ = s ,
(3)

where the ionic currents are now defined as

I ¯ Na = g ¯ Na m 3 h ( v e Na ) , I ¯ K = g ¯ k n 4 ( v e k ) , I ¯ L = g ¯ L ( v e L ) , I ¯ M = g ¯ m w ( v e k ) , I ¯ syn = g ¯ i s ( v e i ) .
(4)

The dimensionless steady state functions and time-scale functions are given as x ˜ (v)= x (v k v ) and τ ˜ x (v)= τ x (v k v )/ τ x , respectively. Note that the function τ x (v k v ) has been scaled by

τ x = min e K v e Na τ x (v k v )

in order to yield the dimensionless time-scale function τ ˜ x (v). Hence, the right-hand sides in system (3) are all of order O(1) and we can identify roughly the typical time-scales on the left-hand side:

τ v : = c m / g max = 0.01 ms , τ m = 0.03 ms ( 0.05 ms ) , τ h = 0.25 ms ( 0.5 ms ) , τ n = 0.3 ms ( 0.7 ms ) , τ w = 38 ms ( 46 ms ) , τ s [ 3 , 100 ] ms .
(5)

Note that the τ x , x=m,h,n,w values given in parentheses refer to the time-scale functions for subthreshold v0.5. The time-scale on which s evolves is given directly by the value of τ s .

It is now possible, based on the above numbers, to propose a hierarchy of variables according to the time-scales on which they evolve. We classify the variables v and m superfast, h and n fast, and w and s slow. Recall the value of τ s directly determines the time-scale on which s evolves. Accordingly, care should be taken in the following singular perturbation analysis since it may no longer be valid to consider s slow as τ s becomes sufficiently small. (In [18], s is considered a fast variable when τ s is small.)

Another way to check time-scales separation is to consider the maximum magnitude of the time derivative of each state variable over the course of a full spiking trajectory; see Fig. 16. For τ s =8, we observe that the max{| d x d t |} for x=v,m,h,n,w,s is roughly given by (10.2,6.4,3.3,2.1,0.025,0.089); compare these values with the inverse of the time-scales given in (5). This suggests that our proposed time-scale hierarchy with w and s slow is reasonable. As alluded to earlier, the specific insight of this model is the interplay between the slow M-current and the inhibitory synaptic current contributing to rebound spiking.

We set k v =100 mV as a typical voltage scale for a spiking neuron and k t =1 ms as a typical fast time-scale of the “fast” variables (n,h). Thus, system (3) becomes

δ d v d τ = I ¯ Na I ¯ K I ¯ L I ¯ M I ¯ syn + I ¯ app = : f 1 ( v , m , h , n , w , s ) , δ d m d τ = c m τ m g max ( m ˜ ( v ) m ) τ ˜ m ( v ) = : f 2 ( v , m ) , d h d τ = k t τ h ( h ˜ ( v ) h ) τ ˜ h ( v ) = : f 3 ( v , h ) , d n d τ = k t τ n ( n ˜ ( v ) n ) τ ˜ n ( v ) = : f 4 ( v , n ) , d w d τ = ε ( w ˜ ( v ) w ) τ ˜ w ( v ) = : ε g 1 ( v , w ) , d s d τ = ε ( τ w τ s s ) = : ε g 2 ( s ) ,
(6)

where δ:= c m /( k t g max )=0.011 measures the time-scale separation between the fast variables (h,n) and the superfast variables (v,m) while ε:= k t / τ w =0.0251 measures the time-scale separation between the fast variables (h,n) and the slow variables (w,s). Thus, the dimensionless system (6) is a singularly perturbed problem with singular perturbation parameters δ1 and ε1. This suggests an inherent three time-scales problem. The specific insight of the model and, therefore, the corresponding analysis should focus on the interplay of the slow M-current and the inhibitory synaptic current contributing to rebound spiking. We therefore group the fast and superfast variables together into one “fast” pool and consider system (6) a 4-fast/2-slow (v,m,h,n)/(w,s) problem. In doing so, we choose ε as the main singular perturbation parameter and keep δ fixed.

4 Geometric Singular Perturbation Theory

The identification of distinct fast and slow time-scales in the modified propofol model allows us to utilize the methods of GSPT to identify the key differences in geometric structure that account for the differences in spiking dynamics between static and dynamic inhibition. Thus, we proceed to identify two lower dimensional subsystems that govern the fast and slow dynamics in order to give us insight into the behavior of the full higher dimensional system.

4.1 Layer Problem

In the limit as ε0 of system (6), we derive the layer problem, which approximates the fast components (v,m,h,n) of the full system dynamics:

δ v ˙ = f 1 ( v , m , h , n , w , s ) , δ m ˙ = f 2 ( v , m ) , h ˙ = f 3 ( v , h ) , n ˙ = f 4 ( v , n ) , w ˙ = 0 , s ˙ = 0 ,
(7)

where x ˙ denotes the derivative of x with respect to the fast time-scale, t fast =τ. Note that the slow variables (w,s) are considered parameters in the layer problem. An important geometric object is the critical manifold, which is defined as the set of equilibria of the layer problem:

S 0 = { ( v , m , h , n , w , s ) R 6 | f 1 = f 2 = f 3 = f 4 = 0 } = { ( v , m , h , n , w , s ) R 6 | f 1 ( v , m ˜ ( v ) , h ˜ ( v ) , n ˜ ( v ) , w , s ) = 0 } .
(8)

Note here that f 2 = f 3 = f 4 =0 can be solved for x= x ˜ (v), x=m,h,n, yielding a single expression for the critical manifold in (v,w,s)-space. This allows for a projection of the critical manifold into three-dimensional space, shown in Fig. 4. Under variation of w and s, the two-dimensional critical manifold contains an upper and lower set of fold or saddle-node bifurcations, F + and F , respectively. These fold lines contain equilibria with a single zero eigenvalue and, in effect, partitions S 0 into three surfaces of equilibria. The lowermost surface of S 0 is denoted S a and is an attracting surface of equilibria, i.e., all equilibria have four eigenvalues with negative real parts. The surface bound between the two fold curves, S r , is a repelling surface of equilibria with a single eigenvalue with positive real part. Above the upper fold lies a second repelling surface of equilibria, S r + , containing equilibria with two eigenvalues with positive real part.

Fig. 4
figure 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 τ s

Moving further up on the uppermost surface there exists a line of Hopf bifurcations, H + , indicating the boundary between S r + and a second attracting surface of equilibria, S a + . Stable limit cycles emanate from the Hopf line through supercritical Hopf bifurcations and terminate in homoclinic cycles, which are homoclinic to the lower fold curve, F (this bifurcation structure is known as a saddle node on invariant cycles, or SNIC, bifurcation). The Hopf line lies outside the physiological range of interest (i.e., w<0); however, the associated limit cycles emanate from outside this region of interest and terminate just within. Thus, the layer problem contains two stable attractors; the lower branch of the critical manifold and the set of stable limit cycles.

4.2 Reduced Problem

The reduced problem describes the dynamics of the slow variables (w,s) on the critical manifold. Time is first scaled by ε to obtain the following slow system:

ε δ v = f 1 ( v , m , h , n , w , s ) , ε δ m = f 2 ( v , m ) , ε h = f 3 ( v , h ) , ε n = f 4 ( v , n ) , w = g 1 ( v , w ) , s = g 2 ( s ) ,
(9)

where x denotes the derivative of x with respect to the slow time-scale, t slow =ετ. In the limit as ε0 of system (9), we derive the reduced problem, which approximates the slow components (w,s) of the full system dynamics

0 = f 1 ( v , m , h , n , w , s ) , 0 = f 2 ( v , m ) , 0 = f 3 ( v , h ) , 0 = f 4 ( v , n ) , w = g 1 ( v , w ) , s = g 2 ( s ) .
(10)

The reduced problem is a two-dimensional differential algebraic equation. Within (10), the first four equations dictate that the dynamics occur on the critical manifold S 0 , and the last two equations describe the dynamics of the slow variables w and s thereon. Note the algebraic equation defining S 0 , f 1 (v, m ˜ (v), h ˜ (v), n ˜ (v),w,s)=0, can be solved for w=w(v,s) respectively s=s(v,w), but not for v reflecting the folded geometry of the manifold; see Fig. 4. Hence, it suffices to study the flow on the critical manifold S 0 in one single coordinate chart, either the (v,s)-chart or the (v,w)-chart where S 0 is defined as a graph.

By definition, the reduced vector field is in the tangent bundle of the critical manifold S 0 . This condition is encoded in the total time derivative of f 1 =0, i.e.,

f 1 ( v , m ˜ ( v ) , h ˜ ( v ) , n ˜ ( v ) , w , s ) =Φ(v,w,s) v + f 1 w w + f 1 s s =0,
(11)

where, once again, a prime denotes a derivative with respect to slow time and

Φ(v,w,s):= f 1 v + f 1 m d m ˜ d v + f 1 h d h ˜ d v + f 1 n d n ˜ d v .
(12)

Projecting the reduced problem onto the (v,s)-chart gives

Φ ( v , w ( v , s ) , s ) v = f 1 w g 1 ( v , w ( v , s ) ) + f 1 s g 2 ( s ) , s = g 2 ( s ) ,
(13)

while projecting the reduced problem onto (v,w)-chart gives the equivalent system

Φ ( v , w , s ( v , w ) ) v = f 1 w g 1 ( v , w ) + f 1 s g 2 ( s ( v , w ) ) , w = g 1 ( v , w ) .
(14)

Both systems are singular along the fold curves, that is, where Φ(v,w,s)=0. The singularity can be removed by rescaling time by Φ(v,w,s). This yields the desingularized reduced system

v = f 1 w g 1 ( v , w ( v , s ) ) + f 1 s g 2 ( s ) , s = Φ ( v , w ( v , s ) , s ) g 2 ( s ) ,
(15)

respectively,

v = f 1 w g 1 ( v , w ) + f 1 s g 2 ( s ( v , w ) ) , w = Φ ( v , w , s ( v , w ) ) g 1 ( v , w ) .
(16)

The flow described by (15), respectively (16), is equivalent to that of (13), respectively (14), on the attracting surface, S a , and the repelling surface, S r + , while reversed on the repelling surface, S r . This is due to the rescaling of time by Φ(v,w,s), which is positive on S a and S r + and negative on S r .

4.2.1 Canards Form a Separatrix for Solutions of the Reduced Problem

Here, we aim to use the geometry of the reduced system in order to identify a manifold that separates the trajectories of the reduced problem into two distinct behaviors: those that return to an equilibrium or rest state and those that proceed to the fold curve where the fast dynamics again become important in establishing a spiking solution. Important to the identification of this manifold is finding folded singularities. We find a folded saddle equilibrium; the canard solution of which is the relevant manifold that separates the behavior of the reduced trajectories.

The desingularized system (15), respectively (16), possesses, in general, two types of equilibria: ordinary singularities and folded singularities,

Q O = { ( v , m , h , n , w , s ) S 0 F | g 1 ( v , w ) = g 2 ( s ) = 0 } , Q F = { ( v , m , h , n , w , s ) F | f 1 w g 1 ( v , w ) + f 1 s g 2 ( s ) = 0 } .
(17)

Ordinary singularities ( g 1 = g 2 =0) are equilibria of the desingularized flow (15), respectively (16), of the reduced flow, (10), and of the original system, (9). Folded singularities, (Φ=0, f 1 w g 1 + f 1 s g 2 =0) on the other hand are generally not equilibria of the reduced flow or the original system.

Within (15), respectively (16), we find three stable nodes, eq 1 , eq 2 and eq 3 , which constitute the set of ordinary singularities, Q O . Note that eq 1 lies on S r + , eq 2 lies on S r and eq 3 lies on S a . The position and stability of these equilibria are independent of the value of the parameter τ s .

We also find folded singularities on both the upper and lower folds of S 0 . In contrast to the set of ordinary singularities, the position and stability properties of the folded singularities are highly dependent on τ s . Only singularities on the lower fold curve are expected to influence the dynamics significantly within the reduced problem as only this fold lies near a physiologically relevant attracting sheet of equilibria. For all values of τ s here considered, there exists a single folded saddle, p fs , on the lower fold curve of S 0 in the physiologically relevant range, s[0,1]. The local geometry of a folded saddle singularity is illustrated within Fig. 5. Recall that we reverse the direction of the desingularized flow on the repelling middle branch of S 0 to obtain the correct reduced flow (Fig. 5b). We identify singular folded saddle canards as trajectories on S 0 , which pass through the folded singularity, p fs , from the attracting critical manifold surface, S a , and follow the repelling critical manifold surface, S r , for a considerable amount of time. Here, we find two canard trajectories; the true canard, ξ t , and the faux canard, ξ f . The faux canard is the particular solution trajectory that lies initially on the repelling surface of the manifold, passes through the folded singularity and onto the attracting surface of the manifold (cf. Figs. 5 and 6). These canards correspond to the global invariant manifolds of the folded singularity, i.e., ξ t = W s ( p fs ) and ξ f = W u ( p fs ). A diagram of the reduced problem, around F , is given for τ s =15 in Fig. 6. The basic structure of the reduced problem, under variation of the parameter τ s , is stretched further along the fold curve in the positive s-direction for larger τ s and contracts toward s=0 for smaller τ s values.

Fig. 5
figure 5

Folded saddle geometry and associated trajectories. The geometry of a generic folded saddle. The folded saddle (purple) is denoted p 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
figure 6

The desingularised reduced problem. The desingularized reduced problem, (15), near the lower fold is projected onto (v,s)-space. Here, τ s =15, although the basic structure shown here is common to all values of τ s analyzed. The lower fold curve, F (gray dashed), is indicated. The folded saddle, p fs (purple), lies on the fold curve and gives rise to two invariant manifolds: a stable invariant manifold, W s ( p fs ) (red), and an unstable invariant manifold, W u ( p 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 eq 2 on S a (orange) and the unstable manifold terminates at 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

The folded saddle canard W s ( p fs ) is a separatrix and effectively organizes the solution trajectories of the reduced problem. Depending on which side of W s ( p fs ) the trajectory is initially, solution trajectories of (15) travel along W s ( p fs ) and either meet the fold curve for s> s fs or move toward, later traveling along, W s ( p fs ). In the former case, the trajectory is no longer accurately approximated by the reduced problem at the fold due to a finite time blow up of system (15) and subsequent dynamics are dictated by the layer problem. In the latter case, the solution trajectory terminates at eq 3 , prevented by W s ( p fs ) from approaching F . Compare Figs. 5 and 6.

5 Firing Threshold Manifolds and Dynamic Inhibition

Using information from the reduced and layer flows, we are able to give a singular limit prediction of post-inhibitory rebound spiking as observed in system (1). Recall that the critical manifold is given by the set of equilibria of the layer problem. Trajectories of the layer problem approach this set along so-called fast fibres. In this context, the critical manifold forms the set of base points of these fast fibers. Hence, base points allow for a connection between the flows of the layer and reduced problems. In particular, the relationship between the base point of the fast fibre through the initial condition due to inhibition and the singular canard solution determines whether a trajectory goes on to spike or not; see Fig. 9.

Singular global trajectories are constructed as continuous concatenations of reduced problem and layer problem solution trajectories. These concatenated dynamics are presented in (v,w,s)-space in Fig. 7. In modeling post-inhibitory rebound spiking, the state point is initially set at eq 3 . This stable node equilibrium simulates the resting membrane potential of the neuron. To model a pulse of inhibitory synaptic input, the trajectory undergoes an instantaneous shift in the s-direction such that s 0 =0.714. Recall that this shift simulates a GABAergic IPSC within a network of neurons. From this point, the dynamics of the system take effect. This synaptic input necessarily translates the trajectory (instantaneously) to a position not on S 0 due to the shape of the manifold. Hence, the layer problem describes the initial motion from this point. Here, the trajectory rapidly approaches an attracting surface of the critical manifold, invariably S a , along a fast fibre through this initial condition. This initial fast motion of the singular solution trajectory toward the manifold is best seen in Fig. 7c.

Fig. 7
figure 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 τ 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

The specific base point of the fast fibre through the initial condition acts then as an initial condition for the reduced system dynamics. This base point is independent of τ s due to the lack of s dynamics within the layer problem. As determined above, there exists a general organizing structure within the reduced problem that determines the subsequent motion of a given trajectory. This structure, termed the canard separatrix, is the portion of the singular canard, ξ t , which lies on S a . The position of the separatrix varies according to the value of τ s while the specific base point on S 0 remains constant. If the base point of the solution trajectory falls within the region on S a bound by ξ t and F , i.e., the shaded region in Figs. 5 and 6, then the trajectory evolves toward the fold curve at which point the dynamics of the trajectory are described by the layer problem. Here, the layer problem dictates that the trajectory undergoes a (fast) oscillation. We observe that the trajectory oscillates only once, during which the reduced problem dictates an average shift in the direction of positive w, and positive s. This shift causes the trajectory to “fall off” the oscillatory regime and return to S a on the other side of the canard separatrix. According to the reduced problem, the trajectory finally comes to rest at eq 3 . This singular global trajectory corresponds to a successful post-inhibitory rebound spike. If, however, the base point falls within the region on S a bound by ξ t and ξ f and characterized by smaller s-values, the solution trajectory is guided along the true and then the faux canard, coming to rest at eq 3 without any fast oscillatory behavior. This singular global trajectory corresponds to a failure to spike after inhibition. The upper transition of spiking behavior is illustrated in Fig. 8.

Fig. 8
figure 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 τ s values for which rebound spiking occurs; τ s [5,24]. The layer problem dictates that the trajectory has a base point on S a independent of τ 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 τ 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 τ s =25 in (v,w)-space. This trajectory lies to the left of the separatrix and evolves in time toward eq 3 (green). This singular prediction corresponds to an unsuccessful rebound spike. An animation of this figure under variation of τ s is given within Additional file 1

The Singular Limit Predicts a Window of Rebound Spiking with Increasing GABAa Inhibition We obtain a singular limit prediction of post-inhibitory rebound spiking by simply identifying the position of the base point of the initial condition relative to the position of the singular canard under variation of τ s . Figure 9 shows that rebound spiking is predicted for τ s [5,24]. This is a good approximation to the range of rebound spiking within the full six-dimensional system; τ s [8,21]. An additional movie file illustrates this shifting of the folded-saddle canard and its effect on spiking behavior in the singular limit [see Additional file 1]. Looking at Fig. 9, the corresponding shift in the s-value falls within the order of the singular perturbation parameter ε, as expected by a singular perturbation analysis.

Fig. 9
figure 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 τ s . For integer values of τ 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 τ s for which there exists rebound spiking. An animation of this figure under variation of τ s is given within Additional file 1

5.1 Nonsingular Canards: The True Threshold Manifolds for Spiking Activity

An important result from GSPT is that a singular canard of folded saddle type perturbs to a nearby canard for the full model problem [21]. Consequently, the perturbed canard forms the true firing threshold manifold of the full system (1). In the following, we confirm numerically this firing threshold manifold by calculating the canard. We note here that due to numerical challenges in calculating canards for systems in R n , n>3 with more than one fast variable, we consider a reduction of the six-dimensional modified propofol model. Resolving this issue is left as a point of focus for future work. By setting each of m, h and n to their respective steady-state values, we derive a three-dimensional reduction of our original system. This reduction is valid as we have determined that m, h and n evolve on a fast time-scale; this system simply approximating this fact by having these processes act instantaneously. Within this reduced system the range of τ s values for which spiking occurs is similar to the original with τ s [9,21]; however, the spiking mechanism is partially disabled and the trajectory is unable to fully reset. At least 2 fast variables are necessary to provide the necessary repolarization mechanism. We note that this does not affect sub or perithreshold dynamics, and consequently does not affect the onset of spiking. Hence, a model reduction to 1 fast variable locally near the fold is justified.

Using the continuation package AUTO [7], and closely following the techniques outlined in [5], the nonsingular canard, Ξ t , is found for τ s =15 and is then continued in τ s . The nonsingular canards of the modified propofol model lie at the intersection of the attracting and repelling slow manifolds, S ε a and S ε r , respectively. These manifolds correspond to S a and S r , respectively, for ε>0. The attracting and repelling manifolds of the perturbed system are calculated using the homotopy continuation of solution trajectories to a suitable boundary value problem. Here, we make use of the normal hyperbolicity of the critical manifold S 0 . Namely, for small ε>0, the slow manifolds S ε a and S ε r are smooth perturbations of the critical manifold S 0 , away from the fold curve F where normal hyperbolicity is lost. A detailed description of slow invariant manifold calculations and canard detection and continuation is provided for the self-coupled FitzHugh–Nagumo system within the AUTO manual, demo fnc.

Within Fig. 10, we find that the nonsingular canard accurately forms the boundary of the spiking and nonspiking behavior within the three-dimensional modified propofol model. We thus confirm that the canard does in fact form the threshold manifold of spiking activity. This nonsingular analysis confirms our geometric understanding of the singular limit prediction on both interval boundaries.

Fig. 10
figure 10

Nonsingular canards and global trajectories. The true nonsingular canard (blue) for the three-dimensional modified propofol model is located for τ s =15. This canard is then continued in τ s using a suitable boundary condition problem. We thus obtain the true canard for any value of τ 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 eq 3 (green) are indicated. a The nonsingular trajectory and associated canard for τ s =8 in (v,w)-space. b The nonsingular trajectory and associated canard for τ s =9 in (v,w)-space

6 Firing Threshold Manifolds and the Classical Step Protocol

We have thus far shown that it is necessary to consider the effects of dynamic inhibition in order to understand the cessation of spiking for large τ s . Here, we analyze the current step protocol and seek to understand why only a single transition of activity can be found. We again use singular perturbation analysis to identify the spiking threshold manifold. In the singular limit, this manifold is given by the concatenation of the middle branch of the one-dimensional critical manifold and the fast fibre through a fold point when I app =0 (see Fig. 11).

Fig. 11
figure 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 syn =3.5× 10 4 to simulate inhibition. We compare the critical manifold for each value of I syn . On each manifold, the lower branch (magenta) is linearly stable, the upper branch (gray) linearly unstable. The stable node equilibria (green), eq 3 and eq ¯ 3 , and saddle equilibrium (orange), eq 3 are indicated. Initially, the system starts at rest on the node equilibrium of the I 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 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

6.1 Singular Perturbation Analysis

In a singular perturbation analysis similar to that above, we consider four fast variables, (v,m,h,n), and a single slow variable, w. Note that s dynamics are no longer considered in a current step protocol; I syn , and thus dependence on s, is here replaced by a constant value. In the singular limit, the layer problem defines a one-dimensional cubic-shaped critical manifold in (v,w)-space. As per the above analysis, only the lower branch of the critical manifold is an attracting branch; the upper two branches being unstable. While the shape and stability properties of the critical manifold are independent of the value of I syn , the position of the critical manifold is shifted in the direction of negative w as the synaptic current is set to I syn =3.5× 10 4 during inhibition (otherwise I syn =0; see Fig. 11). Note that here the critical manifold for I syn =0 is precisely the section defined by {s=0} through the critical manifold derived above, i.e., the critical manifold for I syn =0 is here given by S 0 {s=0}. Once on the critical manifold, dynamics are described by the reduced problem; here formulated as (10) with s =0 and I syn adjusted as per the step protocol. The reduced problem reveals, as before, the existence of three stable node equilibria, one on each branch of the critical manifold. We label the equilibria as before for I syn =0, whereas the equilibria on the shifted critical manifold, i.e., for I syn =3.5× 10 4 , are denoted eq ¯ i , i=1,2,3.

Initially, the state point is held at rest at eq 3 with I syn =0. At the onset of synaptic inhibition the critical manifold shifts, and so defines a new stable equilibrium, eq ¯ 3 . In the singular limit (Fig. 11), the trajectory falls instantaneously, traveling along a vertical fast fibre, onto the stable branch of the shifted critical manifold. This base point is thus a projection of eq 3 onto the shifted critical manifold. The singular trajectory is then described by the reduced problem, slowly moving along this lower branch toward eq ¯ 3 . Once synaptic inhibition is removed, the critical manifold shifts back to its original position. At this point, the layer problem dictates a rapid vertical ascent; the trajectory no longer remains on the critical manifold. The subsequent dynamics are thus dependent on the length of the applied constant synaptic inhibition. If the duration of inhibition is relatively short, at the removal of inhibition, the trajectory simply returns to the now overhead lower branch of the critical manifold. This corresponds to an unsuccessful post-inhibitory rebound spike. If; however, the duration of inhibition allows the trajectory time enough to sufficiently approach eq ¯ 3 then, at the removal of inhibition, the trajectory undergoes a spike in v. This behavior corresponds to a successful post-inhibitory rebound spike. The associated threshold manifold is given as the concatenation of the middle branch of the critical manifold and the fast fiber through the lower fold point, shown in Fig. 11.

6.1.1 The Singular Limit Predicts Only One Spiking Transition with Increasing Duration of Inhibition

We here determine that the singular limit predicts the minimum duration of inhibition such that a trajectory is still able to spike is 6 ms. Since the layer problem acts instantaneously, this value is determined solely according to the dynamics within the reduced problem. We simulate the singular reduced problem along the lower branch of the critical manifold, noting the time at which the singular trajectory passes the w-value of the lower fold. Note here we use the singular problem to avoid the distortion of time within the desingularized problem due to a position dependent rescaling of time in the process of desingularization.

6.2 Nonsingular Firing Threshold Manifold

Here, as expected, the singular limit accurately mimics the behavior of the nonsingular system (compare Figs. 11 and 12). As per the singular limit prediction, the duration of inhibition determines whether or not a given trajectory spikes in the nonsingular system (Fig. 12). We observe that for short step protocols, the corresponding trajectory makes only a small excursion down toward the shifted critical manifold without reaching it. Once synaptic inhibition is released the trajectory returns to eq 3 shortly after. In a voltage time-trace, this behavior manifests in a short dip in v before a return to the resting membrane potential (Fig. 3a). As the step protocol is lengthened, the trajectory reaches the shifted critical manifold, afterward slowly moving toward eq ¯ 3 (Fig. 12a). For synaptic inhibition lasting longer than 14 ms, when inhibition is removed, the local attractor is found within a fast oscillatory regime. Here, we observe a small dip in v followed by a spike within the voltage time-trace (Fig. 3b).

Fig. 12
figure 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 syn =0 and I syn =3.5× 10 4 . On each manifold, the lower branch (magenta) is linearly stable, the upper branch (gray) linearly unstable. The corresponding stable node equilibria (green), eq 3 and eq ¯ 3 , and saddle equilibrium (orange), eq 2 , are indicated. Initially, the system starts at rest at the node equilibrium of the I 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 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 x ˜ = x ˜ (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

6.2.1 A Canard Solution Forms the Nonsingular Separatrix Between Spiking and Nonspiking Solutions

As per the analysis of dynamic inhibition above, we again find a geometric object which demarcates the boundaries of trajectories with different spiking behavior. Within the nonsingular system the perturbed separatrix is located using a shooting method, in backward time, from the stable node eq 2 . As before, the propofol model used to calculate this separatrix requires the reduction x ˜ = x ˜ (v), for x=m,h,n, in order to allow calculation along an unstable manifold. Here, we see that the nonsingular separatrix forms the boundary of spiking and nonspiking behavior (compare Figs. 12a and 12b). We note here that the threshold manifold is a canard, as per the dynamic inhibition analysis. Accordingly, we find that trajectories, which begin exponentially close to this structure follow it for a significant amount of time, even onto a repelling portion of the associated critical manifold, such as that in Fig. 12b near the fold point. The deviation of the perturbed threshold manifold compared with that of the singular analysis near the critical manifold lower fold is explained by a fold analysis within singular perturbation theory (compare Figs. 11 and 12). Here, we expect, and indeed find, that the nonsingular manifold perturbs a distance which is O( ε 2 / 3 ) from the fold [16, 22]. Hence, we have identified the non-singular threshold manifold for the current step protocol.

6.2.2 Explanation of a Finite Number of Rebound Spikes as Inhibition Time Is Increased

As the duration of the step protocol increases, the number of rebound spikes increases steadily to a maximum of three spikes (compare Figs. 3d and 12b). This observation is explained by a limited reset in the direction of w during each rebound spike. After a single spike, as seen in Fig. 12b, the corresponding trajectory remains to the right of the threshold manifold. Three consecutive rebound spikes are required before the trajectory returns to an excitable resting state at eq 3 , having moved past the threshold manifold. At this point, we see that further prolonging the duration of synaptic inhibition has no additional effect. Once the trajectory is sufficiently close to eq ¯ 3 during synaptic inhibition, regardless of the actual duration of inhibition, the subsequent dynamics remain unchanged. This results in a maximum number of rebound spikes. This maximum spike number is encoded within the time-scale separation between the fast and slow dynamics, i.e., within ε. For fixed ε, the maximum spike number is set according to the average shift in the slow dynamics per (fast) spike event. As ε0, we observe that the maximum spike number increases as the reset in w reduces per spike (work not shown). The singular limit picture confirms this finding, here showing no net shift in w per spike for ε=0.

6.3 The Folded Saddle Canards of the Dynamic Protocol Converge Toward the Spiking Threshold Manifold of the Current Step Protocol as τ s 0

In the modified propofol model, we are able to explain both spiking transitions through identification of the firing threshold manifold in a 4-fast/2-slow time-scale separation setting—the threshold manifold being the folded saddle canard. Using the software package AUTO, we continue this folded saddle canard toward small τ s values. Figure 13 shows the firing threshold of the modified propofol model for 1 τ s 15 (solid purple) as well as the firing threshold manifold of the current step protocol (dashed orange). Clearly, the folded saddle canards converge toward the firing threshold manifold of the current step protocol as τ s 0. We note that as τ s 0 we consider s a fast variable within a GSPT analysis. Given that the s dynamics are decoupled from the rest of the system, we can easily identify the stable equilibrium state; s=0. Thus, in the limit as τ s 0, the system dynamics and, therefore, the firing threshold manifold more closely represents that of a system where the s-dynamics rapidly equilibrate to s=0; i.e., the step protocol where I syn =0. Hence, we show the intimate relationship between the two identified distinct firing threshold manifolds in the modified propofol model, the folded saddle canard, and the stable manifold of the saddle equilibrium.

Fig. 13
figure 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 τ s between 1–15 (solid purple). The firing threshold manifold for static inhibition is the stable manifold of eq 2 (orange dashed) shown earlier in Fig. 12. For decreasing values of τ 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 x ˜ = x ˜ (v), for x=m,h,n, to allow calculation along repelling manifolds

7 Comparison of the Original and Modified Propofol Models

The propofol network model presented here is a modified form of the propofol neuron network model in [18]. The propofol model is the same 6–dimensional system of equations (Eqs. (1) and (2)) as the modified propofol model but with two parameter changes. The modified propofol model is phenomenologically similar to the original propofol model, which shows rebound spiking only for values of τ s between 8 and 48 ms (compared with the modified model: 8–21 ms) and was also examined using geometric singular perturbation theory [18]. In particular, a canard of folded saddle type was identified as a firing threshold manifold. In Sects. 7.1–7.3, we detail our reasons for modifying the original propofol model.

7.1 Modifying the Propofol Model: Increasing the Time-Scale Separation Allows for a more Accurate Singular Limit Trajectory Approximation

Both models, the original and the modified propofol model, have a global stable equilibrium that plays the role of the resting membrane potential. We note that this stable equilibrium (v-value of −63.6 mV) of the original propofol model lies quite close to the fold curve, F . Additionally, we observe that the fold-curve F has an almost constant v-value (v63 mV along F ). Given the fairly uniform structure of the critical manifold S 0 in the direction of s near F , the position of the stable equilibrium ensures a proximity of the post-inhibitory initial condition and thus of the post-inhibitory dynamics to F .

This proximity to a fold structure is known to affect the time-scale separation of variables within a singularly perturbed problem. More specifically, near normally hyperbolic equilibria, fast and slow dynamics within a singularly perturbed system evolve as O(1) and O(ε), respectively. Within a neighborhood of a fold structure, at which normal hyperbolicity is lost, the original fast/slow splitting does not persist and the fast and slow dynamics evolve predominantly on an intermediate time-scale O( ε 1 / 3 ) and a slow time-scale O(ε) [22]. This explains why we find an appreciable deviation between the singular and nonsingular trajectories in the original propofol model (Fig. 14a).

Fig. 14
figure 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 τ s =20 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

Thus, we propose the modified propofol model of post-inhibitory rebound spiking. This modified propofol model retains the basic phenomenological and geometric features of the propofol model while the position of the stable equilibrium has been moved away from the fold curves. This modification allows for a more marked separation of time-scales during the early dynamics, and thus describes a system more accurately approximated by a singular limit prediction (Fig. 14b).

To formulate the modified propofol model the steady-state function w (v) is shifted by 3 mV in the direction of negative v. This results in a lower global stable equilibrium point (v-value of −65.8 mV), and thus a lower resting membrane potential. However, this modification also makes it more difficult to generate a rebound spike. In order to counter this effect, the value of the maximal synaptic conductance g i is increased so as to effectively increase the strength of synaptic inhibition. Here, we set g i to 4 mS, i.e., g ¯ i =0.04.

As previously noted, the behavior of the modified propofol model is similar to the original propofol model in that post-inhibitory rebound spiking occurs only for intermediate values of τ s (Fig. 2). Thus, by our modifications, we have not lost the rebound-spiking “window” and we have gained a more accurate geometrical representation of nonsingular trajectories by their singular limit. However, even more importantly, the modified model allows us to use the same fast/slow decomposition to predict the spiking transition at both interval boundaries ( τ s small and τ s large). This was not possible in the original propofol model. We describe this in more detail next.

7.2 Discrepancies Between the Singular and Nonsingular Trajectories in the Propofol Model

Although the singular limit of the original propofol model provides a good approximation for cessation of spiking when τ s is large (52 ms for the singular limit versus 48 ms for nonsingular trajectories), when we compare singular global trajectories with the nonsingular trajectories of the propofol model, we find that the geometry of the singular limit trajectories does not accurately predict the geometry of the nonsingular trajectories (Fig. 14). In particular, the early dynamics of the singular limit trajectories do not accurately mimic those of the nonsingular trajectories using the proposed slow/fast splitting. As mentioned before, this poor predictive power of geometric singular perturbation theory can be explained by the proximity of the trajectory near the fold F , throughout its evolution. The time-scale splitting of the “fast” (v,m,h,n) and “slow” (w,s) variables away from the fold F does not hold anymore in a neighborhood of the fold F . Without this time-scale splitting, and thus the identification of a singularly perturbed problem, regular GSPT analysis does not yield a reliable prediction of system dynamics.

We thus conclude that while the singular limit analysis still provides a good prediction for the large τ s spiking transition in the original propofol model, the particular geometry of this system distorts the time-scale separation, and thus requires an alternative approach of geometric singular perturbation theory using a blow-up analysis along the fold-curve [22] followed by a blow-up of the folded saddle singularity [21, 24]. This two-step approach is left for future work.

7.3 In the Original Propofol Model, Considering s Slow or Fast Does not Explain the Spiking Transition when τ s Is Small

In the original propofol model, considering s as a slow variable cannot predict the spiking transition when τ s is small; see Fig. 10 in [18]. Similarly, considering s as a fast variable cannot predict the spiking transition when τ s is small either; see explanation on p. 13 in [18]. Hence, using geometric singular perturbation theory in the context of either a 4-fast/2-slow time-scales or a 5-fast/1-slow time-scales separation does not predict the transition between spiking and nonspiking when τ s is small in the original propofol model. The answer to this “riddle” about the appropriate dynamics of s for small τ s in the original propofol model lies in an intermediate time-scale of order O( ε α ), 0<α<1. Such an intermediate time-scale reveals itself when using a (cylindrical) blow-up analysis of the fold, F [22, 24]. As mentioned above, we leave this blow-up analysis of the original propofol model for future work.

7.4 Comparison of Geometric Structures in the Original and Modified Propofol Model when τ s Is Small

In [18], the authors use a 5-fast/1-slow time-scale separation to examine the geometry of the spiking transition when τ s is small. Interestingly, this results in the same resting state ( I syn =0) singular limit geometry as the 4-fast/1-slow time-scale problem used here to analyze the current step protocol (Fig. 12). This makes sense because considering s fast allows the manifold {s=0} to be reached on a fast time-scale, and thus the system reduces to the same 4-fast/1-slow time-scale problem used here. As shown in [18], when s is considered one of the fast variables, as in the original propofol model, the stable manifold of eq 2 in Fig. 12 approximates well the threshold for spiking. Recall that this is shown in Fig. 13, i.e., as τ s 0 the canard separatrix is well approximated by the firing threshold manifold of the current step protocol. Hence, by changing parameters to allow the construction via geometric singular perturbation theory of canards, this paper clarifies the underlying geometry of the original propofol model.

8 Concluding Remarks

An important feature of most physiological systems is that they evolve on multiple time-scales. The theory of differentiable dynamical systems for two time-scales (slow/fast) has a successful history in explaining a wide range of physiological behavior such as electrical spiking and bursting in neurons [13]. Complex pattern generation in such slow/fast systems is almost exclusively related to loss of normal hyperbolicity of invariant critical manifolds, which is associated with bifurcation sets in the fast subsystem. In particular, folded critical manifolds are ubiquitous in such systems (see, e.g., Fig. 4). Canard theory deals exactly with these slow/fast time-scales systems where loss of normal hyperbolicity occurs, and its theory is applicable to problems with arbitrary dimensions [2, 3, 8, 15, 16, 21, 23, 24]. A recent success story of canard theory is that it provides an explanation for mixed-mode oscillations (MMOs), a frequently observed mix of small and large amplitude oscillation patterns in slow/fast time-scale physiological models. Canards of folded node and folded saddle-node type play a key role in explaining these patterns. The interested reader is referred to the current review on MMOs [5], and the extensive reference list to applications therein.

Canard theory provides also a new direction for understanding transient dynamics of biological systems that have multiple time-scales. The propofol model studied here is a prime example. We demonstrate here the use of canard theory to explain the dynamics of rebound spiking for a specific range of propofol doses. In the context of neuronal excitability, we identify canards of folded saddle type as firing threshold manifolds. It is remarkable that dynamic information such as the temporal evolution of an external drive (GABAergic inhibition in this study) is encoded in the location of an invariant manifold—the canard. It is the variable positioning of the canard separatrix that explains the observed rebound spiking for a specific range of propofol doses (Fig. 9). The same role of folded saddle canards as firing threshold manifolds was recently identified in a class of Morris–Lecar/FitzHugh–Nagumo type models with dynamic external drive [25]. Since mathematical models of physiological phenomena (both neuronal and nonneuronal) frequently show abrupt transitions in behavior and have dynamics which are encoded by multiple time-scales, the methods of GPST and canard theory are likely applicable to a much broader range of problems within the biological sciences.

Appendix:  Steady State and Time-Scale Functions

The formulations for the steady state functions ( m (V), h (V), n (V)) and time-scale functions ( τ m (V), τ h (V), τ n (V)) are taken from [18]. The steady state function and time-scale function for w are here modified. These functions take the respective forms

x ( V ) = α x ( V ) α x ( V ) + β x ( V ) , τ x ( V ) = α x ( V ) α x ( V ) + β x ( V )
(18)

for x=m,h,n,w. The constituent functions α x (V) and β x (V) are given by

α m ( V ) = 0.32 ( V + 54 ) 1 exp [ ( V + 54 ) / 4 ] , β m ( V ) = 0.28 ( V + 27 ) 1 exp [ ( V + 27 ) / 5 ] , α h ( V ) = 0.128 exp [ ( V + 50 ) / 18 ] , β h ( V ) = 4 1 + exp [ ( V + 27 ) / 5 ] , α n ( V ) = 0.032 ( V + 52 ) 1 exp [ ( V + 52 ) / 4 ] , β n ( V ) = 0.5 exp [ ( V + 57 ) / 40 ] , α w ( V ) = 3.209 10 4 ( V + 33 ) 1 exp [ ( V + 33 ) / 9 ] , β w ( V ) = 3.209 10 4 ( V + 33 ) 1 exp [ ( V + 33 ) / 9 ] .
(19)

The gating variable steady-state and time-scale functions are plotted in Fig. 15. The time derivative of each variable over the course of a single spike is given in Fig. 16. This figure indicates the actual time-scale separation and is used, in part, to justify a time-scale splitting as detailed in Sect. 3.

Fig. 15
figure 15

Gating variable steady-state and time-scale functions. a The four gating variable steady-state functions, m (V) (red), h (V) (orange), n (V) (green), and w (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, τ m (V) (red), τ h (V) (orange), τ n (V) (green), and τ 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 τ m (V), τ h (V) and τ n (V) evolve on much smaller time-scales than τ w (V)

Fig. 16
figure 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 d w d t and d s d t can be seen

Electronic Supplementary Material

References

  1. Bai D, Pennefather PS, MacDonald JF, Orser BA: The general anesthetic propofol slows deactivation and desensitization of GABA A receptors. J Neurosci 1999, 19(24):10635–10646.

    Google Scholar 

  2. Benoît E:Systèmes lents-rapides dans R 3 et leur canards. Astérisque 1983, 109–110: 159–191.

    Google Scholar 

  3. Benoît E, Callot J, Diener F, Diener M: Chasse au canard. Collect Math 1981, 31–32: 37–119.

    Google Scholar 

  4. Brøns M, Kaper T, Rotstein H: Introduction to focus issue: mixed mode oscillations: experiment, computation, and analysis. Chaos 2008., 18: Article ID 015101 Article ID 015101

    Google Scholar 

  5. Desroches M, Guckenheimer J, Krauskopf B, Kuehn C, Osinga H, Wechselberger M: Mixed-mode oscillations with multiple time-scales. SIAM Rev 2012, 54: 211–288. 10.1137/100791233

    Article  MathSciNet  Google Scholar 

  6. Desroches M, Krupa M, Rodrigues S: Inflection, canards and excitability threshold in neuronal models. J Math Biol 2013, 67(4):989–1017. 10.1007/s00285-012-0576-z

    Article  MathSciNet  Google Scholar 

  7. Doedel EJ: AUTO: a program for the automatic bifurcation analysis of autonomous systems. Congr Numer 1981, 30: 265–284.

    MathSciNet  Google Scholar 

  8. Dumortier F, Roussarie R: Canard cycles and center manifolds. Mem Am Math Soc 1996., 121: Article ID 577 Article ID 577

    Google Scholar 

  9. Fenichel N: Geometric singular perturbation theory. J Differ Equ 1979, 31: 53–98. 10.1016/0022-0396(79)90152-9

    Article  MathSciNet  Google Scholar 

  10. FitzHugh R: Mathematical models of threshold phenomena in the nerve membrane. Bull Math Biophys 1955, 7: 252–278.

    Google Scholar 

  11. FitzHugh R: Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1961, 1(6):445–466. 10.1016/S0006-3495(61)86902-6

    Article  Google Scholar 

  12. Hodgkin AL, Huxley AF: A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 1952, 117: 500–544.

    Article  Google Scholar 

  13. Izhikevich EM: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press, Cambridge; 2007.

    Google Scholar 

  14. Jones CKRT: Geometric singular perturbation theory. Springer Lecture Notes Math. 1609. Dynamical Systems 1995, 44–120.

    Chapter  Google Scholar 

  15. Krupa M, Szmolyan P: Relaxation oscillations and canard explosion. J Differ Equ 2001, 174: 312–368. 10.1006/jdeq.2000.3929

    Article  MathSciNet  Google Scholar 

  16. Krupa M, Szmolyan P: Extending geometric singular perturbation theory to nonhyperbolic points—fold and canard points in two dimensions. SIAM J Math Anal 2001, 33(2):286–314. 10.1137/S0036141099360919

    Article  MathSciNet  Google Scholar 

  17. McCarthy MM, Brown EN, Kopell N: Potential network mechanisms mediating electroencephalographic beta rhythm changes during propofol-induced paradoxical excitation. J Neurosci 2008, 28(50):13488–13504. 10.1523/JNEUROSCI.3536-08.2008

    Article  Google Scholar 

  18. McCarthy MM, Kopell N: The effect of propofol anaesthesia on rebound spiking. SIAM J Appl Dyn Syst 2012, 11(4):1674–1697. 10.1137/100817450

    Article  Google Scholar 

  19. Rush ME, Rinzel J: The potassium A-current, low firing rates and rebound excitation in Hodgkin–Huxley models. Bull Math Biol 1995, 57: 899–929.

    Article  Google Scholar 

  20. Schwartz RS, Brown EN, Lydic R, Schiff ND: General anesthesia, sleep, and coma. N Engl J Med 2010, 363(27):2638–2650. 10.1056/NEJMra0808281

    Article  Google Scholar 

  21. Szmolyan P, Wechselberger M:Canards in R 3 . J Differ Equ 2001, 177: 419–453. 10.1006/jdeq.2001.4001

    Article  MathSciNet  Google Scholar 

  22. Szmolyan P, Wechselberger M: Relaxation oscillations in R 3 . J Differ Equ 2004, 200: 69–104. 10.1016/j.jde.2003.09.010

    Article  MathSciNet  Google Scholar 

  23. Wechselberger M: Existence and bifurcation of canards in R 3 in the case of a folded node. SIAM J Appl Dyn Syst 2005, 4: 101–139. 10.1137/030601995

    Article  MathSciNet  Google Scholar 

  24. Wechselberger M: À propos de canards (Apropos canards). Trans Am Math Soc 2012, 364: 3289–3309. 10.1090/S0002-9947-2012-05575-9

    Article  MathSciNet  Google Scholar 

  25. Wechselberger M, Mitry J, Rinzel J: Canard theory and excitability. Random and Nonautonomous Dynamical Systems in the Life Sciences 2013. in press in press

    Google Scholar 

  26. Wieczorek S, Ashwin P, Luke CM, Cox PM: Excitability in ramped systems: the compost-bomb instability. Philos Trans R Soc A, Math Phys Eng Sci 2011, 467(2129):1243–1269. 10.1098/rspa.2010.0485

    Article  MathSciNet  Google Scholar 

Download references

Acknowledgements

This work was supported by NSF Grant DMS-0602204 (MM), NSF Grant DMS-071670 (NK and MM) and ARC Grant FT-120100309 (MW).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to John Mitry.

Additional information

Competing Interests

The authors declare that they have no competing interests.

Authors’ Contributions

JM, MM, NC, and MW conceived of the analysis and of the overall goals of the project. JM carried out the numerical analysis and wrote the manuscript with assistance from MM, NC, and MW. All authors read and approved the final manuscript.

Electronic supplementary material

Additional file 1: Here we present Fig. 8 and 9 in animated form. As tauS varies, the folded-saddle canard is shifted. This canard forms the firing threshold manifold and so the excitability of the cell is thus altered. The singular limit predicts a range of τ s values for which rebound spiking occurs; τ s [5,24]. a Singular limit (i.e., ε=0) solution trajectories. The singular canard (red) and singular faux canard (blue) organize trajectories and govern excitability. The singular canard forms the firing threshold manifold. b Singular limit prediction of spiking. The base point remains independent of τ s while the position of the singular canard varies with τ s . This alters excitability and gives rise to a window of τ s values for which spiking occurs (MOV 264 kB) (MOV 265 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Mitry, J., McCarthy, M., Kopell, N. et al. Excitable Neurons, Firing Threshold Manifolds and Canards. J. Math. Neurosc. 3, 12 (2013). https://doi.org/10.1186/2190-8567-3-12

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-8567-3-12

Keywords