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

Noise is ubiquitous in various systems. In systems with multiple timescales, noise can induce various coherent behaviors. Self-induced stochastic resonance (SISR) is a typical noise-induced phenomenon identified in such systems, wherein noise acting on the fast subsystem causes stochastic resonancelike boundary crossings. In this paper, we analyze the stochastic periodic orbits caused by SISR in fast-slow systems. By introducing the notion of the mean first passage velocity toward the boundary, a distance matching condition is established, through which the critical transition position of boundary crossing can be calculated. The theoretical stochastic periodic orbit can be accordingly obtained via gluing the dynamics along the slow manifolds. It is shown that the theoretical predictions are in excellent agreement with the results of Monte Carlo simulations for a piecewise linear FitzHugh-Nagumo system even for large noise. Furthermore, the proposed method is extended to the original FitzHugh-Nagumo system and also found to exhibit consistent accuracy. These results provide insights into the mechanisms of coherent behaviors in fast-slow systems and will shed light on the coherent behaviors in more complex systems and large networks.

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 noiseinduced 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.

II. PIECEWISE LINEAR FITZHUGH-NAGUMO SYSTEM
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 ε is assumed small, which is key to the realization of SISR. We set ε = 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 The potential value U(x; y) of the fast subsystem x of system (1) can be easily calculated as U(x; y) = − f 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  and U r = − y 2 20 + 3y 2 − 25 4 denote the potential for the left, middle and right branches, respectively). The inset shows the potential differences between the middle branch and the left (dU ml = U m − U l ) or right (dU mr = U m − U r ) branch, respectively. potential at the tips of the x-nullcline (i.e., x = ±1) to be zero. Figure 2

III. STOCHASTIC PERIODIC ORBIT OF PWL-FHN SYSTEM
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 where ξ(t) is Gaussian white noise with zero mean satisfying ξ(t) = 0 and ξ(t)ξ(τ ) = δ(t − τ ). The parameter σ determines the strength of the noise. Figure 3 illustrates the phase diagram and the timeseries for σ = 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 ∆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 ∆W char-acterizes 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., σ → 0 and ε → 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 ε 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 e (y), we define the mean first passage velocity (MFPV) V 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., t(y * ) 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 m (x y ) and U 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 ′′ m (x y ) and U ′′ l (x y ) can be easily obtained as −5 and 10, respectively. Because the fast subsystem is one-dimensional, the quasipotential difference ∆W = 2∆U, where Fig. 2).
The distance can also be readily obtained by 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 pwl for simple notification.
By substituting the distance function S(y) and the MFPT Eq. (7) into Eq.(6) and by using that dy = ε(x + a)dt and x = f −1 (y), we can rewrite the condition for the critical transition position y * as  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 noiseinduced transition. The degree of coherence is mainly influenced by the timescale separation parameter ε. The smaller ε 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 ε → 0, as the other parts of the integrand in the LHS of Eq.  strengths. We compute the results for σ = 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 σ = 5, it is practically impossible to achieve the corresponding transition positions since for σ = 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 σ is increased, and at σ c ≈ 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 ≈ 0.51, which is different from y * c = 0 in the case of the vanishing timescale separation (ε → 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].

IV. STOCHASTIC PERIODIC ORBIT OF FITZHUGH-NAGUMO SYSTEM
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 ε = 1e − 4 as in Ref. [17] and a = 0.8. The Gaussian white noise ξ(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 τ = εt, which yields the noise strength multiplied with a coefficient √ ε, 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.

V. CONCLUSIONS AND DISCUSSIONS
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 in the former is assumed to be instantaneous while it is considered to be continuous in the  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. 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].

ACKNOWLEDGMENTS
We thank the anonymous reviewers for valuable and insightful comments, which have In the main text, it is explained that whether the system is excitable or oscillatory is not important for our prediction. However, since the original work by Muratov et al. [16] considered the excitable case, it is interesting to see that the proposed distance matching condition works also for this situation. We consider the PWL-FHN system (3) and set the bifurcation parameter at a = 1.05 while keeping the other parameters unchanged as in the main text. The stochastic periodic orbits predicted by utilizing the distance matching condition (8) are shown in Fig. 11, which are also in good agreement with the Monte Carlo simulation results.
Unlike the oscillatory case, for the excitable system, when noise is small enough, the state will converge to the equilibrium before the transition. There is a lower bound of noise strength(σ ≈ 0.01, see Fig. 12), below which the stochastic orbit is no longer coherent or even difficult to initiate spikes(see Fig. 11 for σ = 0.001). We here calculate these self-consistent transition points and compare the results with our approximation.
To obtain the final critical transition positions on both branches, we can apply a convergent method as follows. First, we start at y 0 = 5 on the left branch and obtain the candidate transition position y * l . We utilize this value as the starting position y 0 on the right branch and calculate the candidate transition position on the right branch y * r . Then, y * r will be chosen as the starting position on the left branch to compute the next candidate transition position. We continue this process until the difference between adjacent candidate transition positions on each branch is less than a tolerance value (e.g., 0.001).  Appendix C: Influence of the timescale separation parameter on the degree of coherence It is intuitive that decreasing the timescale separation parameter ε increases the degree of coherence of the SISR oscillator. Figure 14 illustrates the influence of ε on the integration curve, i.e., LHS of Eq.(8) of the PWL-FHN system (3) for the left branch. It is shown that smaller ε corresponds to a larger slope (absolute value) at the intersection point. For the same range of distance ∆S, a larger slope implies a smaller range ∆y. This means that a smaller range of y determines the transition process, therefore the stochastic orbits exhibit higher coherence. Similar results also hold for the right branch of Eq. (8). Note that this mechanism is for fixed noise strength. The influence of noise strength on the degree of coherence is complex as explained in the main text. σ ln(ε −1 )). It can be seen that the distance matching condition used in the main text is more accurate than the others.