Stochastic periodic orbits in fast-slow systems with self-induced stochastic resonance

Nature is intrinsically noisy. While noise is often considered to be irrelevant to the realization of organized phenomena, noise does play constructive roles in many systems, e.g., gene expression [1], climate dynamics [2], neural spiking [3], etc. In particular, resonance induced by noise is one of the most intriguing phenomena and is not only of theoretical interest but also of practical importance. Among them, stochastic resonance (SR) [4] attracts broad attention due to its counterintuitive behavior, which was first observed in the climate system [5]. In stochastic resonance, moderate intensity of noise can enhance the signal-to-noise ratio for periodically perturbed random systems. It has been extensively investigated both theoretically and experimentally in all kinds of disciplines [4, 6, 7, 8]. In excitable systems, even without the external periodic signal, the coherence resonance (CR) [9, 10, 11] manifests itself as coherent oscillations for intermediate noise. The mechanism of CR has been clarified by Pikovsky and Kurths via the first passage statistics [11]. Later, Lindner and Schimansky-Geier [12] considered an analytical approximation in the limit of large timescale separation and quantitatively described CR through the marginal probability density of the reduced system. They applied similar ideas to the piecewise linear FitzHugh-Nagumo system for both CR and SR [13]. More accurate results, such as the rate, coefficient of variation and diffusion coefficient, were obtained in the leaky integrate-and-fire model [14]. For detailed descriptions and quantitative approaches to the coherent behaviors in excitable systems, the readers are referred to Ref.[15]. On the other hand, Muratov et al. [16, 17] found that noise on excitable systems with largely different timescales can provide distinct mechanisms of coherence. For noise acting on the fast variable, the resulting coherent oscillation was named self-induced stochastic resonance (SISR) [16, 17]. SISR has been shown to account for the mixed-mode oscillations in a relaxation oscillator [18]. Besides, SISR may also explain the mechanism of anticoherence resonance for large noise [19, 20]. Recently, Yamakou et al. investigated SISR in multiplex neural networks and showed its potential application to optimal information processing [21, 22].

Determining the oscillator’s periodic orbit is important for capturing its dynamical behaviors under perturbations. For deterministic limit-cycle oscillators, the asymptotic phase can be correspondingly defined within the basin of attraction of the limit cycle and the phase reduction approach [23, 24] can be applied to reduce the dimensionality of the system for weak perturbations. However, it is not an easy task to determine the periodic orbit for stochastic oscillators. Some efforts have been made on the CR oscillator [25]. The noise-induced escape from the fixed point can be approximated by a jump process. The remaining part of the stochastic trajectory will be dominated by the noise-free system. Through this approximation, the phase reduction approach can be correspondingly established on this hybrid system [25]. Compared with the CR oscillator, the SISR phenomenon is more robust to parameter variations as the latter does not require the system to be close to bifurcation [16]. This gives SISR broader coherence ranges with respect to the noise strength. To depict the stochastic periodic orbit, Muratov et al. [16, 17] applied the large deviation theory to find the critical transition position via the timescale matching condition. Recently, Yamakou and Jost [26] further clarified the allowed interval for the potential difference and discussed the connection between inverse stochastic resonance and SISR. These analyses required both the slow timescale and the noise strength to approach zero.

However, when the system has different timescales on different branches of slow manifolds as in the FitzHugh-Nagumo (FHN) model, there can be asymmetry in the noise-induced transition as shown later. In this study, to accurately obtain the stochastic periodic orbit in SISR oscillators, we define the mean first passage velocity and propose a condition that determines the critical transition position under the assumption that the transition process is continuous. This assumption differs from the previous studies, where the transition is considered as instantaneous and determined via the timescale matching condition [16, 17]. We first employ this condition in a simplified piecewise linear FHN system and then extend it to the original FHN model. The theoretical results exhibit good consistency with the numerics even for large noise strength. Thus, the method to determine the periodic orbit proposed in this paper is robust and may have a wide range of applications.

The structure of this paper is as follows: In Sec. II, the piecewise linear FitzHugh-Nagumo system is introduced and its dynamical feature is briefly discussed. Next, the self-induced stochastic resonance of this system is investigated in Sec. III and the distance matching condition is proposed through which the critical transition positions on the left and right branches can be identified. Further, the theoretical stochastic periodic orbit is verified by Monte Carlo simulations. In Sec. IV, the proposed method is extended to the original FitzHugh-Nagumo system. Finally, conclusions and discussions are given in Sec. V.

To illustrate the main idea of this paper, we choose the piecewise linear FitzHugh-Nagumo (PWL-FHN) system, also known as the McKean model [27]. This model preserves the essential characteristics of neuronal dynamics while assures computational simplicity. The governing equations are as follows:

and

where $x$ and $y$ represent the membrane potential and recovery variable, respectively. The timescale separation parameter $\varepsilon$ is assumed small, which is key to the realization of SISR. We set $\varepsilon=0.05$ and $a=0.95$ fixed in this paper. Since the bifurcation parameter $a<1$, system (1) has a stable limit cycle. This is different from the previous studies [16, 17] in which the bifurcation parameter was taken in the excitable regime. Actually, whether the system is excitable or oscillatory is not important as the SISR phenomenon occurs away from the equilibrium or the deterministic limit cycle [16, 17]. Considering the oscillatory situation allows SISR to be observed for small enough noise, which is helpful for the validation of the theoretical prediction (results for the excitable situation are also given in Appendix A, which are similar to those for the oscillatory situation discussed in the main text). The timescale separation can be observed in the vector field and the deterministic trajectories as shown in Fig. 1. The trajectory quickly converges to the left or right branches of the $x$-nullcline. The period of the limit cycle is dominated by the slow motion along them, which is determined by the slow subsystem $y$ in system (1). We denote the left, middle and right branches of the $x$-nullcline as $B_{l}$, $B_{m}$ and $B_{r}$, respectively. For fixed $y\in(-5,5)$, the fast subsystem $x$ has three equilibrium points, two of which on $B_{l}$ and $B_{r}$ are stable and the other on $B_{m}$ is unstable.

The potential value $U(x;y)$ of the fast subsystem $x$ of system (1) can be easily calculated as $U(x;y)=-\int f_{\rm pwl}(x)-y\,dx+C$ since it is piecewise linear, where $y$ is considered as a parameter and $C$ is a constant to be determined. Without loss of generality, we set the potential at the tips of the $x$-nullcline (i.e., $x=\pm 1$) to be zero. Figure 2(a)-(b) illustrate the potential $U(x;y)$ for fixed $y$ values and on the three branches of the $x$-nullcline, respectively. The potential difference between the middle branch and the left or right branch is monotonically decreasing with the system’s evolution. They are symmetric to the line $y=0$ reflecting the symmetry of the $x$-nullcline.

As shown by Muratov et al. [16, 17], the SISR is observed when the noise acts on the fast variable. Therefore, we consider a stochastic version of the PWL-FHN system given as follows:

where $\xi(t)$ is Gaussian white noise with zero mean satisfying $\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(\tau)\rangle=\delta(t-\tau)$. The parameter $\sigma$ determines the strength of the noise. Figure 3 illustrates the phase diagram and the timeseries for $\sigma=0.1$, which clearly shows the coherent behavior with a nearly deterministic period. It is notable that the stochastic periodic orbit caused by the noise is smaller than the deterministic limit cycle and that the orbit is not symmetric, namely, the transition to the other branch occurs earlier on the left branch than on the right branch. This asymmetry cannot be captured by the conventional timescale matching condition [16, 17, 26] as we will discuss below. The symmetric transition can only be realized for $a=0$, where the timescales on the two branches are the same owing to the symmetric vector field.

To determine the critical transition position on each branch, the conventional timescale matching condition assumes that the mean first passage time (MFPT) of boundary crossing and the recovery timescale matches with each other at the transition position. Following Muratov et al. [16], the timescale matching condition is given by:

where $\Delta W$ stands for the quasipotential difference between the middle and the left or right branches. The quasipotential $W$ defined by Freidlin and Wentzell [28] is similar to the potential function in the deterministic system, and the quasipotential difference $\Delta W$ characterizes the difficulty for the state to escape from a stable equilibrium. The left hand side (LHS) of Eq.(4) denotes the timescale of noise-induced transition of the fast variable while the right hand side (RHS) represents the timescale of recovery of the slow variable. By calculating the quasipotential difference of the fast subsystem by considering $y$ as a constant parameter, the critical position $y_{*}$ can be obtained, which predicts the transition position on each branch. Equation (4) is accurate in the vanishing noise limit and the slow relaxation limit, i.e., $\sigma\to 0$ and $\varepsilon\to 0$. Thus, according to this criterion, the critical positions would be symmetric with the same absolute value $|y_{*}|$ on the left and right branches. However, because $\varepsilon$ takes a finite non-zero value, this prediction is inconsistent with the observation as in Fig. 3.

To predict the noise-induced transition more accurately, we propose the following criterion. Using the MFPT $T_{\rm e}(y)$, we define the mean first passage velocity (MFPV) $V_{\rm e}(y)$:

where $S(y)$ represents the distance of the left or right branch on which the state exits from the boundary (the middle branch). We expect that the transition occurs at $y_{*}$ at time $t(y_{*})$ where the integration of the MFPV gives $S(y_{*})$, i.e.,

where $y_{0}$ is the starting position at time $t(y_{0})$ of the trajectory on the branch under consideration. Equation (6) could serve as a more precise distance matching condition than Eq.(4) for solving the critical transition position $y_{*}$. This expression takes into account the slow variation in $y$ during the transition process; the LHS represents the approximate accumulation effect on the displacement of the fast variable $x$ as the slow variable $y$ moves from $y_{0}$ to $y_{*}$, and we assume that this accumulated displacement matches the distance $S(y_{*})$ from the boundary at the transition position $y=y_{*}$.

We employ this idea for the PWL-FHN system. For simplicity, we consider the left branch only. The procedure for solving the critical position on the right branch is similar. The mean first passage time can be obtained from the well-known Kramers rate[29, 30]:

where $U_{\rm m}(x_{y})$ and $U_{\rm l}(x_{y})$ represent the potential function at the middle and left branches of the deterministic fast subsystem of system (1), respectively. The double prime denotes the second-order derivative with respect to the variable $x$, where the subscript $y$ indicates the $y$-position for the calculation of the potential function. Since system (1) is piecewise linear, $U^{\prime\prime}_{\rm m}(x_{y})$ and $U^{\prime\prime}_{\rm l}(x_{y})$ can be easily obtained as $-5$ and $10$, respectively. Because the fast subsystem is one-dimensional, the quasipotential difference $\Delta W=2\Delta U$, where $\Delta U=U_{\rm m}(x_{y})-U_{\rm l}(x_{y})=dU_{\rm ml}$ (see the inset in Fig. 2).

The distance can also be readily obtained by $S(y)=f_{\rm m}^{-1}(y)-f_{\rm l}^{-1}(y)$ for the left branch ($S(y)=f_{\rm m}^{-1}(y)-f_{\rm r}^{-1}(y)$ for the right branch), where the subscript indicates the specific branch (i.e., m, l and r for middle, left and right, respectively). Here, we have omitted the original subscript $\rm pwl$ for simple notification.

By substituting the distance function $S(y)$ and the MFPT Eq.(7) into Eq.(6) and by using that $dy=\varepsilon(x+a)dt$ and $x=f^{-1}(y)$, we can rewrite the condition for the critical transition position $y_{*}$ as

where we approximated $y_{0}$ by $y_{l}=5$ for the left branch. (A more appropriate, self-consistent choice of $y_{0}$ is the transition position on the other branch as obtained in Appendix B via the convergent method. However, the critical transition positions obtained by this method are nearly the same as those obtained by using the above two fixed values of $y_{0}$ as shown in Fig. 13. For simplicity, we fix $y_{0}=5$ and $y_{0}=-5$ on the left and right branches, respectively, in the rest of the main text.) By solving Eq.(8), we can obtain the critical transition position $y_{*}$. For illustration, we set $\sigma=0.1$. The values of the LHS and RHS of Eq.(8) are shown in Fig. 4(a) as functions of $y$. From the intersection of the two curves, the critical transition position on the left branch is evaluated as $y_{*}\approx-3.983$. Similarly, the critical transition position on the right branch is $y_{*}\approx 4.392$ (see Fig. 4(b)).

It can be easily observed that the MFPV (thus the accumulated displacement in the direction of the fast variable $x$) remains nearly zero from the beginning ($y=y_{l}$) and quickly increases when the state approaches the transition position. This explains the mechanism of very coherent oscillations for SISR: before the transition position, the MFPV is nearly zero; after the transition position, the MFPV increases rapidly, which results in the noise-induced transition. The degree of coherence is mainly influenced by the timescale separation parameter $\varepsilon$. The smaller $\varepsilon$ is, the quicker the distance (LHS of Eq.(8)) increases with $y$, and this leads to steeper curves in Fig. 4(a). Therefore, a small range of $y$ will determine the transition phenomenon. Furthermore, the left branch is closer to the $y$-nullcline than the right branch, which makes the evolution along the former slower than that along the latter. This causes the slower increase of LHS of Eq.(8) on the right branch, thus the transition on the right branch occurs later than on the left branch as shown in Fig. 4 (this can also be quantitatively understood from the term $f^{-1}_{l}(y)+a$ in Eq.(8)). In the limiting situation, i.e., for $\varepsilon\to 0$, as the other parts of the integrand in the LHS of Eq.(8) remains finite, $\varepsilon\,{\rm exp}\left(\frac{\Delta W}{\sigma}\right)$ should be of order $O(1)$, which is in fact the result obtained by Muratov et al. [16, 17] as in Eq.(4).

Once we obtain the critical transition positions, the other parts of the stochastic periodic orbit can be determined by the deterministic slow dynamics along the left and right branches. Therefore, by connecting the transition positions, the total stochastic periodic orbit is obtained as shown in Fig. 5. It can be seen that the predicted stochastic periodic orbit is in excellent agreement with the stochastic trajectory obtained by the Monte Carlo simulation. Besides, it also predicts the asymmetric transitions on different branches. Following the same procedure, the critical transition position $y_{*}$ can be calculated for other noise strengths. We compute the results for $\sigma=0.01,0.02,0.05,0.1,0.2,0.5,1,2$ and $5$, which are displayed in Fig. 6. The curves for the left branch are steeper compared with those on the right branch, which predicts earlier transition and higher coherence on the left branch. Note that, although we calculated the result for $\sigma=5$, it is practically impossible to achieve the corresponding transition positions since for $\sigma=5$ the critical position $y_{*}$ on the left branch is larger than that on the right branch.

To be more intuitive, the theoretical stochastic periodic orbits and the results of Monte Carlo simulations for various noise strengths are displayed in Fig. 7 (similar results for the excitable PWL-FHN system are given in Appendix A). It can be seen that the theoretical predictions are in good agreement with the simulation results even for large noise, although the degree of coherence is deteriorated for increasing the noise strength. The stochastic periodic orbit becomes smaller as the noise strength $\sigma$ is increased, and at $\sigma_{c}\approx 2.733$, the critical transition positions on both branches become equal to each other. In this case, the stochastic periodic orbit will shrink into a line segment and the period will approach zero. The corresponding critical transition position is $y_{*c}\approx 0.51$, which is different from $y_{*c}=0$ in the case of the vanishing timescale separation ($\varepsilon\to 0$) [26]. This prediction may account for the large-noise asynchrony in interacting excitable systems reported in Ref. [31].

The period of the stochastic oscillation, a fundamental quantity for its characterization, can also be predicted by using the proposed criterion. Because of the timescale separation, the period is dominated by the slow motion along the left and right branches. Therefore, the period of the stochastic periodic orbit can be approximated by the following equation:

where $y_{l}$ and $y_{r}$ represent the predicted critical transition positions on the left and right branches, respectively. The theoretical period of the stochastic periodic orbit for different noise strengths together with the results of Monte Carlo simulations are illustrated in Fig. 8. We can see that the prediction is in excellent consistency with the simulation results even for large noise. Besides, the small standard deviation demonstrates that the coherent oscillations persist over a broad range of noise strength, which is a significant feature of SISR that cannot be observed in CR [16].

To verify the validity of our theory in a more realistic model, we apply the proposed condition to the original FitzHugh-Nagumo system, i.e.,

where parameters are $\varepsilon=1{\rm e}-4$ as in Ref. [17] and $a=0.8$. The Gaussian white noise $\xi(t)$ is the same as that used in system (3). Following the previous procedure, the theoretical stochastic periodic orbits for different noise strengths can be numerically obtained via Eq.(8). They are illustrated in Fig. 9, which also shows excellent consistency with the results of Monte Carlo simulations. In the numerical simulations, to reduce the computational cost, we rescaled the equations by introducing a slow time $\tau=\varepsilon t$, which yields the noise strength multiplied with a coefficient $\sqrt{\varepsilon}$, i.e.,

In a similar way to the PWL-FHN system, the period of the stochastic periodic orbit of the FHN system can be calculated as in Fig. 10. Again, the theoretically predicted period of the stochastic orbit in the FitzHugh-Nagumo system (11) shows excellent agreement with the results of Monte Carlo simulations.

In summary, we have investigated the stochastic periodic orbit in the SISR oscillator. By introducing the notion of the mean first passage velocity, the distance matching condition is established, through which the critical transition position on each branch can be determined. The theoretical stochastic periodic orbit is obtained by gluing the motions along the slow manifolds. By comparison with the results of Monte Carlo simulations, the theoretical predictions of the critical transition positions and periods of stochastic oscillations are proved to have good accuracy even for large noise. One of the main differences between the timescale matching condition by Muratov et al. [16] and our proposed condition is that the transition in the former is assumed to be instantaneous while it is considered to be continuous in the latter. As observed in Figs. 7 and 9, the fluctuations of the stochastic trajectories along $x$-axis gradually increase as they approach the transition point. This is due to slow change in the slow variable during the transition process and our condition approximately incorporates this effect.

From the derived distance matching condition Eq.(8), we can see that the timescale separation parameter $\varepsilon$ is significant to the observation of the SISR phenomenon. Because it controls the slope of the integration curve given by the LHS in Eq.(8), the smaller $\varepsilon$ corresponds to the larger slope, thus induces more coherent stochastic oscillations. This mechanism is briefly discussed in Appendix C. However, for the effect of noise strength on the degree of coherence, the slope of the integration curve does not give similar conclusions. This can be noticed by comparing the results in Fig. 6 and Fig. 8, where larger slopes don’t imply higher coherence. This is because increasing noise strength will not only increase the slope of the integration curve (this can be inferred from Eq.(8), similar to $\varepsilon$), but also increase the amount of fluctuations. The competition between these two impacts determines the effect of noise strength on the degree of coherence.

The definition of the mean first passage velocity in Eq.(5) is a natural extension of the mean velocity in the deterministic system. Instead of Eq.(8), we may also consider different but similar conditions. For example, the distance function $S(y)$ can be replaced by a constant 1 in Eq.(8), then the mean first passage velocity will be replaced by the mean passage frequency and the corresponding matching condition will give the critical position for the first transition. Another possibility is to replace $S(y)$ by the quasipotential difference $\Delta W(y)$, though the physical meaning may be not so intuitive. We have tested these two matching conditions and obtained similar predictions (although not as good as the distance matching condition in Eq.(8) as shown in Appendix D, wherein the comparison with the results via the timescale matching condition (4) is also displayed).

For higher dimensional systems, the computation of the quasipotential would be a major obstacle. Recent advances in the numerical methods for solving the quasipotential [32, 33] would offer the possibility for the application of the proposed method to more complex situations, such as the prediction of inverse stochastic resonance [34].