Weak-noise-induced transitions with inhibition and modulation of neural oscillations

TL;DR

This paper investigates how weak noise influences neural oscillations in a bistable FitzHugh-Nagumo model, revealing noise-induced inhibition and modulation of spiking activity depending on initial conditions and system parameters.

Contribution

It provides a detailed analysis of noise-induced transitions in neural dynamics, highlighting conditions for inhibition and modulation of oscillations in a bistable neuron model.

Findings

Weak noise can inhibit neural spiking activity.

Increasing noise strength can first inhibit then enhance spiking activity.

The phenomenon depends on initial conditions and system parameters.

Abstract

We analyze the effect of weak-noise-induced transitions on the dynamics of the FitzHugh-Nagumo neuron model in a bistable state consisting of a stable fixed point and a stable unforced limit cycle. Bifurcation and slow-fast analysis give conditions on the parameter space for the establishment of this bi-stability. In the parametric zone of bi-stability, weak-noise amplitudes may strongly inhibit the neuron's spiking activity. Surprisingly, increasing the noise strength leads to a minimum in the spiking activity, after which the activity starts to increase monotonically with increase in noise strength. We investigate this inhibition and modulation of neural oscillations by weak-noise amplitudes by looking at the variation of the mean number of spikes per unit time with the noise intensity. We show that this phenomenon always occurs when the initial conditions lie in the basin of attraction of the stable limit cycle. For initial conditions in the basin of attraction of the stable fixed point, the phenomenon however disappears, unless the time-scale separation parameter of the model is bounded within some interval. We provide a theoretical explanation of this phenomenon in terms of the stochastic sensitivity functions of the attractors and their minimum Mahalanobis distances from the separatrix isolating the basins of attraction.

Equations

$$\begin{split}\left\{\begin{array}[]{lcl}dv_{\tau}&=&\frac{1}{\varepsilon}f(v_{\tau},w_{\tau})d\tau+\frac{\sigma}{\sqrt{\varepsilon}}dW_{\tau},\\ dw_{\tau}&=&g(v_{\tau},w_{\tau})d\tau,\end{array}\right.\end{split}$$

$$\begin{split}\left\{\begin{array}[]{lcl}dv_{t}&=&f(v_{t},w_{t})dt+\sigma dW_{t},\\ dw_{t}&=&\varepsilon g(v_{t},w_{t})dt,\end{array}\right.\end{split}$$

$$\begin{split}\left\{\begin{array}[]{lcl}f(v,w)&=&v(a-v)(v-1)-w,\\ g(v,w)&=&bv-cw,\end{array}\right.\end{split}$$

$$(v_{e},w_{e}):=\Big{\{}(v,w)\in\mathbb{R}^{2}:f(v,w)=g(v,w)=0\Big{\}}.$$

$$\displaystyle\begin{split}\left\{\begin{array}[]{lcl}\frac{b}{c}v&=&-v^{3}+(a+1)v^{2}-av,\\ w&=&\frac{b}{c}v,\end{array}\right.\end{split}$$

$$v_0=0\text{ and }{v}_{1,2}=\frac{a+1}{2}\pm \sqrt{\frac{(a-1{)}^2}{4}-\frac{b}{c}},$$

$$\frac{(a-1{)}^2}{4}\ge \frac{b}{c}.$$

$$v_0<a<v_{1,2}<1,$$

$$\frac{(a-1{)}^2}{4}=\frac{b}{c},$$

$$\frac{b}{c}=-a.$$

$$-v^3+v=\frac{b}{c}v,$$

$$\begin{split}\left\{\begin{array}[]{lcl}\varepsilon\frac{d\varphi}{d\tau}&=&-3v_{\star}^{2}\varphi+2(a+1)v_{\star}\varphi-a\varphi-\eta,\\ \\ \frac{d\eta}{d\tau}&=&b\varphi-c\eta.\end{array}\right.\end{split}$$

$$J_{ij}=\left(\begin{array}[]{cc}\frac{1}{\varepsilon}(-3v_{\star}^{2}+2(a+1)v_{\star}-a)&\>\>\>\>\>-\frac{1}{\varepsilon}\\ \\ b&\>\>\>\>\>-c\end{array}\right).$$

$$-3v_{\star }^2+2(a+1)v_{\star }-a<0.$$

$$v_1<-\frac{1}{\sqrt{3}},$$

$$v_1=-\frac{1}{\sqrt{3}},$$

$$\begin{split}\left\{\begin{array}[]{lcl}\frac{1}{\varepsilon}(-3v_{\star}^{2}+2(a+1)v_{\star}-a)+c-2\sqrt{\frac{b}{\varepsilon}}<0,\\ \frac{1}{\varepsilon}(-3v_{\star}^{2}+2(a+1)v_{\star}-a)+c+2\sqrt{\frac{b}{\varepsilon}}>0,\end{array}\right.\end{split}$$

$$\frac{1}{\epsilon }(-3v_{\star }^2+2(a+1)v_{\star }-a)-c=0.$$

$$c^2<\frac{b}{\epsilon },$$

$$v_{\star \star }=\frac{a+1}{3}\pm \sqrt{\frac{(a+1{)}^2}{9}-\frac{a+\epsilon c}{3}},$$

$$3\epsilon c\le a^2-a+1.$$

$$v_{\pm }=\frac{a+1}{3}\pm \sqrt{\frac{(a+1{)}^2}{9}-\frac{a}{3}},$$

$$a=-\epsilon c.$$

$$a=-\frac{b}{c},$$

$$2(a+1)\left(\frac{a+1}{2}\pm \sqrt{\frac{(a-1{)}^2}{4}-\frac{b}{c}}\right)-a$$

$$=3{\left(\frac{a+1}{2}\pm \sqrt{\frac{(a-1{)}^2}{4}-\frac{b}{c}}\right)}^2+\epsilon c,$$

$$-\frac{a^2}{2}+a-\frac{1}{2}+3\frac{b}{c}-\epsilon c=\pm (a+1)\sqrt{\frac{(a-1{)}^2}{4}-\frac{b}{c}},$$

$$\frac{1}{4}(a-1{)}^4+(3\frac{b}{c}-\epsilon c{)}^2-(a-1{)}^2(3\frac{b}{c}-\epsilon c)$$

$$=\frac{1}{4}(a-1{)}^2(a+1{)}^2-(a+1{)}^2\frac{b}{c},$$

$$-a^3+2a^2-a-2a^2\frac{b}{c}+8a\frac{b}{c}-2\frac{b}{c}+a^2\epsilon c$$

$$-2a\epsilon c+\epsilon c+9\frac{b^2}{c^2}-6b\epsilon +{\epsilon }^2{c}^2=0,$$

$$(a-1{)}^2(-a-2\frac{b}{c})+4a\frac{b}{c}+9\frac{b^2}{c^2}+{\epsilon }^2{c}^2$$

Full text

††journal: X - Journal

1 Introduction

Fixed points, periodic, quasi-periodic or chaotic orbits are typical solutions of deterministic nonlinear dynamical systems. Multi-stability is common, as a dynamical system typically possesses two or more mutually exclusive stable solutions (attractors). For a given set of the parameters, coexistent stable states represented by different (or identical) types of attractors in the phase space of the system are by topological necessity separated by some unstable states. In neurodynamics for example, spiking neurons may possess coexistent quiescent (fixed point) and tonic spiking states (limit cycle) Paydarfar , or distinct periodic and chaotic spiking states G. Cymbalyuk . A given state can be reached if the system starts from a set of initial conditions within the state’s basin of attraction. Otherwise, an external perturbation can be used to switch the system from one stable attractor to another. When noise is introduced into the system, random trajectories can visit different stable states of the system by jumping over the unstable ones.

Important and challenging problems in multi-stable systems are to find the residence times of random trajectories in each stable state and its statistics, and the critical value of the noise amplitude and control parameters at which noise-induced jumping becomes significant. The analytical treatment of such problems based on the Fokker-Planck Equation (FPE) becomes complicated for $n$-dimensional dynamical systems, $n\geq 2$, and therefore, various approximations were developed and are now commonly used N. G. Van Kampen ; M.I Freidlin .

The quasi-potential method gives exponential asymptotics for the stationary probability density. In the vicinity of the deterministic attractor, the first approximation of the quasi-potential is a quadratic form G. Mil'shtein , leading to a Gaussian approximation of the stationary probability density of the FPE. The corresponding covariance matrix characterizes the stochastic sensitivity of the deterministic attractor: its eigenvalues and eigenvectors define the geometry of bundles of stochastic trajectories around the deterministic attractors. The Gaussian distribution centered on an attractor can be viewed as a confidence ellipsoid, while a minimal distance from this ellipsoid to the boundary separating the basins of attraction is proportional to the escape probability I. Bashkirtseva . The appropriate measure for this distance is the so-called Mahalanobis distance P. C. Mahalanobis , the distance from a point to a distribution.

The residence time of random trajectories in a basin of attraction depends on two factors. The first factor is the geometry of the basin of attraction, e.g., the larger the distance is between an attractor and the separatrix isolating its basin of attraction, the longer is the residence time of phase trajectories in the basin. Second, the attractors are sensitive to random perturbations: the higher the stochastic sensitivity function (SSF) is, the higher is the probability to escape from the basin of attraction, and thus the shorter is the residence time M.I Freidlin . Therefore, considering only the geometrical arrangement of stable attractors and the separatrix (the Euclidean distance between them) might not be sufficient for a theoretical explanation of a stochastic phenomenon, like the one to be investigated in the present work, and so the sensitivity of the attractors to random perturbations must also be taken into account. The Mahalanobis distance, which combines geometrical and stochastic sensitivity aspects of the dynamics, allows for a proper theoretical explanation of the transitions between attractors.

The effects of noise in neurobiological dynamical systems have been intensively investigated, for both single neurons and neural networks. Some of the most studied noise-induced phenomena are: stochastic resonance (SR) Lindner ; Longtin ; Collins , coherence resonance (CR) Pikovsky , and noise-induced synchronization Kim . During SR, the neuron’s spiking activity becomes more closely correlated with a sub-threshold periodic input current in the presence of an optimal level of noise. In $1997$, Pikovsky and Kurths showed that CR is basically SR in the absence of a periodic input current. During CR, noise can activate the neuron by producing a sequence of pulses which can achieve a maximal degree of coherence for an optimal noise amplitude if the system is in the neighborhood of its Andronov-Hopf bifurcation. We notice in these phenomena that noise has a facilitatory effect and leads only to increased responses.

More recently, it was discovered both experimentally Paydarfar and theoretically Gutkin et al. Gutkin1 ; Gutkin2 (see also Jost ) that noise can also turn off repetitive neuronal activity . Gutkin1 ; Gutkin2 ; Jost used the Hodgkin-Huxley equations in bistable regime (fixed point and limit cycle) with a mean input current consisting of both a deterministic and random input component, to computationally confirm the inhibitory and modulation effects of Gaussian noise on the neuron’s spiking activity. They found that there is a tuning effect of noise that has the opposite character to SR and CR, which they termed inverse stochastic resonance (ISR). Very recently (August 2016), the first experimental confirmation of ISR and its plausible functions in local optimal information transfer was reported in Buchin , where the Purkinje cells that play a central role in the cerebellum are used for the experiment. During ISR, weak-noise amplitudes may strongly inhibit the spiking activity down to a minimum level (thereby decreasing the mean number of spikes to a minimum value), after which the activity starts and continuously increases with increasing noise amplitude (thereby monotonically increasing the mean number of spikes with increasing noise amplitude). In Gutkin1 ; Gutkin2 , it was shown that ISR occurred and persisted regardless of which basin of attraction the initial conditions are chosen from, provided the deterministic input current component is above its Andronov-Hopf bifurcation value.

In the present work, we investigate ISR in a theoretical neuron model in the absence of a deterministic input current component. We consider the model with only a random input component, which is in a state of bi-stability consisting of a stable fixed point and a stable unforced limit cycle. We show that ISR occurs as well in this case and greatly depends not only on the location of the initial conditions, but also on the time-scale separation parameter of the model. More precisely, we show that ISR always occurs when the initial conditions are chosen from the basin of attraction of the stable limit cycle. When the initial conditions are in the basin of attraction of the stable fixed point, we show that ISR disappears, except interestingly when the time-scale separation parameter of the model lies within a certain interval. A theoretical explanation of this phenomenon is given in terms of the SSFs of the stable attractors and their Mahalanobis distances from the separatrix.

This paper is organized as follows: In Sect.2, we present the theoretical neuron model used to analyze ISR. In Sect.3, we make explicit deterministic bifurcation and slow-fast analyses of the model equation and show how bi-stability consisting a stable fixed point and a stable unforced limit cycle establishes itself. In Sect.4, we make a stochastic sensitivity analysis of the stable attractors. In Sect.5, we investigate ISR through numerical simulations and provide a theoretical explanation of the phenomenon using the results in Sect.4. In Sect.6, we have concluding remarks.

2 Model equation and phenomenon

In this paper, we consider a stochastic perturbation of a version of the Fitzhugh-Nagumo (FHN) neuron model FitzHugh . We consider the resulting stochastic differential equation both in the slow time-scale $\tau$ (Eq.(1)) and in the fast time-scale $t$ (Eq.(2))

[TABLE]

[TABLE]

with the deterministic velocity vector fields given by

[TABLE]

where $(v,w)\in\mathbb{R}^{2}$ represent the activity of the action potential $v$ and the recovery current $w$ that restores the resting state of the model. We have as constant parameters $b>0$, $c>0$, and $a$ is often confined to the range $0<a<1$, but the case $a<0$ will be examined in this work.

We have a singular parameter, $0<\varepsilon:=\tau/t\ll 1$, i.e., the time-scale separation ratio between the slow time-scale $\tau$ and the fast time-scale $t$. We note that Eq.(1) preserves the sense of the dynamics on the trajectories of Eq.(2). In other words, the phase trajectories of both systems of dynamical equations have exactly the same dynamical behavior. The only difference is the speed of these trajectories in the phase space. Because the speeds of the trajectories do not affect in any way our analysis, we will work on both time-scales. The slow time-scale equation at some points allows for quicker conclusions in bifurcation analysis while the fast time-scale equation has an advantage in numerical simulations as it avoids the division by the very small parameter $\varepsilon$.

$dW_{t}$ is standard white noise, the formal derivative of Brownian motion with mean zero and unit variance, and $\sigma$ is the amplitude of this noise. The random term in Eq.(1) is rescaled in Eq.(2) according to the scaling law of Brownian motion. That is, if $W_{t}$ is a standard Brownian motion, then for every $\lambda\geq 0$, $\lambda^{-1/2}W_{\lambda t}$ is also a standard Brownian motion, i.e., the two processes have the same distribution Durrett .

Fig.1 shows the time series produced by the dynamics of the action potential variable $v$. In the deterministic case ($\sigma=0$), and for $a=-0.05$, $b=1.0$, $c=2.0$, and $\varepsilon=0.02785$, Eq.(2) can result in two different dynamics. In Fig.1a with initial conditions at \big{(}v(0),w(0)\big{)}=(0.001,0.001), the neuron has only sub-threshold oscillations with $v$ converging to zero and remaining at this value as the time $t$ increases. In Fig.1b with initial conditions now at \big{(}v(0),w(0)\big{)}=(-0.4,0.2), the neuron shows self-sustained supra-threshold oscillations. Thus, the system is bistable.

Fig.1c-e show a stochastic behavior ($\sigma>0$), with again $a=-0.05$, $b=1.0$, $c=2.0$, $\varepsilon=0.02785$, and the initial conditions all at \big{(}v(0),w(0)\big{)}=(-0.4,0.2). We count a spike when $v$ is greater than or equal to the threshold value $v_{th}=0.25$. In Fig.1c with a weak-noise amplitude, i.e., $\sigma=1.5\times 10^{-9}$, we have supra-threshold oscillations for a certain time length with $21$ spikes after which $v$ starts to converge to zero and the spiking eventually stops. In Fig.1d, when the noise amplitude is increased (but still relatively weak) to $\sigma=1.2\times 10^{-6}$, we have an even faster inhibition of the spiking with a smaller number of spikes. In this realization, we have only $3$ spikes. In Fig.1e, with a stronger noise amplitude, $\sigma=1.0\times 10^{-4}$, the number of spikes increases again up to $38$. Intuitively, it is surprising that weak-noise amplitudes inhibit the spiking activity of the neuron with the occurrence of a minimum in the number of spikes as the noise amplitude increases even though the initial conditions are exactly the same as in Fig.1b. We shall investigate in detail the mechanisms behind this phenomenon.

3 Bifurcation and slow-fast analysis

We now consider the deterministic dynamics corresponding to Eq.(1) when $\sigma=0$ and perform an explicit bifurcation and slow-fast analysis through which we find the parametric conditions for the establishment of a bi-stability regime consisting of a stable fixed point and a stable unforced limit cycle. We note again that our model has no deterministic input current. At the fixed points $(v_{e},w_{e})$ (rest states of the neuron), the variables $v(t)$ and $w(t)$ reach a stationary state, while the set of fixed points is defined by the intersection of nullclines as

[TABLE]

From Eq.(4), we obtain the fixed point equations

[TABLE]

which has the solutions for $v$ as

[TABLE]

where the solutions $v_{1,2}$ exist if

[TABLE]

We always assume $b,c>0$. When $0<a<1$, we then have

[TABLE]

but when $a<0$, this no longer holds.

A bifurcation occurs when some of these fixed points coincide. We have $v_{1}=v_{2}$ if and only if

[TABLE]

and we have $v_{0}=v_{1}$ if

[TABLE]

When $a=-1$ and $\frac{b}{c}=1$, all three fixed points coincide. The case $a=-1$ is somewhat simpler because the cubic polynomial in Eq.(5) reduces to

[TABLE]

and much of the subsequent analysis will only be carried out for that particular case. For the moment, however, we return to the general case. The linearization of Eq.(1) at such a fixed point $v_{\star}$ is

[TABLE]

In order to determine the bifurcation behavior at such a fixed point, we need to investigate the eigenvalues of the Jacobian matrix given by

[TABLE]

The stability of the fixed points $v_{\star}$ will depend on the signs of the trace and determinant of $J_{ij}$. For a fixed point $v_{\star}$ to be stable, we should have $\mathrm{tr}J_{ij}<0$ and $\mathrm{det}J_{ij}>0$. Since $\varepsilon,c>0$, we have $\mathrm{tr}J_{ij}<0$ and $\mathrm{det}J_{ij}>0$ only if

[TABLE]

Eq.(14) means that the fixed point $v_{\star}$ has to be on the part of the cubic polynomial $-v^{3}+(a+1)v^{2}-av$ that has a negative derivative (see the red curve in Fig.3b). When we have three distinct fixed points $v_{\star}$ in Eq.(6), this could hold for the leftmost and the rightmost of them, and these two would then be stable, while the middle one would be unstable. The fixed point $v_{0}=0$ is stable for $a>0$, and unstable for $a<0$. We should point out, however, that Eq.(14) is a sufficient, but not a necessary condition for the stability of a fixed point.

When we vary the parameters so that two of the fixed points merge, we obtain a saddle-node type bifurcation. In the case $a=-1$, the fixed points $v_{1,2}$ (which are symmetric in this case, that is, $v_{1}=-v_{2}$) are stable according to Eq.(14) if and only if

[TABLE]

that is, if $v_{1}$ is to the left of the local minimum of the cubic polynomial $w=-v^{3}+v$ (equivalently, $v_{2}>\frac{1}{\sqrt{3}}$ is to the right of the local maximum). The limiting case

[TABLE]

of Eq.(15) corresponds to $\frac{b}{c}=\frac{2}{3}$, by Eq.(6). In this case, $v_{1}$ is precisely the minimum of the cubic polynomial $-v^{3}+v$ in Eq.(11). Likewise, $v_{2}$ then is the maximum of that polynomial. This case will reoccur below in Eq.(29).

When Eq.(7) is not satisfied, $v_{0}=0$ is the only fixed point. For $v_{0}$, the determinant of Eq.(13) is $\frac{1}{\varepsilon}(ac+b)$, and the trace is $-\frac{a}{\varepsilon}-c$. Thus, $v_{0}$ is stable when $a$ is not too negative, but in the limit $\varepsilon\to 0$, stability only persists for $a\geq 0$. (Recall here that Eq.(14) was sufficient, but not necessary for the stability of a fixed point $v_{\star}$.)

We obtain complex conjugate eigenvalues if $0>\frac{1}{4}(\mathrm{tr}J_{ij})^{2}-\mathrm{det}J_{ij}$, and the real part of these eigenvalues vanishes when $\mathrm{tr}J_{ij}=0$. For Eq.(13), we thus have the condition for complex conjugate eigenvalues as

[TABLE]

and the condition for a vanishing real part

[TABLE]

In order to be able to satisfy Eq.(17) and Eq.(18) simultaneously (and therefore have a Andronov-Hopf bifurcation), the coefficients $b,c\ (>0)$ need to satisfy

[TABLE]

which is easily satisfied for small $\varepsilon>0$, since we assume $b>0$.

When we solve Eq.(18) for $v_{\star}$, we obtain

[TABLE]

which we can solve as long as

[TABLE]

We can then check when a solution $v_{\star\star}$ of Eq.(20) coincides with one of the points given by Eq.(6), in order to get a Andronov-Hopf bifurcation at one of those fixed points. We had already observed above that a fixed point on the decreasing part of the cubic $v$-nullcline curve $w=-v^{3}+(a+1)v^{2}-av$ (the red curve in Fig.3b) is stable as the eigenvalues of the linearization then have negative real parts. But this was a sufficient, but not a necessary condition. Therefore, when we vary the slope $\frac{b}{c}$ of the $w$-nullcline and move a fixed point to the middle increasing part of the cubic polynomial, that fixed point may eventually lose its stability. When the parameter regime just identified for a Andronov-Hopf bifurcation is right, that loss of stability will occur through a Andronov-Hopf bifurcation.

In order to see the significance of Eq.(21), we observe that the local extrema of the cubic polynomial $-v^{3}+(a+1)v^{2}-av$ are at

[TABLE]

and since $\varepsilon c>0$, $\frac{a+1}{3}-\sqrt{\frac{(a+1)^{2}}{9}-\frac{a+\varepsilon c}{3}}>\frac{a+1}{3}-\sqrt{\frac{(a+1)^{2}}{9}-\frac{a}{3}}$, that is, the left Andronov-Hopf bifurcation is to the right of the minimum of our cubic polynomial, hence on its ascending branch. Whenever a fixed point is on the left descending branch, that is, to the left of that minimum, it is stable, and stability persists a little into the ascending branch, but when $\varepsilon\to 0$, the Andronov-Hopf bifurcation point, that is, where the fixed point looses its stability, moves towards that minimum. For $\varepsilon>0$, the Andronov-Hopf bifurcation occurs to the right of the minimum.

We return to the investigation of Eq.(18). In particular, for $v_{\star}=0$, Eq.(18) does not have a solution under our constraint $c>0$ when also $a>0$. When $a<0$, Eq.(18) is satisfied for $v_{\star}=0$ when

[TABLE]

When $a>0$ and Eq.(21) holds, the solutions $v_{\star\star}$ of Eq.(20) are positive, and we might then tune the parameter $b$ which does not occur in Eq.(20) so that one of those $v_{\star\star}$ coincides with $v_{1}$ or $v_{2}$ from Eq.(6).

When we have equality in Eq.(19), i.e., $c^{2}=\frac{b}{\varepsilon}$, Eq.(23) becomes

[TABLE]

that is, Eq.(10). In that case, we have a limit of a Andronov-Hopf bifurcation co-occurring with a saddle-node bifurcation.

For $v_{\star}=v_{1,2}$ in Eq.(6), we obtain

[TABLE]

that is,

[TABLE]

which implies

[TABLE]

hence

[TABLE]

which is equivalent to

[TABLE]

Thus, Eq.(3) is the condition on the parameters for a Andronov-Hopf bifurcation at one of the equilibria $v_{1,2}$.

In particular, in the limit $\varepsilon=0$, we get the condition

[TABLE]

For the case $a=-1$, we directly get from Eq.(3)

[TABLE]

Thus, for $\varepsilon>0$, in Eq.(28), we have $\frac{b}{c}>\frac{2}{3}$, and therefore the fixed point is to the right of the minimum, that is, on the increasing part of the cubic curve $-v^{3}+v$, in accord with what we had said after Eq.(20).

In the limit $\varepsilon=0$, we obtain

[TABLE]

as the condition for a Andronov-Hopf bifurcation at an equilibrium $v_{1,2}$. This is the case of Eq.(16) (recalling Eq.(6)). Here, $v_{1}$ and $v_{2}$ lose their stability, and below in the slow-fast analysis, these are also the points where the critical manifold will not be normally hyperbolic, and where a switch from slow to fast dynamics will occur, generating an unforced limit cycle. Since we had identified Eq.(29) as the parameter constellation for a Andronov-Hopf bifurcation, this is a singular limiting situation for a Andronov-Hopf bifurcation.

Now, we use slow-fast techniques to understand how a stable unforced limit cycle emerges from Eq.(1) (with $\sigma=0$). For a more detailed introduction to multiple time-scale dynamics see Kuehn . In the singular limit $\varepsilon=0$, we define the critical manifold $\mathcal{M}_{0}$ of Eq.(1) which coincides with our cubic polynomial nullcline:

[TABLE]

$\mathcal{M}_{0}$ can be viewed as the algebraic constraint of the differential-algebraic slow subsystem Eq.(34) whose initial conditions must satisfy this constraint for solutions to exist. From Eq.(30), we have

[TABLE]

with

[TABLE]

Thus, $\mathcal{M}_{0}$ naturally splits into three parts: two decreasing stable branches and a middle increasing unstable branch, see the red curve in Fig.2 or Fig.3b. $\mathcal{M}_{0}$ looses its normal hyperbolicity at the two singular points $v_{\pm}=\frac{a+1}{3}\pm\frac{1}{3}\sqrt{a^{2}-a+1}$, where it changes its stability property. These two points are the local extrema of $\mathcal{M}_{0}$. In fact, at $v_{\pm}=\frac{a+1}{3}\pm\frac{1}{3}\sqrt{a^{2}-a+1}$, the existence and uniqueness theorems for ordinary differential equations (ODEs) do not longer apply, and because of this, the solutions of the slow subsystem in Eq.(34) are forced to leave $\mathcal{M}_{0}$ at these singular points.

In the singular limit $\varepsilon=0$, the dynamics of the fast variable $v$ on $\mathcal{M}_{0}$ (i.e., a 1-D dynamical system of the variable $v$ whose phase space is $\mathcal{M}_{0}$), is obtained from Eq.(30) by implicit differentiation

[TABLE]

and using the algebraic constraint on $w$ in Eq.(30), we eliminate this variable to get the slow flow for the variable $v$ as

[TABLE]

which, of course, becomes singular at the points $v_{\pm}$. In fact, at these points, since the critical manifold $\mathcal{M}_{0}$ loses its stability, the slow flow in Eq.(34) should detach from $\mathcal{M}_{0}$ and should become fast, that is, satisfy $\frac{dw}{d\tau}=0$ and horizontally jump to another stable branch of $\mathcal{M}_{0}$.

In Fig.2, all the black trajectories (with single and double arrows) represent the solution of the slow flow of Eq.(34) on $\mathcal{M}_{0}$. Because of the failure of the existence and uniqueness theorems of ODEs at $v_{-}$, the lower horizontal (fast) trajectory (in black with double arrow) leaves the left stable part of $\mathcal{M}_{0}$ at its local minimum at $v_{-}$ to the right stable part of $\mathcal{M}_{0}$. For the same reason, the upper horizontal (fast) trajectory (in black with double arrow) leaves the right stable part of $\mathcal{M}_{0}$ at its local maximum at $v_{+}$ to the left stable part of $\mathcal{M}_{0}$. In this same figure, the non-horizontal (slow) trajectories (all in black with a single arrow) of Eq.(34), evolve on $\mathcal{M}_{0}$ towards the singular points at $v_{\pm}$, where they eventually leave $\mathcal{M}_{0}$.

The solution (or more precisely, the singular solution with $\varepsilon=0$) of the slow subsystem in Eq.(34) is related to the solution of the full system Eq.(1) with $\varepsilon>0$ by Fenichel’s theorem Feni1 ; Feni . By that theorem, for $0<\varepsilon\ll 1$, the slow manifold $\mathcal{M}_{\varepsilon}$ is a perturbation of the critical manifold $\mathcal{M}_{0}$, and it has the following properties:

• •

(F1) $\mathcal{M}_{\varepsilon}$ is diffeomorphic to $\mathcal{M}_{0}$.

• •

(F2) $\mathcal{M}_{\varepsilon}$ has distance $\mathcal{O}(\varepsilon)$ from $\mathcal{M}_{0}$.

• •

(F3) The flow on $\mathcal{M}_{\varepsilon}$ converges to the slow flow on $\mathcal{M}_{0}$ as $\varepsilon\rightarrow 0$.

• •

(F4) $\mathcal{M}_{\varepsilon}$ is $C^{r}$-smooth for any $r<\infty$ (as long as $f,g\in C^{\infty}$).

• •

(F5) $\mathcal{M}_{\varepsilon}$ is normally hyperbolic, with the same stability properties w.r.t. the fast variable $v$ as $\mathcal{M}_{0}$ (attracting, repelling or saddle-type).

• •

(F6) $\mathcal{M}_{\varepsilon}$ is usually not unique. Manifolds satisfying (F1)-(F5) lie at distance $\mathcal{O}(e^{-}{K/\varepsilon})$ from each other $K>0$, $K=\mathcal{O}(1)$.

• •

(F7) Similar conclusions hold for the stable/unstable manifolds of $\mathcal{M}_{0}$.

Most importantly, the flow of Eq.(1) (with $0<\varepsilon\ll 1$) on $\mathcal{M}_{\varepsilon}$ will follow the slow flow of Eq.(34) (with $\varepsilon=0$) on $\mathcal{M}_{0}$. By the definitions in Eq.(4) and Eq.(30), we see that the fixed points $v_{\star}$ in Eq.(6) also lie on $\mathcal{M}_{0}$. We have the fixed points $v_{0}=0$ and $v_{2}$ on the decreasing part of $\mathcal{M}_{0}$, and they are therefore stable, while the fixed point $v_{1}$ is located between $v_{0}=0$ and $v_{2}$ and unstable.

In Fig.2, the blue trajectories with arrows pointing in the direction of the flow on the slow manifold $\mathcal{M}_{\varepsilon}$ (not shown, but at a distance $\mathcal{O}(\varepsilon)$ from $\mathcal{M}_{0}$), represent solutions of Eq.(1) with $\varepsilon=0.1$. They converge towards to the stable fixed points $v_{0}$ and $v_{2}$ located respectively on the left and right stable decreasing parts of the slow manifold $\mathcal{M}_{\varepsilon}$.

For a better visualization, we can henceforth discuss the dynamics with respect to $\mathcal{M}_{0}$, since Fenichel’s theorem tells us that the same dynamical behavior takes place on $\mathcal{M}_{\varepsilon}$. With the fixed point configuration: $v_{0}<v_{1}<v_{2}$ with $v_{0}$ and $v_{2}$ stable and located respectively on the left and right decreasing parts of $\mathcal{M}_{0}$, $v_{1}$ unstable and located on the increasing part of $\mathcal{M}_{0}$, trajectories of Eq.(1) cannot exhibit a spiking behavior (i.e., a limit cycle solution cannot emerge). This is because trajectories converge and stay at the stable fixed points $v_{0}$ and $v_{2}$ whenever they encounter them on decreasing parts of $\mathcal{M}_{0}$. Hence, these trajectories cannot evolve and reach the singular points $v_{\pm}$ located at the extrema of $\mathcal{M}_{0}$, where they can jump from the left to the right and then from the right to the left part of $\mathcal{M}_{0}$ to produce a limit cycle solution. In the blue trajectories in Fig.2, they stick at the fixed points $v_{0}=0$ and $v_{2}$.

Thus, it depends on the relative positions of the singular points at $v_{\pm}$ and the fixed points at $v_{0}$ and $v_{2}$ on the critical manifold $\mathcal{M}_{0}$ whether the flow of Eq.(1) will first reach a stable fixed point (and stay there with no possibility for a limit cycle) or first reach a singular point $v_{\pm}$ and detach from $\mathcal{M}_{0}$ with the emergence of a limit cycle.

When, say, $v_{0}$ is the smallest fixed point and is located to the right of $v_{-}$ and $v_{2}$, the largest fixed point, to the left of $v_{+}$, we expect a periodic cycle for $0<\varepsilon\ll 1$. For instance, when we start close to the left stable branch of $\mathcal{M}_{0}$, by Fenichel’s theorem, the flow closely and slowly follows $\mathcal{M}_{0}$ until we get into the vicinity of $v_{-}$. There, $\mathcal{M}_{0}$ becomes unstable, and the flow will move away from it and become fast and therefore move not perfectly horizontally this time, but with some inclination (because $\varepsilon>0$) to the right until it comes into the vicinity of the right stable branch of $\mathcal{M}_{0}$. It will then again move slowly and close to $\mathcal{M}_{0}$ until it gets into the vicinity of $v_{+}$. Hence it moves away fast, also not perfectly horizontally but with some inclination to the left, until it encounters the left stable branch again, and the cycle repeats.

We therefore see that even in the absence of a deterministic input current, the slow-fast structure of Eq.(1) (with $\sigma=0$) can naturally induce a limit cycle solution (spiking) if the flow in Eq.(1) leaves from the neighborhood of the stable parts of $\mathcal{M}_{0}$ at the singular points $v_{\pm}$ before it encounters a stable fixed point. This can only happen if the fixed point is located to right of $v_{-}$ or to the left of $v_{+}$.

Comparing Eq.(14) and Eq.(31), we see that the fixed point $v_{\star}$ is stable for all $\varepsilon>0$ if it is to the left of the minimum of $\mathcal{M}_{0}$. In that case, it is on the left stable branch of $\mathcal{M}_{0}$, and the dynamics on that stable branch will therefore converge towards $v_{\star}$, without a limit cycle emerging. Analogously, a fixed point $v_{\star}$ is stable if it is on the right stable branch of $\mathcal{M}_{0}$ with no possibility for a limit cycle as well. Again, however, these are sufficient, but not necessary conditions for stability of the fixed points. For $\varepsilon>0$, stability persists a little into the middle (increasing) part of $\mathcal{M}_{0}$ as we have found when we discussed the possibility of an Andronov-Hopf bifurcation. Therefore, to have a bi-stability regime consisting of a stable fixed point and a stable unforced limit cycle, the fixed point should be located on the middle unstable branch of $\mathcal{M}_{0}$ and has be to stable. This can be done by choosing the parameters such that the fixed point is located between the minimum $v_{-}$ of $\mathcal{M}_{0}$ and its Andronov-Hopf bifurcation value.

For the purpose of this work, we shall henceforth consider the situation where $v_{0}=0$ is the only fixed point, that is we choose $a$, $b$, and $c$ such that Eq.(7) is not satisfied. The persistence of stability on the middle unstable part of $\mathcal{M}_{0}$ also occurs for $v_{0}=0$. For $a<0$, $v_{0}=0$ is to the right of the minimum $v_{-}=\frac{a+1}{3}-\sqrt{\frac{(a+1)^{2}}{9}-\frac{a}{3}}$ of $\mathcal{M}_{0}$ (i.e., $v_{-}<v_{0}$) and therefore in the region where it eventually loses its stability through an Andronov-Hopf bifurcation when $\varepsilon$ increases. It now suffices to choose specific values of $a$ ($a<0$) and $c$ (both values of $a$ and $c$ also not satisfying Eq.(7)) such that $v_{0}<-\frac{a}{c}=\varepsilon$ (see Eq.(23)) to have a stable fixed point $v_{0}$ to the right of $v_{-}$.

With the stable fixed point $v_{0}=0$ on the middle part of $\mathcal{M}_{0}$, and from the slow-fast analysis above, we also have a stable limit cycle surrounding this fixed point. For topological reasons these attractors should be separated from each other by a repeller, in this case an unstable limit cycle. In fact, the unstable limit cycle is the boundary of the basin of attraction of the stable fixed point $v_{0}$. This immediately indicates that the Andronov-Hopf bifurcation of the fixed point $v_{0}$ which eventually occurs as $\varepsilon$ increases should be sub-critical.

We choose and maintain throughout this work the values of the parameters as: $a=-0.05$, $b=1.0$, and $c=2.0$. For these values, we have $v_{-}=-0.25305<v_{0}=0<-\frac{a}{c}=0.025$ and therefore $v_{0}$ is located on the middle part of $\mathcal{M}_{0}$ and it is stable. The Andronov-Hopf bifurcation value of $\varepsilon$ is computed from Eq.(23) as $\varepsilon_{hp}=0.025$.

For these values of the system parameters, see Fig.3a, we computed the bifurcation diagram by selecting the maximum values of the action potentials $v$ as a function of the bifurcation parameter $\varepsilon$. For $0.024\leq\varepsilon<0.025$, the fixed point $v_{0}=0$ is unstable as $\mathrm{det}J_{ij}=\frac{0.9}{\varepsilon}>0$ and $\mathrm{tr}J_{ij}=\frac{0.05}{\varepsilon}-2>0$ and surrounded by a stable limit cycle, and therefore no bi-stability. At $\varepsilon=\varepsilon_{hp}=0.025$, a sub-critical Andronov-Hopf bifurcation occurs and the unstable fixed point $v_{0}=0$ changes its stability. For $0.025<\varepsilon<0.027865$, the fixed point $v_{0}=0$ is stable ($\mathrm{det}J_{ij}>0$ and $\mathrm{tr}J_{ij}<0$ for those values of $\varepsilon$) and co-exist with the stable limit cycle. Thus, for $0.025<\varepsilon<0.027865$ we have a bi-stability regime. At $\varepsilon=\varepsilon_{sn}=0.027865$, we have a saddle-node bifurcation of limit cycles, in which case the stable limit cycle surrounding the stable fixed point $v_{0}=0$ shrinks and eventually collides with the boundary of the basin of attraction of $v_{0}$ (i.e., the unstable limit cycle). In this saddle-node bifurcation, the stable and the unstable limit cycle annihilate each other leaving behind the stable fixed point $v_{0}$. The fixed point $v_{0}$ maintains its stability in the interval $0.027865<\varepsilon\leq 0.029$, within which we have no bi-stability as there is only one attractor in the entire phase space.

In the bi-stability regime $\varepsilon_{hp}<\varepsilon<\varepsilon_{sn}$, depending on whether the initial conditions are chosen in the basin of attraction of the fixed point or in that of the limit cycle, the dynamics will converge to either the fixed point or to the limit cycle. This behavior is seen in Fig.3b (also already seen in the time-series in Fig.1a and b) which shows a phase portrait of one trajectory with initial conditions in the basin of attraction of the stable fixed point at $(v_{0},w_{0})=(0,0)$ (the blue dot at the origin) and two other trajectories with initial conditions in the basin of attraction of the stable limit cycle (the blue closed curve). In this work, we will focus on weak-noise effects on the spiking dynamics of Eq.(1) with $\varepsilon_{hp}<\varepsilon<\varepsilon_{sn}$.

4 Stochastic sensitivity analysis and the Mahalanobis metric

We now introduce noise to the neuron model (i.e., $\sigma>0$) and perform a stochastic sensitivity analysis of the stable attractors. In this section, we analyze the neuron model on the fast time-scale $t$. The stochastic sensitivity matrix associated to a stochastic dynamical system is an asymptotic characteristic of the random attractors of the system Bashkirtseva1 . For our model equation Eq.(2) with $0<\sigma\ll 1$, it allows us to approximate a spread of random trajectories around the stable fixed point $(v_{0},w_{0})=(0,0)$ and stable limit cycle which we now denote by \big{[}\bar{v}(t),\bar{w}(t)\big{]}. The random trajectories in the basin of attraction of the stable fixed point $(v_{0},w_{0})$ evolve according the evolution of the probability density of the FPE corresponding to Eq.(2) Risken .

Suppose that a stationary solution, P\big{[}v(t),w(t)\big{]}, of this FPE exists. Generally, for $n$-dimensional systems with $n\geq 2$, one usually cannot find such a stationary probability density analytically Risken . This is the situation with Eq.(2). When $0<\sigma\ll 1$, the constructive asymptotics and approximations based on a quasi-potential function, $\varphi$, given in Eq.(35) are frequently used M.I Freidlin .

[TABLE]

A quadratic form of the quasi-potential gives a Gaussian approximation of P_{g}\big{[}v(t),w(t)\big{]} in the vicinity of the fixed point $(v_{0},w_{0})$,

[TABLE]

where $Z$ is the normalization constant and $\Omega_{ij}$ is the covariance matrix of random trajectories around the stable fixed point $(v_{0},w_{0})$, i.e., $\Omega_{ij}$ plays the role of the stochastic sensitivity matrix of this stable fixed point and it is determined by the algebraic equation

[TABLE]

J_{ij}=\left(\begin{array}[]{cc}-a&-1\\ \varepsilon b&-\varepsilon c\end{array}\right) is the Jacobian matrix at $(v_{0},w_{0})=(0,0)$ of the deterministic neuron equation Eq.(2). G_{ij}=\left(\begin{array}[]{cc}1&0\\ 0&0\end{array}\right) is the diffusion matrix of system Eq.(2) and $\top$ denotes the transpose.

As the fixed point $(v_{0},w_{0})$ is exponentially stable (that is, all the eigenvalues of $J_{ij}$ at $(v_{0},w_{0})=(0,0)$ have strictly negative real parts), the matrix equation in Eq.(41) has as its unique solution the stochastic sensitivity matrix $\Omega_{ij}\overline{}$ of the fixed point Bashkirtseva1 . The eigenvalues $\lambda_{k}(\varepsilon)$, $k=\{1,2\}$ of the stochastic sensitivity matrix $\Omega_{ij}$ define the variance of the random trajectories around the fixed point $(v_{0},w_{0})$. The largest eigenvalue (the largest SSF) $\lambda_{max}=\text{max}\{\lambda_{k}(\varepsilon)\}$ for each value of the singular parameter $\varepsilon$ indicates the sensitivity of the stable fixed point $(v_{0},w_{0})$ to the random perturbation. As $\lambda_{max}$ increases, the sensitivity of the $(v_{0},w_{0})$ to noise also increases. This means we have a higher probability of escape (i.e., shorter residence time) from the basin of attraction of the stable fixed point $(v_{0},w_{0})$ which we now denote for short as $\mathcal{B}(v_{0},w_{0})$.

The matrix equation Eq.(41) for Eq.(2) reduces to the system of algebraic equations

[TABLE]

In the parametric zone of bi-stability: $a=-0.05$, $b=1.0$, $c=2.0$, $\varepsilon_{hp}<\varepsilon<\varepsilon_{sn}$, the stochastic sensitivity matrix $\Omega_{ij}$ of the stable fixed point at $(v_{0},w_{0})=(0,0)$ of Eq.(2) is given by

[TABLE]

with the eigenvalues (the SSFs) given by

[TABLE]

where $\lambda_{max}=\lambda_{2}(\varepsilon)$ . The corresponding generalized eigenvectors are

[TABLE]

For a fixed noise strength $\sigma$, the difference between $\lambda_{1}$ and $\lambda_{2}$ reflects a spatial non-uniformity of the dispersion of the random trajectories around the fixed point $(v_{0},w_{0})$ in the direction of the eigenvectors $U_{1}$ and $U_{2}$ respectively. The dependence of $\lambda_{1}$ and $\lambda_{2}$ on the singular parameter $\varepsilon$ is shown in Fig.4.

Firstly, we observe that the SSFs diverge as we approach the Andronov-Hopf bifurcation value at $\varepsilon=\varepsilon_{hp}=0.025$. This means that the fixed point $(v_{0},w_{0})$ becomes more and more sensitive to noise as we approach the Andronov-Hopf bifurcation value and therefore the highest probability of escape (i.e., shortest residence time) from $\mathcal{B}(v_{0},w_{0})$ when $\varepsilon\approx\varepsilon_{hp}$.

Secondly, we observe that $\lambda_{2}$ diverges faster than $\lambda_{1}$ as $\varepsilon\rightarrow\varepsilon_{hp}=0.025$. This shows that the eigenvector $U_{2}$ localizes the main direction for deviations of random trajectories from the fixed point $(v_{0},w_{0})$, providing the direction in which the intersection with the unstable limit cycle at $[v(t),w(t)]$ is most probable.

The Mahalanobis metric is a widely used metric in cluster and discriminant analyses G. McLachlan . Basically, it measures the distance between a point $x$ and a distribution. This metric is a natural tool for the quantitative analysis of noise-induced transitions as it combines both the geometric distance from a random attractor to a point and the stochastic sensitivity of this attractor. The metric allows us to estimate a preference of the stable fixed point $(v_{0},w_{0})$ or the stable limit cycle \big{[}\bar{v}(t),\bar{w}(t)\big{]} in the stochastic dynamics of Eq.(2), when the random trajectory passes from one attractor to another. For $0<\sigma\ll 1$, the Mahalanobis distance from the unstable limit cycle \big{[}v(t),w(t)\big{]} (separatrix) to the stable fixed point $(v_{0},w_{0})$ or to the stable limit cycle \big{[}\bar{v}(t),\bar{w}(t)\big{]} is related to the residence time of trajectories in the corresponding basin of attraction: the larger the Mahalanobis distance, the longer is the residence time (i.e., lower probability of escape) in the corresponding basin.

In the stochastic sensitivity analysis of our fixed point $(v_{0},w_{0})$, we approximate the probability density by a Gaussian distribution in Eq.(4) centered at the stable fixed point at $(v_{0},w_{0})=(0,0)$. The Mahalanobis distance D_{m}\Big{\{}\big{[}v(t),w(t)\big{]};(v_{0},w_{0})\Big{\}} from a point

\big{[}v(t),w(t)\big{]} (i.e., a point on the unstable limit cycle) to the distribution of random trajectories around the stable attractor at $(v_{0},w_{0})$ is given by

[TABLE]

where $\Omega_{ij}$ is the stochastic sensitivity matrix of the fixed point $(v_{0},w_{0})$, and so the Gaussian approximation in Eq.(4) can be written in terms of the Mahalanobis distance as

[TABLE]

For Eq.(2), we calculate the Mahalanobis distances from the stable fixed point at $(v_{0},w_{0})=(0,0)$ to points on the unstable limit cycle at \big{[}v(t),w(t)\big{]}, and then we choose the minimal distance. We calculate coordinates of the unstable limit cycle numerically. This is done by assigning \big{[}v(t),w(t)\big{]} to the limiting values of the initial conditions \big{(}v(0),w(0)\big{)} such that infinitesimal perturbations (to the right and to the left) of these initial conditions will lead to the convergence of the trajectories either to the stable fixed point at $(v_{0},w_{0})$ or to the stable limit cycle at \big{[}\bar{v}(t),\bar{w}(t)\big{]} depending on which side the infinitesimal perturbation is made.

We have

[TABLE]

and using Eq.(4), we calculate the Mahalanobis distance from the stable fixed $(v_{0},w_{0})=(0,0)$ to the unstable limit cycle \big{[}v(t),w(t)\big{]} as

[TABLE]

Because of the dominance of $\lambda_{2}$ over $\lambda_{1}$, we numerically calculate the minimum Mahalanobis distance $D_{m}$ from the fixed point at $(v_{0},w_{0})=(0,0)$ to all points on the unstable limit cycle at \big{[}v(t),w(t)\big{]} by approximating the Mahalanobis distance in Eq.(4) along the eigenvector $U_{2}$, i.e.,

[TABLE]

Fig.5 shows the variation of the minimal Mahalanobis distance $D_{m}$ from the fixed point $(v_{0},w_{0})$ to the unstable limit cycle with the singular parameter $\varepsilon$. The Mahalanobis distance vanishes at the sub-critical Andronov-Hopf bifurcation value $\varepsilon_{hp}=0.025$, and increases with increasing $\varepsilon$. This means as $\varepsilon$ increases from $\varepsilon_{hp}$, the basin of attraction of the fixed point increases in the direction of the eigenvector $U_{2}$, and therefore a lower and lower probability of escape (i.e., longer residence time) from $\mathcal{B}(v_{0},w_{0})$ results.

From Figs. 4b and 5, we see that the two factors (stochastic sensitivity of an attractor and distance of the attractor to the separatrix) determining the length of the residence time of random trajectories in $\mathcal{B}(v_{0},w_{0})$ are not competing. Approaching $\varepsilon_{hp}$ from above increases the SSF of the fixed point and at the same time, decreases its Mahalanobis distance to the unstable limit cycle. This has the combined effect of considerably reducing the residence time of trajectories in $\mathcal{B}(v_{0},w_{0})$. In other words, there is a higher probability (lower probability) that the random trajectories escape from $\mathcal{B}(v_{0},w_{0})$ when $\varepsilon\rightarrow\varepsilon_{hp}$ ($\varepsilon\rightarrow\varepsilon_{sn}$).

We now apply a similar analysis to our randomly perturbed stable limit cycle. In the deterministic system ($\sigma=0$) of Eq.(2), with $\varepsilon_{hp}<\varepsilon<\varepsilon_{sn}$, we have an exponentially stable limit cycle defined by a $T$-periodic solution, \big{[}\bar{v}(t),\bar{w}(t)\big{]}=\big{[}\bar{v}(t+T),\bar{w}(t+T)\big{]}. For the transversal hyperplane $\varSigma_{t}$ in the neighborhood of any point \big{[}\bar{v}(t),\bar{w}(t)\big{]} on the stable limit cycle, the Gaussian approximation of the probability density reads

[TABLE]

Here, the stochastic sensitivity matrix is periodic in time, $\Theta_{ij}(t)=\Theta_{ij}(t+T)$. For an exponentially stable limit cycle, the largest Lyapunov exponent is [math] and the others are negative. Consequently, the matrix $\Theta_{ij}(t)$ is the unique solution of the Lyapunov equation Bashkirtseva1 ,

[TABLE]

with the conditions

[TABLE]

where $P_{ij}(t)$ is a matrix of the orthogonal projection onto the Poincaré section $\varSigma_{t}$ at the point $[\bar{v}(t),\bar{w}(t)]$ on the stable limit cycle, which is symmetric for our model equation Eq.(2), and whose entries are given by

[TABLE]

$J_{ij}(t)$ is the Jacobian of the deterministic neuron at a point $[\bar{v}(t),\bar{w}(t)]$ on the stable limit cycle and given by

[TABLE]

and the constant diffusion matrix $G_{ij}$ is the same as before.

The Mahalanobis distance from a point \big{[}v(t),w(t)\big{]} on the unstable limit cycle to the distribution of random trajectories around the stable limit cycle at \big{[}\bar{v}(t),\bar{w}(t)\big{]}, D_{m}\Big{\{}\big{[}v(t),w(t)\big{]};\big{[}\bar{v}(t),\bar{w}(t)\big{]}\Big{\}}, is also a periodic function of time and given by

[TABLE]

Because $\Theta_{ij}(t)$ is singular for Eq.(2), “+” means a pseudo-inverse in this case.

For 2-D systems, $\Theta_{ij}(t)$ can also be written in the form Bashkirtseva2

[TABLE]

and the Mahalanobis distance is given by

[TABLE]

Here, $\mu(t)=\mu(t+T)>0$ is the unique solution of the boundary problem

[TABLE]

with $T$-periodic coefficients

[TABLE]

$q(t)$ is a normalized vector orthogonal to the velocity vector field \left(\begin{array}[]{c}f\big{[}\bar{v}(t),\bar{w}(t)\big{]}\\ g\big{[}\bar{v}(t),\bar{w}(t)\big{]}\end{array}\right) and for Eq.(2) is given by

[TABLE]

We note that because our model is a $2D$ system, the hyperplane given by $\varSigma_{t}$ is a tangent line to the stable limit cycle solution which is normal to $q(t)$ at $[\bar{v}(t),\bar{w}(t)]$. The functions $\alpha(t)$ and $\beta(t)$ for Eq.(2) are given by

[TABLE]

[TABLE]

The explicit solution of Eq.(71) is given by

[TABLE]

Because $\mu(t)$ is $T$-periodic, we write

[TABLE]

and use the periodic property: $\alpha(t)=\alpha(t+T)$ with $\int_{t}^{t+T}\alpha(s)ds=\int_{0}^{T}\alpha(t+s)ds$, for a fixed $t$, to get the constant $C$ as

[TABLE]

With Eq.(70) and the numerical value of $\mu(t)$ in Eq.(79), the Mahalanobis distance is computed as in Eq.(4) and the minimal Mahalanobis distance is calculated by taking the minimum value of Eq.(4) over $t\in[0,T)$ and $(v,w)\in[v(t),w(t)]$. See Fig.6a.

[TABLE]

This set, we obtain the entries of $\Theta_{ij}(t)$ in Eq.(69) using the numerical value of $\mu(t)$ in Eq.(79). As in the case of the stable fixed point $(v_{0},w_{0})$, the eigenvalues $\lambda_{k}(t)$, $k=\{1,2\}$, of $\Theta_{ij}(t)$ characterize the distribution of random trajectories in the Poincaré section $\varSigma_{t}$ near a point \big{[}\bar{v}(t),\bar{w}(t)\big{]} of the stable limit cycle. The maximum of the largest eigenvalue indicates the SSF of the stable limit cycle. See Fig.6b.

5 Simulation results and discussion

In this section, numerical simulations are carried out with our model equation in the bistable regime to understand the ISR we observed in Fig.1. We want to see how ISR depends on which basin of attraction the initial conditions are located in, and how the singular parameter $\varepsilon$ affects ISR. We provide a theoretical explanation of the numerical results in terms of the results obtained in the stochastic sensitivity analysis of our model equation. We recall that the differences in the SSFs and Mahalanobis distances of our stable attractors define the direction of noise-induced transitions between them.

Using the fourth-order Runge-Kutta algorithm for stochastic processes Kasdin , simulations are carried out for $200$ realizations of the noise and for $7500$ unit time intervals for each realization, a sufficiently long time interval for convergence of solutions for $0<\sigma\ll 1$. In Fig.7, we depict the variations of the mean number of spikes $\langle N\rangle$ with the noise amplitude $\sigma$. The set of numerical results are for different values of Mahalanobis distances and SSFs of the stable attractors (encoded in the value of the singular parameter $\varepsilon$ as in Fig.6).

Sub-threshold responses are not counted as spikes. Again, we count a spike (supra-threshold response) when the action potential variable $v$ is greater than or equal to the threshold value of $v_{th}=0.25$. We show simulation results for six values of the singular parameter $\varepsilon\in(\varepsilon_{hp},\varepsilon_{sn})$, namely: $\varepsilon=0.02501$ which is in the vicinity of the Andronov-Hopf bifurcation value $\varepsilon_{hp}$ of the fixed point, $\varepsilon=0.02559$, $\varepsilon=0.0260$, $\varepsilon=0.0266$ at which both attractors have equal Mahalanobis distances, $\varepsilon=0.027673$ at which both attractors have equal SSFs, and $\varepsilon=0.02785$ which is in the vicinity of the saddle-node bifurcation value $\varepsilon_{sn}$ of the limit cycles.

The initial conditions \big{(}v(0),w(0)\big{)} are fixed in every simulation. In Fig.7, the red curves correspond to when \big{(}v(0),w(0)\big{)}\in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}. In this case, when $\sigma=0$, there are $106$ spikes. The inhibitory effect of the spiking activity begins as soon as $\sigma>0$, where we see the mean number of spikes $\langle N\rangle$ decreasing to a minimum value before increasing monotonically as $\sigma$ increases. We have ISR always occurring when \big{(}v(0),w(0)\big{)}\in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}.

The blue curves correspond to the situation where \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}). In this case, when $\sigma=0$, we have no spike, $\langle N\rangle=0$, and as soon as $\sigma>0$, $\langle N\rangle$ only increases monotonically and ISR does not occur. However, in Fig.7a and b (blue curves), interestingly, ISR does actually occur although \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}). We now explain these behaviors theoretically in terms of the Mahalanobis distances and the SSFs of the attractors.

In Fig.7a (blue curve), $\varepsilon=0.02501\gtrapprox\varepsilon_{hp}$, in which case the Mahalanobis distance of the fixed point, $D_{m}(fp)$, is small and smaller than the Mahalanobis distance of the limit cycle, $D_{m}(lc)$. At this same value of $\varepsilon$, the SSF of the fixed point, $\lambda_{max}(fp)$, is high and far higher than the SSF of the limit cycle, $\lambda_{max}(lc)$. Therefore, with \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}), very weak-noise amplitudes are already capable of kicking the random trajectories out of $\mathcal{B}(v_{0},w_{0})$ into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}, thereby increasing $\langle N\rangle$. For $0\leq\sigma\leq 0.25\times 10^{-5}$, $\langle N\rangle$ increases from [math] to a maximum of $79.8$. This same interval of $\sigma$ is incapable of causing the reverse event, i.e., not strong enough to kick random trajectories back into $\mathcal{B}(v_{0},w_{0})$ because of a large $D_{m}(lc)$ and a low $\lambda_{max}(lc)$ and therefore, no inhibitory effect of the spiking activity. As a result, $\langle N\rangle$ can only increase monotonically for $0\leq\sigma\leq 0.25\times 10^{-5}$. As soon as $\sigma>0.25\times 10^{-5}$, it becomes strong enough to kick the random trajectories it previously kicked into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} back to $\mathcal{B}(v_{0},w_{0})$, thereby decreasing $\langle N\rangle$ (inhibiting the spiking activity) down to a minimum value of $\langle N\rangle=69.3$ at $\sigma=1.0\times 10^{-5}$, before increasing monotonically with $\sigma$.

From a series of simulations carried out for several different values of $\varepsilon\in(\varepsilon_{hp},\varepsilon_{cr}]$ (not all shown except for $\varepsilon=0.02501$, $\varepsilon=0.02559$, and $\varepsilon_{cr}=0.0260$ in Fig.7a-c respectively), ISR remarkably persisted with \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}), except at the critical value $\varepsilon_{cr}$, where it just disappeared. That is, as $\varepsilon$ increased in the interval $(\varepsilon_{hp},\varepsilon_{cr})$ (with increasing $D_{m}(fp)$ and residence time in $\mathcal{B}(v_{0},w_{0})$), ISR became less and less pronounced and eventually disappeared when $\varepsilon\geq\varepsilon_{cr}$. This is so because as $D_{m}(fp)$ becomes larger and larger with increasing $\varepsilon\geq\varepsilon_{cr}$, trajectories stay longer and longer in $\mathcal{B}(v_{0},w_{0})$ (and therefore there is no spike), and as $\sigma$ increases, it becomes at each time just strong enough to kick the trajectories into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} and hence $\langle N\rangle$ can only increase monotonically from $\langle N\rangle=0$ with $\sigma$. See the blue curves in Fig.7c-f where $\varepsilon\geq\varepsilon_{cr}$, ISR does not occur as opposed to the cases in Fig.7a and b (blue curves) where $\varepsilon\in(\varepsilon_{hp},\varepsilon_{cr})$.

Still in Fig.7a (now the red curve), with \big{(}v(0),w(0)\big{)}\\ \in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}, in the thin interval of very weak-noise amplitudes $0\leq\sigma<0.2\times 10^{-5}$, $\langle N\rangle$ is almost constant, near $106$ (i.e., no considerable drop in $\langle N\rangle$ for this interval of $\sigma$). This “almost constant” value of $\langle N\rangle$ in that interval of $\sigma$ happens because at $\varepsilon=0.02501$, $D_{m}(lc)$ is large and $\lambda_{max}(lc)$ is very low and as \big{(}v(0),w(0)\big{)}\in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}, trajectories have the tendency of staying in this basin of attraction for a very long time and therefore almost no inhibition of the spiking activity occurs for $0\leq\sigma<0.2\times 10^{-5}$. As soon as $\sigma>0.2\times 10^{-5}$, it becomes strong enough to kick trajectories out of the relatively larger \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}. The inhibitory effect of noise then becomes pronounced with a clear decrease in $\langle N\rangle$ from about $106$ to a minimum of $78.3$ at $\sigma=2.24\times 10^{-5}$ before increasing monotonically with $\sigma$.

In Fig.7d, $\varepsilon=0.0266$, $D_{m}(fp)=D_{m}(lc)$, with \big{(}v(0),w(0)\big{)}\in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} (red curve), there is a rapid decrease in $\langle N\rangle$ for weak $\sigma$. $\langle N\rangle$ moves from $106$ to a minimum of $34.6$ within $0\leq\sigma\leq 0.5\times 10^{-5}$ before increasing monotonically with increasing $\sigma$. A quicker decrease with a lower minimum in $\langle N\rangle$ as compared to the cases in Fig.7a-c (red curves) occurs because of a smaller $D_{m}(lc)$ in Fig.7d than in all previous cases, with therefore a shorter residence time in \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}.

Still in Fig.7d with \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}) (blue curve), ISR disappears. Because at $\varepsilon=0.0266$ we have $D_{m}(fp)=D_{m}(lc)$, we explain this disappearance in terms of the other factor determining the residence time in basins of attraction i.e., the SSFs of the attractors. At $\varepsilon=0.0266$, $\lambda_{max}(fp)$ is still sufficiently high (with $\lambda_{max}(lc)<\lambda_{max}(fp)$, which means that the fixed point is more sensitive to noise than the limit cycle) and therefore, even weak-noise amplitudes have the tendency of kicking the trajectories initially in $\mathcal{B}(v_{0},w_{0})$ into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} (thereby increasing $\langle N\rangle$). As in the cases of Fig.7a and b (blue curves), one will expect that $\langle N\rangle$ increases with $\sigma\geq 0$ up to a certain maximum, and then start to decrease through the inhibitory effect of noise. This is not happening in Fig.7d (and also already in Fig.7c blue curve) firstly because $D_{m}(lc)$ is still sufficiently large (even though equal to $D_{m}(fp)$) to keep the trajectories in \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}. Secondly, and mainly because $\lambda_{max}(lc)<\lambda_{max}(fp)$ at $\varepsilon=0.0266$, which means that when trajectories get into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]}, they prefer to stay in this basin. And of course stronger and stronger noise only increases $\langle N\rangle$.

In Fig.7e, $\varepsilon=0.027673$, $D_{m}(lc)$ is much smaller and $\lambda_{max}(lc)$ is much higher than in the previous cases, but $\lambda_{max}(fp)=\lambda_{max}(lc)$. For the case \big{(}v(0),w(0)\big{)}\in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} (red curve), we therefore have a faster drop in $\langle N\rangle$, i.e., from $106$ to a minimum of $6.4$ within $0\leq\sigma\leq 0.25\times 10^{-5}$. With \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}) (blue curve) ISR disappears for basically the same reason as previously given. In this case, for $0\leq\sigma\leq 0.25\times 10^{-5}$, $\langle N\rangle$ remains at zero (since $D_{m}(fp)>D_{m}(lc)$) and as $\sigma$ increases and becomes stronger, random trajectories start to jump into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} (thus increasing $\langle N\rangle$) and remain in this basin for stronger and stronger noise with the immediate consequence of just increasing $\langle N\rangle$.

In Fig.7f, $\varepsilon=0.02785\lessapprox\varepsilon_{sn}$, $D_{m}(lc)$ is the smallest and $\lambda_{max}(lc)$ very high. For \big{(}v(0),w(0)\big{)}\in\mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} (red curve), there is a much more rapid and deeper drop in $\langle N\rangle$ with weak-noise amplitudes as compared to all the previous cases with a well defined minimum value of $\langle N\rangle=4.1$ at $\sigma=0.25\times 10^{-5}$ and then a monotonic increase in $\langle N\rangle$ with increasing $\sigma$. In this case, for some noise realizations, the number of spikes could drop down to zero. That is, weak-noise amplitudes completely terminate the spiking dynamics.

For \big{(}v(0),w(0)\big{)}\in\mathcal{B}(v_{0},w_{0}) (blue curve), because $D_{m}(fp)$ is larger at $\varepsilon=0.02785$, the residence time in $\mathcal{B}(v_{0},w_{0})$ is on average the longest for weak-noise amplitudes compared to all previous cases. We have $\langle N\rangle=0$ for $0\leq\sigma<0.30\times 10^{-5}$. For $\sigma\geq 0.30\times 10^{-5}$, the noise is now sufficiently strong to start kicking trajectories into \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} causing an increase in $\langle N\rangle$. And as $\sigma$ becomes stronger, it keeps driving the neuron and so the trajectories remain \mathcal{B}\big{[}\bar{v}(t),\bar{w}(t)\big{]} with $\langle N\rangle$ increasing monotonically with $\sigma$.

6 Conclusion

The effects of weak-noise amplitudes on the spiking dynamics of the FHN neuron model without a deterministic input current was investigated. Through bifurcation and slow-fast analyses, we determined the conditions on the parameter space for the establishment of a bi-stability regime consisting of a stable fixed point and a stable unforced limit cycle. This bi-stability regime induces a sensitivity to initial conditions in the immediate neighborhood of the separatrix isolating the basins of attraction of the attractors. Introducing noise to the system then causes transitions from the basin of attraction of the fixed point to that of the limit cycle (the neuron gets into the spiking state) and as well, from the basin of attraction of the limit cycle to that of the fixed point (the neuron gets into the quiescent state, no spiking).

We observed that in this bistable regime, weak-noise amplitudes may decrease the mean number of spikes down to a minimum value after which it increases monotonically as the noise strength increases. We showed that this phenomenon always occurred if the initial conditions were chosen from the basin attraction of the stable limit cycle.

For initial conditions in the basin of attraction of the stable fixed point, the phenomenon disappeared, unless the time-scale separation parameter $\varepsilon$ of the neuron model is bounded between $\varepsilon_{hp}=0.0250$, its Andronov-Hopf bifurcation value and $\varepsilon_{cr}=0.0260$. furthermore, the phenomenon became less and less pronounced as $\varepsilon$ increased in the interval $(\varepsilon_{hp},\varepsilon_{cr})$, and disappeared at $\varepsilon\geq\varepsilon_{cr}$.

We point out that this was not the case in Gutkin2 where the neuron model considered had both a deterministic input and a random input current. There, it was shown that the decrease to a minimum and then a monotonic increase in the spiking activity with increase in noise amplitude occur regardless of the basin of attraction from which the initial conditions are chosen from, provided the deterministic input current is above its Andronov-Hopf bifurcation value. The model in the present work has the same underlying dynamical structure as in Gutkin2 except that we do not have a deterministic input current, and only a random perturbation component is considered.

We have seen that the stochastic sensitivity functions of the stable attractors and their Mahalanobis distances from the separatrix, which both determine the length of the residence time of random trajectories in each state (quiescent or spiking state of the neuron), themselves depended on the time-scale separation parameter $\varepsilon$ of the model. From this dependence, we provided a theoretical explanation of the noise-induced phenomenon of ISR in terms of the stochastic sensitivity functions and the Mahalanobis distances of the stable attractors.

Finally, we see that the key to ISR is the multi-stability between fixed points and limit cycles, a characteristic of dynamical systems with sub-critical Andronov-Hopf bifurcations. To obtain bi-stability, in the present work, a careful relative positioning of the fixed point on the critical manifold was made. We can see in Yamakou1 how a small change in the relative position of fixed points brings about a completely different dynamical behavior in the same weak-noise limit. Plausible implications of ISR in information processing and transmission in neurons are discussed in Buchin .

Acknowledgements.

This work was supported by the International Max Planck Research School Mathematics in the Sciences (IMPRS MiS), Leipzig-Germany.