Coherent neural oscillations induced by weak synaptic noise

AbstractWe analyze the effect of synaptic noise on the dynamics of the FitzHugh–Nagumo (FHN) neuron model. In our deterministic parameter regime, a limit cycle solution cannot emerge through a singular Hopf bifurcation, but such a limit cycle can nevertheless arise as a stochastic effect, as a consequence of weak synaptic noise in a regime of strong timescale separation ($$\varepsilon \rightarrow 0$$ε→0) between the slow and fast variables of the model. We investigate the mechanism behind this phenomenon, known as self-induced stochastic resonance (SISR) (Muratov et al. in Physica D 210:227–240, 2005), by using multiple-time perturbation techniques and by analyzing the escape mechanism of the random trajectories from the stable manifolds of the model equation. Even though SISR occurs only in the limit as the singular parameter $$\varepsilon \rightarrow 0$$ε→0, decreasing $$\varepsilon $$εdoes not increase the coherence of the oscillations in the FHN model, but rather increases the interval of the noise amplitude $$\sigma $$σ for which coherence occurs. This is in contrast to the dynamical system studied in Muratov et al. 
(2005). Moreover, the phenomenon is robust under parameter tuning. Numerical simulations exhibit the results predicted by the theoretical analysis. In fact, our analysis together with that in Yamakou and Jost (Weak noise-induced transitions with inhibition and modulation of neural oscillations, 2017) reveals that the FHN model can support different stochastic resonance phenomena. While it had been known (Pikovsky and Kurths in Phys Rev Lett 78:775–778, 1997) that coherence resonance can occur when the slow variable is subjected to noise, we show that, when noise is added to the fast variable, two other types, inverse stochastic resonance and SISR, may emerge in the same weak noise limit and that the transition between these phenomena can be induced by varying a simple parameter.


Introduction
Since electrical information in the nervous system is encoded, processed, and transmitted by trains of neuronal action potentials [4,5], a major goal in neuroscience is to understand how neurons generate action potentials both spontaneously and in response to possibly random synaptic and ion channel inputs. It is a standard electrophysiological observation that neurons with nearly identical physiological features and external stimuli may react differently and therefore exhibit different neurocomputational properties. Conversely, neurons having different physiological features and external stimuli may show nearly identical activities. It is generally believed that these variances in neuronal behavior are caused by differences in the neurons' bifurcation dynamics. The latter, however, can only be understood with the help of sophisticated tools from nonlinear dynamics.
Noise is ubiquitous in neural systems. It may arise from many different sources, and it may have different effects. In fact, depending on the deterministic parameter regime of a neuron model, the addition of noise to either the membrane potential variable and/or the recovery current variable can induce different dynamical effects. Some of these effects can be somewhat counterintuitive, as noise, instead of being a nuisance, in some settings actually plays a constructive role in signal detection.
The sources of neuronal noise include synaptic noise, that is, the quasi-random release of neurotransmitters by synapses or random synaptic input from other neurons, and channel noise, that is, the random switching of ion channels. Since contributions of the inherently stochastic nature of voltage-gated ion channels to neuronal noise levels are widely assumed to be minimal, because typically a large number of ion channels are involved and fluctuations should therefore average out, the channel noise is frequently ignored in the mathematical modeling. Synaptic noise on the other hand is believed to be the dominant source of membrane potential fluctuations in neurons, and it can have a strong influence on their integrative properties [6]. In this paper, we shall investigate the effects of synaptic noise on the fast membrane potential variable of the FHN neuron model. Here, we do not consider channel noise, which, however, could be included by simply adding it to the slow current recovery variable equation of this neuron model.
Of course, we are not the first to study noise-induced phenomena in the context of computational neuroscience. Let us therefore briefly discuss four such phenomena treated in the literature. The oldest such phenomenon is stochastic resonance (SR). It was first proposed in 1981 by Benzi et al. [7] in order to explain the periodic recurrence of ice ages. The idea behind SR is that a system in a noisy environment becomes more sensitive to a subthreshold periodic perturbation at some optimal finite level of noise intensity. In 1993, Douglass et al. [8] reported the first experimental demonstration of SR in crayfish sensory neurons. This, coupled with the development of a more general characterization of SR as a threshold phenomenon independent of dynamics, led to a widespread interest in applications of the idea in biological systems, in general, and in computational neuroscience in particular [9,10].
In 1993, Longtin [11] demonstrated the occurrence of SR in a computational neuron model. Here, the 2dimensional FHN neuron model is considered and perturbed by a subthreshold external periodic stimulus and a noisy stimulus governed by an Ornstein-Uhlenbeck process. As predicted by the SR scheme (see also [12]), it was shown that the oscillation of the FHN neuron becomes more closely correlated with the subthreshold periodic input current at some optimal level of noise intensity. The stochastic analysis of the FHN neuron and other neuron models has been extensively studied in recent years with various types of noises, multistability, and even how SR can be controlled; see, for example, [13][14][15][16][17][18][19][20][21][22][23].
In 1997, Pikovsky and Kurths [3] reported on another noise-induced phenomenon in an excitable system. They considered the FHN neuron model with the noise term added to the slow recovery variable equation containing the Hopf bifurcation parameter. For this neuron model, they then showed that, even in the absence of an external periodic stimulus (unlike SR), the noise could activate the neuron by producing a sequence of pulses that could achieve a maximal degree of coherence for an optimal noise amplitude, provided the bifurcation parameter is in the neighborhood of its Hopf bifurcation value. They called this phenomenon coherence resonance (CR). Following [3], many research papers have investigated CR in a variety of systems ranging from multiple timescale systems [24][25][26] to systems with no multiple timescales [27,28]. CR has also been extensively studied [29][30][31][32][33][34][35][36] and the relation between CR and chimeras states in network of excitable systems has also been recently reported [37].
In 2005, Muratov et al. discovered and developed the fundamental theory (based on large deviation theory [38,39]) of a new form of coherence behavior. They called this new phenomenon self-induced stochastic resonance (SISR); here, a weak noise amplitude could induce coherent oscillations (a limit cycle behavior) that the deterministic model equation cannot exhibit [1]. SISR has been investigated theoretically and numerical in different systems including Brown-ian ratchets and cancer model [40][41][42][43]. A natural question is: in what ways is SISR different from SR and CR? First, unlike SR and CR, which do not necessarily require a system to evolve on a multiple timescale [27,44], SISR requires the system to possess fast and slow timescales. Moreover, in the complete absence of a deterministic perturbation or in the presence of a subthreshold deterministic perturbation, SISR may occur. This contrasts with SR, which always requires a subthreshold periodic perturbation (in fact, SR is precisely the enhancement of this weak low-frequency harmonic signal by some optimal noise strength). Also, in contrast to CR, SISR does not require the model parameter to be in the immediate neighborhood of its bifurcation value. Thus SISR, unlike CR, is robust to parameter tuning.
In the FHN neuron model, coherent spike trains emerge through very different mechanisms depending on whether the noise term is added to the fast variable (SISR) or to the slow variable (CR); the former is analyzed in the present paper, and for the latter, see [3,24,25,45]. There are, however, examples of multiple timescale systems in which CR occurs when the noise term is added to the fast variable equation, which contains the bifurcation parameter (see, e.g., [26]). What is crucial for CR is the proximity of the parameter to the bifurcations (Hopf or saddle-node bifurcations on a limit cycle).
In all of these phenomena (SR, CR, and SISR), noise has a facilitating effect on the oscillations and increases coherent responses. Noise can, however, also have the opposite effect and turn-off repetitive neuronal activity, as was discovered both experimentally [46] and theoretically [47]. In fact, Gutkin et al. [47][48][49] used neural model equations with a mean input current consisting of both a constant deterministic and a random input component, to computationally confirm the inhibitory and modulation effects of noise on the neuron's spiking activity. They found that there is a tuning effect of noise that has a character opposite to SR, CR, and SISR. They called this other noise-induced phenomenon inverse stochastic resonance (ISR). During ISR, when the intensity of the (weak) noise amplitude is increased, it will first strongly inhibit the spiking activity and decrease the mean number of spikes to a minimal value, but a further increase in the noise amplitude will increase the mean number of spikes again, even above what is observed in the noise-free case. Mathematically, ISR has not been extensively studied as compared to other resonance phenomena; see [50][51][52][53], and it was recently confirmed experimentally in [54].
In the current paper, we analyze the mechanism behind SISR in the FHN neuron model under the influence of synaptic noise. We clarify the different conditions leading to either SISR or ISR (analyzed in our previous paper [2]) in the weak synaptic noise limit. For the FHN model, we obtain the interval in which the noise amplitude σ has to lie in order to induce a coherent spike train (SISR) in the limit as ε → 0. We also determine another very small interval for σ in which the FHN neuron generates only a Poisson sequence of spikes in the limit as ε → 0. One of the most important findings of our paper is that, even though SISR occurs only in the limit as the singular parameter ε → 0, decreasing ε does not increase the coherence of the oscillations due SISR in the FHN model, but rather increases the interval of σ for which coherence occurs. That is, making ε → 0 does not improve the coherence of the oscillations (the coefficient of variation, CV, does not turn to 0; it stays low and constant) when the conditions required for SISR are satisfied. This is in contrast to the excitable model considered in [1], where theoretical and numerical analyses show that when the conditions required for SISR are satisfied (which are essentially the same conditions as required for SISR in the FHN model), the coherence of the oscillations becomes perfect (CV → 0) in the limit as the singular parameter turns to 0 (see Fig. 3(b) in [1]). Therefore, the coherence of the oscillations as the singular parameter turns to 0 depends on the excitable model under consideration, and this might be worth to be investigated in more detail.
Furthermore, our analysis also clarifies the crucial difference between the deterministic parameter setting required for SISR to occur and those required for the occurrence ISR (analyzed in [2]) in the FHN model. We find that it depends on the stability of the (unique) fixed point and whether it is located to either the left or right of the fold point of the critical manifold, whether SISR or ISR can emerge in the same weak noise limit. We thereby produce a unified mathematical setting to analyze both ISR and SISR. This allows us to understand how a neuron could effectively switch from one phenomenon to the other (thus encoding different information) without changing the weak synaptic noise limit required for both phenomena to occur. The rest of the paper is organized as follows: Sect. 2 introduces the version of the neuron model used for our investigation. Section 3 analyzes the deterministic dynamics of the model equation: The singular Hopf bifurcation scenario of the fixed point considered is investigated. We define the critical manifold of the model equation on which we determine the reduced equation of the fast variable together with its evolution timescale. In Sect. 4, we analyze the stochastic dynamics: For the model equation, we obtain the order of the stochastic timescale at which random trajectories escape from the basins of attraction of the stable parts of the critical manifold. In Sect. 5, we analyze the scaling asymptotic limit of the deterministic and stochastic timescales and the consequences on resonance. In Sect. 6, we characterize the limit cycle behavior induced by noise. In Sect. 7, we present and discuss numerical results. We have concluding remarks in Sect. 8.

Model equation
The FHN neuron model (for biophysical details, see [55]) has been extensively used to investigate many complex dynamical behaviors occurring in neuroscience [56][57][58][59][60]. For our investigation of SISR, we consider the following version of the FHN neuron model written on the slow timescale τ as with the deterministic velocity vector fields given by where 0 < ε := τ/t 1 is the timescale separation parameter between the slow timescale τ and the fast timescale t; the fast variable v ∈ R represents the membrane potential (or action potential) variable and the slow variable w ∈ R represents the recovery current (or sodium gating) variable which restores the resting state of the neuron model. Biologically, ε accounts for the slow kinetics of the sodium channel in the nerve cell and controls the main morphology of the action potential generated [61]. d ∈ (0, 1) is a constant parameter, and c > 0 is a co-dimension-one Hopf bifurcation parameter. d and c influence the generation of the action potential, and their ranges have been chosen such that the FHN model best captures some excitableoscillatory behaviors of the biophysically more realistic Hodgkin-Huxley neuron model [55,[61][62][63].
dW τ is standard white noise, the formal derivative of Brownian motion with mean zero and unit variance. This term which is added to the membrane potential variable (v) equation mimics the influence of synaptic noise on the dynamics of the neuron [64]. We note that because of the scaling law of Brownian motion, the scaling σ √ ε of the noise term will guarantee that the amplitude of the noise (σ ) measures the relative strength of the noise term compared to the deterministic term f (v, w) irrespective of the value of ε, when we transform Eq. (1) to its equivalent fast timescale equation [see Eq. (22)].
To analyze SISR, we notice that in Eq. (1), the random term is added to the membrane potential variable equation to mimic the effect of synaptic noise. This is different for CR in the FHN model, where the random term is added instead to the recovery variable to mimic the effect of channel noise (see [3,24,25]). We further note that, just as for the phenomenon of CR, our model equation in Eq. (1) has no periodic perturbation, in contrast to the phenomenon of SR [11].

Deterministic dynamics and its timescale on stable manifolds
The deterministic dynamics and the emergence of a limit cycle solution (relaxation oscillation) in the FHN neuron model equation have been studied extensively [65][66][67][68][69][70], and various types of behavior were detected or confirmed. In this section, we shall establish the deterministic setting in which we want our neuron model to be before a random perturbation is introduced into the equation. In particular, we shall clarify the role played by the relative positioning of the fixed point and the local minimum of the cubic nullcline in the occurrence of either SISR or ISR (analyzed in [2]) in the FHN neuron. It seems to be a new discovery that just (carefully) changing the position of the fixed point with respect to the minimum of the cubic nullcline of the FHN neuron model can lead to very different noise-induced phenomena (in particular SISR and ISR) in the same weak noise limit. We shall also determine the timescale of the slow dynamics of the fast variable v on the stable parts of the critical manifold.
We thus consider the deterministic dynamics corresponding to the behavior of Eq. (1) when σ = 0. At the fixed (or equilibrium) points (v e , w e ) (rest states of the neuron), the variables v(τ ) and w(τ ) reach a stationary state, while the set of equilibrium points defined by the intersection of nullclines as depends on the parameters d ∈ (0, 1) and c > 0. The sign of If = 0, we have two fixed points. This happens if and only if φ = ± π 6 . These limiting cases correspond to saddle-node-type bifurcations during which two fixed points in the case < 0 coincide. If φ = π 6 , then v e 1 and v e 3 will coincide and the two fixed points are given by If φ = − π 6 , then v e 1 and v e 2 will coincide and the two fixed points are given by We denote such a fixed point by (v , w ) and study its stability from the linearization of Eq. (1) at the fixed point chosen. The linearization at (v , w ) is given by with the Jacobian matrix The stability of a fixed point (v , w ) depends on the signs of the trace (tr J i j ) and the determinant (det J i j ) of the Jacobian matrix in Eq. (7). For a fixed point (v , w ) to be stable, we need tr J i j < 0 and det J i j > 0. Since ε, c > 0, we have tr J i j < 0 and det J i j > 0 only if It is important to note here that the condition in Eq. (8) is sufficient, but not necessary for the stability of a fixed point (v , w ). In [2], this was crucial for the analysis of ISR, via the consequences that the stability of a fixed point has on the dynamics of the slow-fast neuron model when the condition in Eq. (8) is violated.
In the present paper, we shall choose c and d such that in Eq. (4) is greater than 0, in which case we have a unique fixed point at (v e , w e ) given by Eq. (5). We also choose c and d such that 1−v 2 e < 0, so that we have a unique and stable fixed point located at (v e , w e ). For the FHN neuron model with a unique and stable fixed point, in the absence of any perturbation as in Eq. (1) with σ = 0 (or in general in the presence of a subthreshold perturbation), the neuron cannot maintain a self-sustained oscillation (i.e., no limit cycle solution). One says in this case that the neuron is in the excitable regime [71]. In this regime, choosing an initial condition in the basin of attraction of (v e , w e ) will result in at most one large non-monotonic excursion into the phase space after which the trajectory returns exponentially fast to this fixed point and stays there until the initial conditions are changed again.
However, as we saw in [2], the existence of a unique and stable fixed point in the FHN model with no perturbation (or even with a subthreshold perturbation) by itself does not necessarily mean that the neuron is in an excitable regime, as in some cases it might have just one stable fixed point coexisting with a stable limit cycle. We consider here the case of a bistable regime (a crucial requirement for the occurrence ISR) consisting of a stable fixed point and a stable limit cycle (see also [72,73]) in which the dynamical behavior is globally different from the one in the excitable regime.
With parameters c and d chosen such that condition Eq. (8) holds for the fixed point (v e , w e ) in Eq. (5), a Hopf bifurcation is the only way through which Eq. (1) (with σ = 0) can exhibit a limit cycle solution. As we shall see later why this should be so, the noise-induced phenomenon of SISR requires that the timescale parameter ε → 0. For the deterministic system in Eq. (1), we shall, for ε → 0, calculate a special value of c = c H that will give us the location of the so-called singular Hopf bifurcation of the fixed point (v e , w e ). A singular Hopf bifurcation in Eq. (1) with σ = 0 (and in planar slow-fast dynamical systems in general) is said to occur if the linearized center manifold system has a pair of singular eigenvalues λ(ε; c) at the Hopf bifurcation point c = c H [74], that is, Clearly, for the fixed point (v e , w e ) in Eq. (5), with the associated linearization in the slow timescale τ given by the matrix J i j in Eq. (7), the Hopf bifurcation occurs at tr J i j = 0, that is, at c = 1 ε (1 − v 2 e ). The eigenvalues of J i j (v e , w e ) at the Hopf bifurcation which tend to infinity as ε → 0. Alternatively, we can also look at Eq. (1) on the fast timescale t = τ/ε [see Eq. (17)] with linearization ε J i j (v e , w e ) and eigenvalues at the Hopf bifurcation given by ] which tend to 0 as ε → 0. On both timescales, the eigenvalues at the Hopf bifurcation are "singular." In short, a singular Hopf bifurcation occurs when the eigenvalues become singular as ε → 0. We shall return to the singular Hopf bifurcation later.
We now need to define the critical manifold of Eq. (1), to determine its stability properties and to obtain the reduced equation describing the evolution of the fast variable v and the deterministic timescale at which v evolves on this manifold. The deterministic critical manifold M 0 defining the phase space of the slow subsystem associated with Eq. (1) [that is, the system we obtain from Eq. (1) in the singular limit ε = 0] is, by solving f (v, w) = 0 for the slow variable w, given by We note that M 0 coincides with the v-nullcline of the neuron model; see the red cubic curve in Fig. 1a, b. This manifold is normally hyperbolic away from the fold points at v = ± 1 satisfying dw dv (± 1) = 0 and naturally splits into three sub-manifolds: The linearized stability of points on M 0 as steady states of the fast sub-system [that is, the system we obtain in the singular limit ε = 0 of Eq. (1) when it is written on the fast timescale t, that is, Eq. (17)], is determined by the Jacobian scalar This shows that points with |v| > 1 are stable, while points with |v| < 1 are unstable. It fol- A trajectory will therefore be attracted to v * − (w) or v * + (w) (depending on initial conditions) and move along these stable manifolds toward the stable fixed point or fold points. Along v * 0 (w), motion is not possible as it is a repulsive branch of M 0 . Motion is however possible on v * 0 (w) when conditions for the so-called canard explosion [75] are satisfied. This situation is not of interest in the context of this work.
For simplicity, we shall denote the basins of attraction of the stable parts of the critical manifold v * − (w) and v * + (w) [that is, the sets of initial conditions that plays the role of the separatrix between these basins of attraction.
From Fenichel's invariant manifold theorem [76], we know that, while the dynamics of the slow subsystem takes place on the critical manifold M 0 , the dynamics of the full system [the dynamics of Eq. (1) when ε = 0, i.e., ε strictly different from 0] takes place not on M 0 itself, but on a perturbed M 0 , the so-called slow manifold denoted by M ε . The theorem states that M ε is at a distance of O(ε) from M 0 and the flow on M ε converges to the flow of the slow subsystem on M 0 as ε → 0. This can be seen in Fig. 1. Here, the black trajectories (with arrows) of the slow sub-system in Fig. 1a evolve on M 0 with fast jumps (indicated by the horizontal double arrows) at the fold points of M 0 . In Fig. 1b, we see the trajectories of the full system for 3 different values of the singular parameter ε. We see that with ε = 0, the trajectories of the full system do not move on M 0 , but instead move at a distance from M 0 , that is, on the slow manifold M ε (not shown). In Fig. 1b, trajectories get closer and closer to M 0 and eventually will coincide with the trajectories of the slow sub-system in Fig. 1a in the limit as ε → 0.
Therefore, because the noise-induced phenomenon under investigation in this paper occurs only in the limit as ε → 0, this suggests to consider the approximation where the full dynamics of Eq. (1) (with ε → 0) takes place on the critical manifold M 0 , which, as just explained, is very close to M ε in the limit ε → 0. In other words, the full dynamics of Eq. (1) (with ε → 0) will coincide with the dynamics of the slow sub-system whose equation is given by Eq. (11), itself obtained in the singular limit ε = 0. Moreover, with this approximation, the basins of attraction of the stable parts of M ε can be assumed to coincide with B(v * − (w)) and Fig. 1 In a, black trajectories moving on the critical manifold M 0 (red cubic curve) represent the singular solution of Eq. (11) (which coincides with the solution of Eq. (1) in the limit as ε → 0). The single arrows indicate slow motion, and the double arrows indicate fast motion (horizontal jumps at the fold point). In b, trajectories of the full system with ε > 0 move on M ε (not shown but at a distance of O(ε) from M 0 ) into the stable fixed point located at the intersection of M 0 and the w-nullcline, represented by the green line. As ε becomes smaller and smaller, the dynamics of the full system get closer and closer and eventually coincide with the dynamics of Eq. (11) on the stable parts of M 0 The reduced dynamics along v * − (w) and v * + (w) for the fast variable v is governed by Eq. (11). This reduced equation is obtained by Combining this with the equation for w in Eq. (1), we eliminate the slow variable to get Equation (11) is an ODE with an algebraic constraint w = v− v 3 3 . Therefore, we have a differential-algebraic equation whose initial conditions must satisfy this constraint for solutions to exist. In other words, the dynamics on the slow timescale is described by the reduced equation [in Eq. (11)] defined by the projection of the slow dynamics onto the tangent space of the critical manifold. We note that Eq. (11) becomes singular at v = ± 1 (the fold points) where M 0 loses its normal hyperbolicity. In particular, the existence and uniqueness theorems for ODEs do not apply at these points. This is the reason why a trajectory moving along v * − (w) and v * + (w) jumps away horizontally when it reaches the fold points located at (± 1, ± 2 3 ). See the black trajectory with arrows in Fig. 1a.
In view of the properties of this approximation, we shall henceforth talk of the dynamics of the full system on the critical manifold M 0 . Thus, in the limit ε → 0, the dynamics of the full system in Eq. (1) is described by Eq. (11). For concreteness, we shall choose ε = 0.0001 in the numerical plots, except in Figs. 5 and 8, where we vary ε to see how it affect the phenomenon.
With σ = 0, > 0 and 1 − v 2 e < 0, a trajectory with initial conditions in B(v * − (w)) will be attracted to v * − (w) and slowly (as ε → 0) move along it downwards toward the unique and stable fixed point (v e , w e ) in Eq. (5), where it stops. This behavior is shown in Fig. 2d for the blue trajectory with the initial conditions at (v 0 , w 0 ) = (− 2.0, 0.25). If the initial conditions are now chosen in B(v * + (w)), as in Fig. 2d for the black trajectory with the initial conditions at (v 0 , w 0 ) = (2.0, − 0.25), then this trajectory is also attracted to v * + (w) and moves slowly (as ε → 0) along it upwards toward the fold point at 1, 2 3 , where it is forced to jump to v * − (w) and then moves again down toward the stable fixed point (v e , w e ), where it stops.
We note that as a trajectory moves along on v * − (w), and when the unique fixed point (v e , w e ) satisfies the condition in Eq. (8), then it will stop at this fixed point and not move onto the fold point at (− 1, − 2 3 ), from where it will have jumped to the other stable manifold v * + (w). This is observed again in Fig. 2d with the blue and black trajectories. There, the trajectories stop at the stable fixed point located at (v e , w e ) = (− 1.003988, − 0.666651) and do not move on to the jump point (− 1, − 2 3 ), which is located to the right of (v e , w e ); notice 1 − v 2 e < 0. Thus, the fold point (− 1, − 2 3 ) will not be reached by the trajectories of the full system. Equation (11) together with the equations of the stable manifolds v = v * − (w) and v = v * + (w) specifies the dynamics of the fast variable v arising on the deterministic slow timescale of O(ε −1 ). We note that this timescale becomes slower and slower as ε → 0.
We now return to our discussion of the singular Hopf bifurcation. We state and use the theorem by Krupa and Szmolyan [75] to determine the value of the singular Hopf bifurcation c H for Eq. (1) together with its criticality. Consider a general slow-fast dynamical system on the fast timescale t given in the normal form: where ε is the timescale separation parameter, α is the Hopf bifurcation parameter, and the functions h i are given by We define several computable constants, abbreviating (0, 0, 0, 0) = (0) in the definitions: Note that all a 1 to a 4 only depend on the partial derivatives with respect to x. Next we define another constant: Theorem Suppose (x, y) = (0, 0) is a generic folded singularity for α = 0 with normal form Eq. (12). Then there exist ε 0 > 0 and α 0 > 0 such that for 0 < ε < ε 0 and |α| < α 0 , in a suitable neighborhood of the origin, system Eq. (12) has precisely one equilibrium point (x e , y e ) with (x e , y e ) → (0, 0) as (α, ε) → (0, 0). Furthermore, the equilibrium point (x e , y e ) undergoes a Hopf bifurcation at α H with: (a) The Hopf bifurcation is non-degenerate for A = 0, supercritical for A < 0, and subcritical for A > 0.
To apply the theorem, we write Eq. (1) with σ = 0 on the fast timescale t = τ/ε as The fold points of M 0 are at its extrema: Since we want our unique fixed point (v e , w e ) in Eq. (5) to be stable, we again choose c and d such that Eq. (8) holds. It is useful to see that this will imply that our unique fixed point (v e , w e ) will be located on a decreasing part of the critical manifold M 0 , that is, either on v * − (w) or on v * + (w). It is easy to check that the fold point p 1 satisfies the conditions necessary for a generic folded singularity [77]. We therefore direct our interest to the fold point is sufficient to make it stable. We shall consider this deterministic setting for our analysis of SISR.
It is important to point out that, even though the fold point , the relative position of these points is crucial for the occurrence of either SISR or ISR as analyzed in [2]. To see why this relative positioning (positioning the fixed point (v e , w e ) either to the left or right of the fold point p 1 = (− 1, − 2 3 )) is important for SISR, consider the opposite situation, that is, assume that the fixed point (v e , w e ) lies to the right of the fold point With the fixed point (v e , w e ) lying on v * 0 (w), it can be either stable or unstable (recall that Eq. (8) is only a sufficient but not a necessary condition for stability of a fixed point). If (v e , w e ) is unstable (after losing stability through a Hopf bifurcation), then even in the complete absence of noise, a trajectory with initial conditions in B(v * − (w)) will be attracted to v * − (w) along which it moves toward the fold point p 1 = (− 1, − 2 3 ), at which point it jumps horizontally and avoids the unstable fixed point (v e , w e ) located on the unstable branch v * 0 (w) as it finally drops on the stable branch v * + (w). Along v * + (w), it moves toward the fold point where it again jumps horizontally to drop back on the stable branch v * − (w), from where the whole cycle repeats again. This behavior is shown in Fig. 2b where the blue closed trajectory with initial conditions located at (v 0 , w 0 ) = (− 2.0, 0.25) jumps at the fold points located at the extrema (± 1, ± 2 3 ) of M 0 . This will lead to a deterministic limit cycle (relaxation oscillation) around the unstable fixed point. This is precisely the deterministic setting we do not want to have, as the phenomenon of SISR requires the emergence of a limit cycle behavior due only to the introduction of noise into the system and not because of the presence of a Hopf bifurcation.
On the contrary, as we have shown in [2], the phenomenon of ISR requires that the fixed point (v e , w e ) is located to the right of the fold point p 1 [that is, on the unstable manifold v * 0 (w)] and, most importantly, has to be stable. This deterministic setting requires that ε > 0, a condition which, as we shall see, does not allow SISR to occur, but which is crucial for the occurrence of ISR. Therefore, for SISR to occur in Eq. (1) where α := 1 − d − 2c 3 . We compare Eq. (18) and the standard form in Eq. (12), and we compute the relevant parameters a i defined in Eq. (14) to get We use Eq. (16) to get the singular Hopf bifurcation value c H of our neuron model equation as It is important to note that Eq. (20) gives the location of the Hopf bifurcation only in the limit as ε → 0. We choose d = 0.5 (this value is maintained through out the paper), and we choose ε very small, for example ε = 0.0001, to have  Fig. 2b, that is, when c = 0.745 < c H . Here, the fixed point is located to the right of the p 1 and it is unstable. We recall that the situation where c < c H is not of interest to us because in this regime, we already have a limit cycle due to a supercritical singular Hopf bifurcation at c H .  8). This is precisely the state in which we want our deterministic neuron in Eq. (17) to be before a random perturbation is added to its fast variable (v) equation, and we want to see how this noise can induce a limit cycle behavior (coherent spike train) that the deterministic model cannot exhibit, that is, the phenomenon of SISR.

Stochastic dynamics and timescale of escape processes from the stable manifolds
So far, we have considered the dynamics of Eq. (1) when there is no noise, i.e., when σ = 0. Now we switch on the noise. In order to use the corresponding theory [78,79], we shall derive the necessary formulae for the FHN neuron model in detail. Due to the diffusion that results from the presence of noise, random trajectories may then eventually escape from B(v * − (w)) or B(v * + (w)) before the fixed point or fold points are reached. To understand the escape processes from the basins of attraction and to quantify the stochastic timescale of these escape processes, we transform Eq. (1) to its equivalent fast timescale equation [which evolves on a timescale of O(1)] by rescaling the time as t = τ/ε. We note that the noise term (and not only the variables v τ and w τ ) has to be taken into consideration during the rescaling. The noise term is rescaled according to the scaling law of Brownian motion. That is, if W τ is a standard Brownian motion, then, for every λ > 0, λ −1/2 W λτ is also a standard Brownian motion (i.e., the two processes have the same distribution [80]). We therefore have Eq. (1) written on the fast timescale as follows: We also notice that the noise term is now independent of ε, and hence, σ will measure the relative strength of the noise term compared to the deterministic term irrespective of the value of ε. Furthermore, the form of Eq. (22), unlike Eq. (1), will avoid singularity problems when we are in the limit ε → 0.
In the singular limit ε → 0, the timescale separation between v t and w t becomes larger and Eq. (22) becomes a singularly perturbed ODE, with the equation for w t reducing to dw t = 0. This indicates that w t is frozen on the O(1) fast timescale and, on this timescale, the dynamics is governed by the fast equation in Eq. (22), where w t enters merely as a fixed parameter. Thus, Eq. (22) in this limit becomes a 1-D stochastic differential equation of the form where viewed as a function of v with w nearly constant, is a double-well potential. The two local minima are located at the value of v where v * − (w) and v * + (w) intersect the horizontal line w = k with k ∈ − 2 3 , 2 3 , and a local maximum at the intersection with the unstable branch v * 0 (w); see Fig. 3. We now define the quantities where from Eq. (9), we use Viète's trigonometric expressions of the roots in the three-real-roots case (w 2 ≤ 2/3) to get ⎧ ⎨ ⎩ v * − (w) = 2 cos 2 3 π + 1 3 arccos(− 3 2 w) , v * 0 (w) = 2 cos − 2 3 π + 1 3 arccos(− 3 2 w) , v * + (w) = 2 cos 1 3 arccos(− 3 2 w) .
The graphs of v * − (w) and v * + (w) in Eq. (26) are strictly monotonically decreasing for values of w ∈ − 2 3 , 2 3 . In each case, U ± (w) is the potential difference (energy barrier function) between the local maximum v * 0 (w) and a local minimum v * ± (w), respectively. See Fig. 4.
In order to understand the escape processes of the trajectories of Eq. (23) from the basins of attraction , and to estimate the stochastic timescale at which these escape processes take place, one may picture these escape processes as the motion of a particle in a double-well potential under the influence of a stochastic perturbation. That is, one could view the escape of a trajectory of Eq. (23) from the basin of attraction B(v * − (w)) into B(v * + (w)) as the escape of a particle from the minimum of the left well at v * − (w) into the minimum of the right well at v * + (w) of the double-well potential U (v, w) and vice versa. For the neuron model equation under consideration, there is a perfect analogy between the escape process of a trajectory from the basins of attraction B(v * ± (w)) and the escape process of a particle from those minima of the double-well potential U (v, w).
With this picture in mind, we can conveniently estimate the escape rate E r − of a trajectory from B(v * − (w)) into B(v * + (w)) and hence the stochastic timescale 1/E r − at which the trajectory escapes. The escape rate E r + of a trajectory from B(v * + (w)) back into B(v * − (w)) and the stochastic timescale 1/E r + at which it occurs can be estimated analogously to E r − because of the symmetry between the two escape processes.
To estimate E r − , we consider our particle obeying the dynamics of Eq. (23), initially at the minimum v * − (w) of the left well and wanting to escape to the minimum v * + (w) of the right well of the potential U (v, w). This setting corresponds to the situation in Fig. 3a. For this escape event to occur, the particle should be able to go over the potential barrier U − (w) given in Eq. (25).
With σ = 0, the probability density p(v, t) of finding the particle at position v at time t in the double-well potential U (v, w) obeys the Fokker-Planck equation (FPE) [81] corresponding to the Langevin equation in Eq. (23); that is, with initial and natural boundary conditions chosen, respectively, as where v 0 and t 0 are the initial position and starting time, respectively. We write Eq. (27) in the form of the continuity equation [81] as where the so-called current or probability flux j (v, t) is given by as w is merely a constant here. With our particle at the left minimum v * − (w) of the potential U (v, w), an escape event into the right minimum at v * + (w) occurs if the noise is strong enough to push the particle over the potential barrier U − (w). If the particle is initially (t = 0) at position v * − (w), then from the initial boundary condition in Eq. (28), we have the probability density given by When the quantity U (v) σ 2 in Eq. (30) is large (e.g., in the limit as σ → 0), we expect p(v, t) to be concentrated around the minimum v * − (w) of the left potential well. This implies that the escape over the potential barrier U − (w) becomes a rare event. We are precisely interested in the probability of this rare event. However, if we allow the particle to interact with the noise for a sufficiently long time (in other words, if we integrate Eq. (23) in a sufficiently large time interval), then this weak noise may eventually push the particle far away from v * − (w) into the right minimum at v * + (w) by going over U − (w). Thus, making U (v) σ 2 → ∞, we make the probability flux j in Eq. (30) infinitely small and the probability density p(v, t) becomes almost time-independent. In this case, we obtain a stationary To obtain the escape rate of the particle, we therefore have to integrate Eq. (32) from the initial position of the particle [that is, at the minimum v * − (w)] to some point β 1 beyond the maximum v * 0 (w) (β 1 > v * 0 (w)) of the potential U (v, w), i.e., From the natural boundary condition in Eq. (28), we have that p(β 1 , t) ≈ 0 and we obtain the constant probability flux j 0 from Eq.(33) as When the particle is at the minimum v * − (w), and the quantity U (v) σ 2 → ∞ as σ → 0, and in Eq. (30), we have j ≈ 0 which means that p(v, t) exp 2U (v) σ 2 = β 2 , where β 2 is a constant. We then obtain the invariant probability distribution as The constant β 2 can be given a convenient arbitrary value; we choose , so that the invariant probability density of the particle when it is at the minimum v * − (w) is now given by We now consider two points β 3 < v * − (w) < β 4 in the immediate neighborhood of the minimum at v * − (w). The probability p 0 of finding the particle in the interval β 3 , β 4 is calculated by integrating over the invariant probability density p i (v, t), i.e., As the value of the invariant density p i (v, t) in Eq. (36) becomes small away from the minimum at v * − (w) is very small, then we do not need to know the values of β 3 and β 4 to evaluate the integral in Eq. (37).
The escape rate E r − can then be obtained by noting that Eq. (34) gives the conditional probability of escape per unit time, given that the escaping particle is initially at the minimum v * − (w) of the potential. With and Eqs. (34) and (37), we obtain Since the integrands in the first and second integral of Eq. (39) are peaking at the minimum v * − (w) and the maximum v * 0 (w), respectively, the next approximation consists in Taylor expanding U (v) to second order about these peak values while noting that at the local extrema of the potential U (v), we have Noting that U (v * − (w)) and U (v * 0 (w)) are, respectively, the dominant terms in the two integrals in Eq. (39), we insert the Taylor expansions into Eq. (39) and use the fact that ∞ −∞ exp (−αx 2 )dx = π α to obtain the escape rate E r − as By the same analysis, we get the escape rate E r + in the reverse direction as where the prefactors μ ± are given, respectively, by the square roots of the product of the curvatures of the potential at its local extrema as Here, we could extend the domain of the integration in Eq. (39) to the real line because the contributions of points away from the minimum point v * − (w) and the maximum point v * 0 (w) are exponentially small. For this step, it is important that the original domain of integration extended to some point β 1 beyond the maximum v * 0 (w), so that we integrate over some interval containing v * 0 (w) in its interior. For bistable systems, the inverse of the escape rate from a basin of attraction is the escape time from this basin of attraction [78,82]. The stochastic timescales at which trajectories escape from the basins of attraction B(v * ± (w)) are therefore of O 1

Asymptotic matching of timescales and resonance
We still need to do some work, as even though in the previous section we derived for Eq. (22) the stochastic timescale of the escape processes in 1-D by neglecting the slow equation in the limit ε = 0, ε is only very small, but not exactly 0 in the full dynamics of the neuron [that is, the dynamics of Eq. (22) when ε = 0]. However, as we already pointed out in Sect. 3, from Fenichel's invariant manifold theorem, the limit ε → 0 used to obtain the stochastic timescales in the previous section is a sufficiently good approximation because under this limit, the basins of attraction of the stable parts of the slow manifold M ε coincide with B(v * ± (w)). In this section, where we consider the full dynamics, this approximation will continue to be valid. Now consider 0 < ε 1 (positive but very small), that is, we allow w in Eq. (22) to move as well. We choose c > c H so that no limit cycle can appear due to a singular Hopf bifurcation. If the stochastic timescales of the escape processes are much longer (this can happen with extremely weak noise) than the deterministic timescale of the reduced equation in Eq. (11) on the manifolds v * ± (w), then the trajectory has no time to escape and is most likely to stay in the basins of attraction B(v * ± (w)) for a very long time with no or very rare transitions between B(v * − (w)) and B(v * + (w)). In this case, the trajectory is forced to move "almost" deterministically on v * ± (w) toward the stable fixed point (v e , w e ) where it sticks.
On the other hand, when the stochastic timescales are much shorter (this can happen when the noise is strong) than the deterministic timescale, the noiseinduced transitions are very frequent and incoherent. In this case, the neuron's dynamics can only be captured by its invariant density. This is an immediate consequence of a strong noise with which no coherent structure can emerge.
An interesting case occurs when the escape events from B(v * ± (w)) and motions along v * ± (w) have comparable timescales. With the stochastic timescales given by the inverse of Eqs. (41) and (42), a trajectory with initial conditions in B(v * − (w)) may escape into B(v * + (w)) and conversely at some specific points of w with probabilities close to 1 if some suitable scaling limit conditions are imposed. With these conditions, the reduced equation in Eq. (11) may not be valid along v * − (w) all the way down to the fixed point (v e , w e ) and may neither be valid along v * + (w) all the way up to the fold point p 2 = (1, 2 3 ). This equation will only govern the motions on these stable manifolds until some well-defined points, w − ∈ v * − (w) and w + ∈ v * + (w) where the trajectory, respectively, escapes from B(v * − (w)) and B(v * + (w)). Indeed, as soon as w − is reached, the stochastic timescale becomes comparable to the deterministic timescale of the motion along v * − (w) and an escape event instantaneously happens at w − and likewise for w + . Therefore, we can only expect a coherent spike train with a weak noise when the neuron can match the deterministic timescale to the stochastic timescale only when w = w ± , but not at other points. Thus, we need For this, we must therefore choose a suitable double scaling limit: (σ, ε) → (0, 0), such that for some finite Φ ∈ (Φ min , Φ max ), with Φ min and Φ max to be specified later. Eq. (45) implies that U − (w − ) → U + (w + ), and therefore, if a trajectory escapes from B(v * − (w)) at w − , then at w + = − w − , we should expect an escape from B(v * + (w)). Provided that the equations in Eq. (46) have solutions on the accessible parts of the stable manifolds v * ± (w), the jump points w ± are the unique solutions of The uniqueness follows from the monotonicity of the energy barrier functions U − (w) for w ∈ (− 2 3 , 0) and U + (w) for w ∈ (0, 2 3 ); see Fig. 4. The matching of these timescales implied in Eq. (45) is precisely the resonance mechanism in standard stochastic resonance [12].
We note that we must have w e < w − , that is, when the point w = w − can be reached by the slow flow of Eq. (11) along v * − (w). Otherwise, the trajec-tory of Eq. (11) will stick to the stable fixed point at w e . Since U − (w) is monotone, this will occur when Φ > Φ min = U − (w e ) (since w e is the lowest attainable point on v * − (w)). On the other hand, we should also have w − < w + , which is violated when Φ ≥ Φ max = 3 4 (w − = w + = 0 at Φ = Φ max ); see Fig. 4. The limit cycle behavior (coherent spike train) emerges only if we choose σ and ε sufficiently small such that Φ ∈ (Φ min , Φ max ) = U − (w e ), 3 4 . This interval is not empty because w e = − 0.666651 < 0.
With the asymptotic limits in Eq. (45), that is, Φ ∈ U − (w e ), 3 4 , a truly deterministic limit cycle behavior could emerge even though c > c H (a regime in which the zero-noise dynamics cannot display a limit cycle). The phase portrait of this noise-induced limit cycle behavior is composed of the two portions of the critical manifolds v * ± (w) between the jump points w = w − and w = w + , together with the horizontal lines joining these manifolds at w − and w + . See the phase portraits in Fig. 7; they are different from the one in Fig. 2b, where the jumps from the stable manifolds can only occur at the fold points (± 1, ± 2 3 ). We note also that on v * − (w), w − is a jump point, while w + is a drop point. On v * + (w), w − is a drop point, while w + is a jump point.

Characterization of the limit cycle
First, we note that while the deterministic characteristics (e.g., the period) of the limit cycle behavior obtained are controlled by the value of Φ, the degree of its coherence is controlled by σ and ε. It can thus be made as high as desired by choosing σ and ε very small, provided that Φ ∈ U − (w e ), 3 4 . Since a trajectory spends most of the time on the stable manifolds v * − (w) and v * + (w) (due to the slow motion on the timescale of O(ε −1 ) with ε → 0), asymptotically, the period of the limit cycle will be the sum of the time it spends on these manifolds. That is, it is the sum of the time it takes to go from the drop point w + to the jump point w − on v * − (w) and from the drop point w − to the jump point w + on v * + (w). This period, T (d, c, σ, ε), can be obtained explicitly as the sum of the integrals over the reduced equation in Eq. (11) from v * − (w + ) to v * − (w − ) and from v * + (w − ) to v * + (w + ), that is, The integrals in Eq.(47) certainly exist for any w ± ∈ − 2 3 , 2 3 , thanks to the boundedness of the graphs of v * − (w) and v * + (w) for w ∈ − 2 3 , 2 3 and thanks to the boundedness and continuity of the integrand , v * + − 2 3 = 1, 2 in the first and second integral, respectively.
From the boundaries of integration in Eq. (47), one can see that the period of the limit cycle created by noise has a non-trivial dependence on σ and ε through its dependence on the points w ± , themselves connected to σ and ε as in Eqs. (45) or (46). Thus, the period of the limit cycle can be controlled by σ for a fixed ε without significantly affecting its coherence, provided Eq. (44) holds. From the asymptotic theory above, it is determined [from Eq. (46)] that a trajectory escapes at w − = − w + = − 0.432 for ε = 0.0001 and σ = 0.005, which gives a period of T ≈ 1.6396 for c = 0.760 > c H .
We now investigate what happens in the neighborhood of the singular Hopf bifurcation value, more precisely how the emergence of the limit cycle behavior is affected as c approaches c H from above (c → c + H ). The case in which c H is approached from below (c → c − H ) is not interesting because in this case, a limit cycle already exists due to a supercritical singular Hopf bifurcation, while we are interested in a limit cycle behavior due only to the presence of noise. As before, we take the double limit (σ, ε) → (0, 0) such that The upper bound of this open interval is fixed and does not depend on the choice of the parameter c. But the lower bound depends on w e which in turn depends on the parameter c as in Eq. (5). Therefore, the effect of c → c + H on the emergence of the limit cycle behavior could be "observed" just by looking at the behavior of the following limit: If this limit decreases as we get closer and closer to c H , thus increasing the length of the open interval in Eq. (48), then the emergence of the limit cycle behavior is facilitated in the sense that for a fixed ε, weaker and weaker noise could trigger the limit cycle behavior.
On the other hand, if this limit continuously increases, thus decreasing the length of this open interval, then the emergence of the limit cycle behavior is inhibited in the sense that for a fixed ε, we are now having a thinner and thinner interval of the noise amplitude for which the behavior can occur. And if this interval length shrinks to zero, then the limit cycle behavior simply disappears regardless of the value of σ we choose for a fixed ε.
With d = 0.5, we see [from Eq. (5)] that v e → −1 as c → c + H , making w e → − 2 3 . Then, in Eq. (49), U − (w e (c)) → 0 as w e → − 2 3 ; this limiting behavior is readily seen in the red curve in Fig. 4. Thus, as we approach the singular Hopf bifurcation value, the limit cycle behavior due to noise is facilitated; see Fig. 9 and compare with Fig. 6.

Numerical simulations and discussion
For the purpose of numerical computations, the term 1 2 σ 2 in 1 2 σ 2 log e (ε −1 ) will be absorbed into σ so that the matching of the timescales in Eq. (45) leading to the emergence of a coherent spike train will require that With Eq. (50), we easily calculate the minimum (σ min ) and the maximum (σ max ) of the noise amplitude required to induce a coherent spike train in the neuron's activity when ε → 0 as With the fixed point (v e , w e ) = (− 1.003988, − 0.666651) located to the left of the fold point p 1 = (− 1, − 2 3 ), we therefore must have w e > − 2 3 (even though w e is very close to − 2 3 , the inequality has to be strict). Moreover, the energy barrier function U − (w) at the fold point p 1 is U − (w = − 2 3 ) = 0 (again seen from the red curve in Fig. 4). These facts imply that U − (w e ) > U − (− 2 3 ) = 0. Therefore, in the double limit (ε, σ ) → (0, 0) we have an incoherent spike train if a very small but non-empty interval, since w e > − 2 3 and U − (w) is a monotonically increasing function of w ∈ (− 2 3 , 0). Inserting the value of w e in U − (w e ), we get the interval in Eq. (52) for which only a rare and incoherent spike train emerges in the limit (ε, σ ) → (0, 0) as . It follows from Eq. (51) that choosing the noise amplitude σ such that σ ≤ σ min = 1.46667 × 10 −5 / log e (ε −1 ) for a fixed ε → 0 will lead only to incoherent spike trains. For σ ≥ σ max , we of course have incoherent spike trains as the noise amplitude is already too strong. Now we present and discuss the numerical simulations carried out to illustrate the theoretical analysis. To numerically generate the Gaussian white noise W t , we start from two random numbers b 1 and b 2 which are uniformly distributed on the unit interval [0, 1] and, with the Box-Muller algorithm [83], we generate a standardized Gaussian distributed sequence. We use the Euler algorithm [84] to integrate Eq. (22) with initial conditions at (v 0 , w 0 ) = (− 2.0, 0.25); the numerical scheme is To quantify the regularity of the spiking activity of the neuron, we use the first two moments of the interval between excitations [the interspike interval (ISI)] of [3]. We determine the spike occurrence times by the upward crossing of the membrane potential variable v past the spike detection threshold of v th = 0. Then, we calculate the regularity of the spike train using the coefficient of variation (CV) defined as where ISI and ISI 2 represent the mean and the mean squared interspike intervals, respectively. Biologically, CV is important because it is related to the timing precision of the information processing in neural systems [85]. For a Poissonian spike train (rare and incoherent spiking), CV = 1. If CV < 1, the sequence becomes more coherent and CV vanishes for periodic deterministic spikes, for example in the deterministic limit cycle regime of Eq. (17); see the perfect coherence of the spike train in Fig. 2a. In the case CV < 1, large excursions of trajectories in the phase space can be interpreted as motion on a "stochastic limit cycle" [86]. CV values > 1 correspond to a point process that is more variable than a Poisson process.
In Fig. 5a, we have the mean interspike interval ISI and its standard deviation shown as error bars as a function of the noise amplitude σ for 4 different values of the timescale parameter ε. The mean interspike interval decreases with increase in the noise if σ / ∈ σ min , σ max and stays almost unchanged if σ ∈ σ min , σ max , for each value of ε.
In Fig. 5b, we show the CV as a function of σ for the same 4 values of ε as in Fig. 5a. For very small values of σ (i.e., σ ∼ 10 −7 < σ min for ε = 0.0001, σ ∼ 10 −6 < σ min for ε = 0.0005, σ ∼ 10 −5 < σ min for ε = 0.001, and σ ∼ 10 −4 < σ min for ε = 0.01, with each σ min (ε) given by Eq. (51)), the firing has the character of a Poisson process since CV is close to 1 (see also the large error bars). The activity in these cases represent rare and incoherent spike trains. For σ ∼ 10 −7 and ε = 0.0001 for example, the noise is just too weak to induce any oscillation, and because c = 0.76 > c H , the trajectory just evolves and stops at the stable fixed point at (v e , w e ).
On the other hand, when σ begins to be large, i.e., when σ ∼ 10 −2 > σ max for all 4 values of ε, with σ max (ε) given by Eq. (51), the spike train also starts losing coherence as CV starts to increase rapidly.
When 10 −6 < σ < 10 −2 for ε = 0.0001, 10 −5 < σ < 10 −2 for ε = 0.0005, 10 −4 < σ < 10 −2 for ε = 0.001, and 10 −3 < σ < 10 −2 for ε = 0.01, we have CV ≈ 0.2, indicating a high level of coherence of the spike train consistent with theory. In the value range of ε satisfying Eq. (50), a smaller value does not increase the coherence of the spike train, but instead increases the interval in which the noise amplitude must lie to be able to induce SISR. This is not the case in [1], where CV → 0 as the singular parameter turns to 0 and the interval of the noise for which SISR occurs stays almost unchanged. We can also see in Fig. 5a that in these intervals of σ , for each ε, the short error bars do not almost change in height indicating that the high coherence of the spike train is not affected by noise when σ lies in these intervals for each value of ε.  Fig. 9 having a frequent and coherent spike train with the same noise strength and same integration time interval but with c = 0.756 which is closer to c H than in Fig. 6 10 −6 . As predicted, we have a rare and incoherent spike train. Recall that the emergence of a coherent spike train requires the matching of the timescales such that Eq. (50) holds. One easily checks that this is not the case for Fig. 6, where Eq. (52) holds instead. Therefore, even though σ and ε are small, we have only a rare and incoherent spike train for a very large integration time interval.
In the double limit σ → 0 and ε → 0 such that σ log e (ε −1 ) ∈ U − (w e ), 3 4 , the spike train abruptly changes to a frequent and coherent one (limit cycle behavior; see Fig. 7), consistent with the theoretical analysis. We further note that the simulations in Fig.  7 run for a time interval 5 times shorter than those in Fig. 6. But because of SISR in Fig. 7, we have a more frequent and coherent spiking than in Fig. 6. In Fig. 7a, b, with σ = 0.005, the jump points are at w − ≈ −0.585 ± 0.075 and w + ≈ 0.591 ± 0.025, with ISI ≈ 1.9348. These jump points are a little later than those predicted by the theory and ISI is within ∼ 18% of the period from the theory.
The simulations in Fig. 7 show that the stronger the noise is (of course the noise amplitude always has to satisfy Eq. (50) for a given ε), the further away are the jump points w − from the fixed point (v e , w e ) and w + from the fold point p 2 = (1, 2 3 ). The trajectories never succeed in attaining (v e , w e ) and (1, 2 3 ) as they are systematically kicked out of B(v * − (w)) and B(v * + (w)) before these points are reached. Moreover, increasing the noise amplitude has the effect of decreasing the period of the oscillations while remarkably keeping the high degree of coherence of the spike train.
The simulations in Fig. 8 show how the spiking frequency for a fixed weak noise amplitude varies with the singular parameter ε. As the value of ε increases, the spiking frequency also increases. The reason for this lies in the fact that the deterministic timescale on which trajectories move on the stable manifolds v * − (w) and v * + (w) is of the order of ε −1 . Therefore, as ε becomes larger and larger (for example, 0.0005, 0.001, 0.002), the time spent by the trajectories on these stable manifolds becomes shorter and shorter, thus decreasing the period of oscillations for a given integration time interval.
In Fig. 9, we verify that as c → c + H the limit cycle behavior can be induced by a weaker noise. We choose c = 0.756 > c H , a value at which the deterministic neuron cannot display a limit cycle as in Fig. 2c and at the same time is closer to c H than in the case of Fig.  6, where the same weak noise amplitude (σ = 1.55 × 10 −7 ) cannot induce a coherent spike train. At c = 0.756, this noise intensity produces a coherent spike train; see Fig. 9. This has been predicted theoretically above by Eq. (49) where U − (w e (c)) → 0 as c → c + H . This increases the length of the open interval in Eq. (50) which means that weaker and weaker noise amplitudes have enough strength to induce a coherent spike train. The limit cycle behavior in this case is much closer to that in Fig. 2a due to proximity to c H .
Finally, we verify the robustness of SISR to parameter tuning. We show that the effect does not require tuning of the bifurcation parameter and can be realized for any value c even much farther away from c H . In contrast, for CR to occur in the FHN model, the bifurcation parameter has to be in the immediate neighborhood of the Hopf bifurcation [3,24].
To check the robustness of SISR, we take c = 1.5, which is two times larger than c H . The result is shown in Fig. 10. Note that for this value of c one needs a stronger noise intensity to induce the oscillations. We see a coherent almost periodic spike train, with clearly defined jump points which are significantly higher than the fixed point (v e , w e ) on the left branch of the critical manifold M 0 and lower than the fold point at p 2 = (1, 2 3 ) on the right branch of M 0 .  Fig. 9 Here, we have ε = 0.0001 and c = 0.756 which is closer to c H than in Fig. 6. We observe that the same noise amplitude producing only rare and incoherent spike train in Fig. 6 produces rather a frequent and coherent spike train in Fig. 9 for the same integration time interval. Thus, closeness to the singular Hopf bifurcation value facilitates SISR by extending the interval of σ for which a coherent spike train emerges

Concluding remarks
In this paper, we have analyzed the effects of synaptic noise on the dynamics of a spiking neuron model with a strong timescale separation (ε → 0) between its dynamical variables. First, we analyzed the deterministic dynamics (σ = 0) and determined the locations of the unique and stable fixed point and the singular Hopf bifurcation value c H . From the slow-fast structure of the neuron model, we obtained the deterministic timescale at which the slow trajectories move on the stable parts of the critical manifold.
We set our neuron model in a regime where deterministic dynamics cannot display a limit cycle (c > c H ) and added white noise (σ > 0) to the fast membrane potential variable (to model the effect of synaptic noise). By considering the singular limit on the fast timescale, we obtained the stochastic timescales of the escape processes from the basins of attraction of the stable manifolds of the model equation.
If the deterministic timescale of the motion of a trajectory on the stable parts of the critical manifold is always matched to the stochastic timescale of the escaping trajectory at specific points on these stable manifolds, then weak noise can induce a transition to a coherent spike train (limit cycle behavior) which has no deterministic counterpart, that is, the following conditions should be satisfied: such that σ log e (ε −1 ) ∈ U − (w e ), 3 4 , σ log e (ε −1 ) = O(1).
We saw that, within the range of values of ε for which Eqs. (55) and (56) hold, a smaller value does not significantly increase the coherence of the spike train, but instead increases the interval of the noise amplitude which could induce SISR.
Moreover, the coherent spiking activity due to noise is shown to be robust to parameter tuning as it persists for values of c c H . This effect is quite different from CR, where the bifurcation parameter has to be in the immediate neighborhood of the Hopf bifurca-tion threshold, making it sensitive to parameter tuning. SISR could therefore, from a biological point of view, provide a mechanism through which real neurons may convert irregular signals into a frequency-encoded periodic output irrespective of the state of other parameters.
Summarizing and comparing the findings of [2] and the present paper, we first used Fenichel's theorem to approximate the dynamics for a small timescale separation ε by the dynamics on the stable part of the slow manifold, with jumps to the other stable branch occurring at those points where the slow manifold loses its stability. These fold points are given by the local extrema of the cubic polynomial that defines the slow manifold. The crucial point then is whether the trajectory on the slow manifold first reaches the stable fixed point or the fold point. In the first case, introducing noise can have the effect of causing an escape before the stable fixed point (located to the left of the fold point) is reached. That is self-induced stochastic resonance requires ε → 0 as investigated in this paper. In contrast, in the second case, noise can induce a transition to the basin of attraction of the fixed point (now located to the right of the fold point) before the fold point is reached. This is inverse stochastic resonance, which requires ε > 0, as studied in [2]. In both phenomena, the fold point corresponds to the minimum of the slow manifold.
However, both SISR and ISR require a weak synaptic noise limit to occur. This could mean that in a weak synaptic noise limit, a neuron could spontaneously switch from encoding information generated during SISR to encoding information generated during ISR, just by changing the membrane voltage of its quiescent state (i.e., by changing the position of the fixed point). Our results could explain experimental observations where very distinct spiking behaviors in physiologically identical neurons may emerge even when these neurons are excited by the same strength and type of stimulus.
Finally, it is worth noting that the type of bifurcation leading to deterministic spiking (Hopf bifurcation for a type II neuron or saddle-node bifurcation for a type I neuron) does not affect the mechanisms behind SISR. We see that SISR occurs if an excitable system is driven by small noise, the escape times of the trajectories from some part of the phase space are exponentially distributed, and their transition rate is governed by an activation energy. If the system is placed out of equilibrium, the parameters are chosen such that the system relaxes slowly to the equilibrium, and the activation energy function monotonically decreases as the system relaxes, then at a specific time during the relaxation period, the timescale for escape events matches the relaxation timescale, and the system jumps at this point almost surely. If the jump brings the system back out of equilibrium, the relaxation stage can start over again, and then the scenario repeats itself indefinitely, leading to a rather regular behavior (i.e., coherent spiking) which cannot occur in the noise-free system. This coherent spiking (due to SISR) cannot occur in the noise-free system simply because the system parameters are chosen such that a Hopf or Saddle-node bifurcation leading to deterministic spiking activity cannot occur. What is crucial for SISR is the interplay between different timescales, and not the type of bifurcation leading to deterministic spiking.