Phase dynamics of noise-induced coherent oscillations in excitable systems

Introduction.- Noise is ubiquitous in nature and is generally considered to hinder ordered behaviors of systems. However, counterintuitive phenomena in which noise brings order have also been revealed and aroused much attention in diverse fields. Indeed, noise can facilitate the formation of self-organized structures in nonequilibrium systems in physics, chemistry, and biology Prigogine2017book ; Horsthemke1984book ; Bressloff2013RMP . For example, synchronization of nonlinear oscillators usually requires mutual coupling, but common or correlated noise applied on them can induce synchronization even when the oscillators are uncoupled Zhou2002PRL ; Teramae2004PRL ; Nakao2007PRL . In nonlinear systems, stochastic trajectories under the effect of noise are not necessarily blurred versions of the deterministic trajectories Alexandrov2021PR . Noise may induce coherent trajectories that do not exist without noise, e.g., in stochastic resonance, coherence resonance, noise-induced synchronization, spatiotemporal patterns, etc. Lindner2004PR . In particular, in excitable systems, even if no periodic orbit exists in the absence of noise, coherent oscillations can still occur when noise with appropriate intensity is applied, which resemble deterministic limit-cycle oscillations.

Phase reduction is a powerful tool for reducing the dimensionality of limit-cycling systems under weak perturbations Kuramoto1984 ; Nakao2016 . Due to its simplicity and efficiency, this approach has been widely employed in analyzing various systems of coupled oscillators and also generalized to nonconventional systems such as delay-induced oscillations Kotani2012PRL ; Kotani2020PRR , reaction-diffusion systems Nakao2014PRX , stochastic limit-cycle oscillators Teramae2009PRL ; Nakao2007PRL ; Lai2011PRL , hybrid oscillators Shirasaka2017PRE ; Park2018EJAM , relaxation oscillators Izhikevich2000SIAMJAM ; Zhu2020PRE , quantum nonlinear oscillators Kato2019PRR ; Kato2020PRE , etc. The phase reduction relies on the notion of the asymptotic phase Winfree1980 ; Kuramoto1984 ; Nakao2016 of the limit cycle to characterize the dominant dynamical behaviors. However, it is not an easy task to establish a phase reduction theory for noise-induced coherent excitable systems due to the lack of a reference periodic orbit that characterizes the coherent oscillations. Efforts have so far been made mainly to develop a phenomenological phase model based on numerical simulations and data processing Han1999PRL ; Kralemann2008PRE ; Schwabedal2010PRE , or to define the stochastic version of the asymptotic phase and amplitude by solving the eigenvalue problem of the backward Kolmogorov operator for stochastic oscillators Thomas2014PRL ; Perez-Cervera2021PRL ; Kato2021Mathematics .

In this Letter, we construct a quantitative phase reduction theory for noise-induced coherent excitable systems by (i) finding a reference orbit which plays the role of the limit cycle in deterministic oscillatory cases; (ii) establishing an approximate hybrid system for calculating the phase sensitivity function; and (iii) constructing an effective phase equation and applying it to the analysis of periodically forced or mutually coupled oscillators.

Phase reduction of noise-induced coherent systems.- We consider a FitzHugh-Nagumo (FHN) system perturbed by noise applied on the fast variable as a typical example:

where $x$ and $y$ represent the fast membrane potential and slow recovery variable, respectively, and $y=f(x)=x-\frac{x^{3}}{3}$ is the nullcline of $x$. The Gaussian white noise $\nu(t)$ satisfies $\left<\nu(t)\right>=0$ and $\left<\nu(t)\nu(\tau)\right>=\delta(t-\tau)$, and $D_{\nu}$ represents its intensity. The timescale separation parameter $\varepsilon$ and the bifurcation parameter $a$ are fixed as $\varepsilon=0.0001$ and $a=1.01$. Without noise, the system (1) has a globally stable fixed point ($x_{0}$, $y_{0}$), where $x_{0}=-a$ and $y_{0}=\frac{a^{3}}{3}-a$ (see the inset in Fig. 1). When noise is applied, the system (1) can exhibit noisy but coherent oscillations due to large timescale separation, called self-induced stochastic resonance (SISR) Muratov2005PD . Figure 1 illustrates the SISR oscillations observed at noise intensity $D_{\nu}=0.01$. We refer to the system (1) as the SISR oscillator hereafter.

To determine the reference periodic orbit characterizing the coherent oscillations, we apply the distance matching condition (DMC) that we developed in Ref. Zhu2021PRR to calculate the transition positions (see Supplemental Material (SM) supplement for details). The transition position $y_{l}$ of the stochastic trajectory started from the initial state $y_{0}$ on the left branch can be determined from the following DMC:

where $f_{l}^{-1}(y)$ is the value of $x$ on the left branch at $y$ and $T_{e}(y)$ is the mean first passage time. The left-hand side of Eq. (2) represents the accumulated effect of noise on the displacement of the state away from the stable branch and $S(y)$ on the right-hand side is the distance between the middle and left branches. This condition implies that the transition occurs when the noise-induced displacement and the distance from the left to the middle branch match. As shown in Fig. 1, the transition positions on the left and right branches predicted by Eq. (2) are in good agreement with the Monte Carlo (MC) simulations. It is found that the transition position can also be predicted by considering the first passage time distribution (FPTD) Gardiner1985 ; Lim2010JCN ; Li2019PRE , which can be calculated for the left branch (and similarly for the right branch) as supplement

The peak of FPTD agrees well with the transition position predicted from DMC on each branch as shown in Fig. 1.

It is noted that the stochastic oscillations caused by SISR possess a well-defined orbit and almost keep a deterministic period. These features pave the way for applying the phase reduction approach to this system. This stochastic periodic orbit is completely different from the deterministic orbit in the oscillatory regime of the system and cannot be approximated by some limiting process of the latter. In particular, the transition from one branch to the other happens before reaching the tips of the $x$ nullcline. In order to simplify our analysis, we fix the noise intensity $D_{\nu}=0.01$ in what follows. However, we note that the SISR phenomenon can also be induced at different noise intensities as shown in Ref. Zhu2021PRR and the proposed reduction method is generally applicable to a wide range of parameters.

Considering the fast-slow characteristics of the stochastic periodic orbit, we approximate the SISR oscillator by using the following hybrid (piecewise-continuous) dynamical system:

Here, ${\bm{X}}=(x,y)$, ${\bm{F}}({\bm{X}})$ is the deterministic vector field of the system (1), ${\bm{\Pi}_{i}}$ are switching surfaces on the left and right branches, and ${\bm{\Phi}}_{i}$ are transition functions. That is, we approximate the slow stochastic dynamics along the left and right branches by the deterministic orbit of the original system and the fast dynamics between the branches by instantaneous discontinuous transitions. The transition functions and the switching surfaces can be calculated as ${\bm{\Phi}_{l}}({\bm{X}})=\left[2\cos(\varphi),y\right]^{\top}$, ${\bm{\Phi}_{r}}({\bm{X}})=\left[2\cos\left(\varphi+\frac{2\pi}{3}\right),y\right]^{\top}$, and ${\bm{\Pi}_{i}}=\left\{{\bm{X}}|L({\bm{X}})=y_{i}\right\}$, where $L({\bm{X}})=y$, $\varphi=\frac{1}{3}\arccos\left(-\frac{3}{2}y\right)$, and $y_{l}$ and $y_{r}$ are the transition positions on the left and right branches obtained by DMC supplement . This system has a stable, piecewise-continuous limit cycle, denoted as $\gamma_{s}$, of frequency $\omega_{h}$, where $\omega_{h}$ is nearly equal to the average frequency of the stochastic oscillations.

Through the above approximation, we have transformed the original stochastic system (1) to the hybrid system (4). We now apply the phase reduction method for hybrid systems Shirasaka2017PRE to further reduce the system (4) into the phase equation of the form ${\dot{\theta}(t)}=\frac{\mathrm{d}\Theta({\bm{X}}(t))}{\mathrm{d}t}=\frac{\partial\Theta({\bm{X}})}{\partial{\bm{X}}}\cdot\bm{F}(\bm{X})=\omega_{h}$, where $\omega_{h}$ is the frequency of Eq. (4) and $\theta(t)=\Theta({\bm{X}}(t))$ is the phase of the system. Here, the phase function $\Theta({\bm{X}})$ gives the asymptotic phase Winfree1980 ; Kuramoto1984 ; Nakao2016 of the state ${\bm{X}}$ within the basin of attraction of the limit cycle $\gamma_{s}$. When the SISR oscillator is additionally subjected to a weak perturbation $\bm{P}(\bm{X},t)$, the first-order approximate phase dynamics is given by ${\dot{\theta}(t)}=\omega_{h}+\frac{\partial\Theta({\bm{X}})}{\partial{\bm{X}}}\cdot\bm{P}(\bm{X},t)\approx\omega_{h}+\bm{Z}(\theta)\cdot\bm{P}\left(\theta,t\right)$, where we have approximated the system state $\bm{X}(t)$ by the state $\bm{X}_{0}(\theta(t))$ on $\gamma_{s}$ sharing the same asymptotic phase, and $\bm{Z}(\theta)=\frac{\partial\Theta({\bm{X}})}{\partial{\bm{X}}}\big{|}_{{\bm{X}}={\bm{X}}_{0}(\theta)}$ is the phase sensitivity function of $\gamma_{s}$, which can be obtained via the adjoint method for hybrid limit cycles supplement . The $y$ component of the phase sensitivity functions $Z_{y}(\theta)$ is illustrated in Fig. 2, which has discontinuities at the two transition positions.

In the above derivation of the phase equation, we omitted the stochastic fluctuations of the SISR oscillator. To better describe the stochastic dynamical behaviors and quantitatively evaluate the accuracy of our prediction, we further incorporate the stochasticity of the system into the phase equation as an effective additive noise Schwabedal2010PRE ; Nakao2010Chaos ,

where $\omega_{e}$ denotes the effective frequency of the stochastic oscillations, $\xi(t)$ is the Gaussian-white noise satisfying $\left<\xi(t)\right>=0$ and $\left<\xi(t)\xi(\tau)\right>=\delta(t-\tau)$ , and $D_{e}$ represents the effective noise intensity. The effective frequency and noise intensity are evaluated by the ensemble average Nakao2010Chaos as $\omega_{e}=\left<[\theta(t)-\theta(0)]/t\right>$ and $D_{e}=\left<\left([\theta(t)-\theta(0)]/t-\omega_{e}\right)^{2}\right>t$, respectively, by MC simulations of the original system (1), which are obtained as $\omega_{e}=2.5161$ and $D_{e}=0.0104$. As discussed in SM supplement , the approximate theoretical value of the effective frequency can be obtained from DMC as $\tilde{\omega_{e}}=\omega_{h}=2.5266$ and that of the effective noise intensity can be evaluated from FPTD as $\tilde{D_{e}}=0.0095$, which agree well with the values of $\omega_{e}$ and $D_{e}$ and quantitatively validate the hybrid system (4).

We can also evaluate the effective noise intensity by direct measurement of the phase sensitivity function supplement . The $y$ component of the phase sensitivity functions $Z_{y}(\theta)$ evaluated using several different perturbation intensities is shown in Fig. 2. As discussed in SM supplement , as the perturbation intensity $\delta$ used for the measurement becomes smaller, the mean value of the measured $Z_{y}(\theta)$ approaches the theoretical result for infinitesimal perturbation intensity calculated by the adjoint method, while its standard deviation increases as $\sigma_{g}=\delta^{-1}\sqrt{2D_{e}t}$. From the inset of Fig. 2, the effective noise intensity is evaluated as $D_{e}=0.0106$, which is also consistent with the values obtained by the other methods. In the following analysis, we fix the parameters as $\omega_{e}=2.5161$ and $D_{e}=0.0104$ in the effective phase equation (5). As we will demonstrate, the simple reduced phase equation (5) that we have derived can accurately predict the dynamical behaviors of the SISR oscillator (1) under general weak perturbations, such as the periodic forcing and mutual coupling.

Periodic forcing.- We first consider a periodically forced SISR oscillator described by

where the forcing frequency $\Omega$ is close to the effective frequency $\omega_{e}$ and $\mu$ characterizes the strength of the periodic forcing, which is weak in the sense that the frequency difference is given by $\omega_{e}-\Omega=\mu\Delta$ where $|\mu|\ll 1$ and $\Delta$ is of $O(1)$. By applying the reduction method described above, we obtain the reduced phase equation $\dot{\theta}=\omega_{e}+\sqrt{D_{e}}\xi(t)+\mu Z_{y}(\theta)\sin(\Omega t)$. By further introducing the slow relative phase $\phi(t)=\theta(t)-\Omega t$, we obtain $\dot{\phi}=\omega_{e}-\Omega+\sqrt{D_{e}}\xi(t)+\mu Z_{y}(\phi+\Omega t)\sin(\Omega t).$ Since $\phi(t)$ varies much more slowly than $\Omega t$ because $\omega_{e}-\Omega$ is of $O(\mu)$, we can average this equation via the corresponding Fokker-Planck equation over one period of fast oscillation Kuramoto1984 . This yields a further simplified phase equation

where $\Gamma_{p}(\phi)=\frac{1}{2\pi}\int_{0}^{2\pi}Z_{y}(\phi+\psi)\sin(\psi)\mathrm{d}\psi$ is a phase coupling function representing the effect of the periodic forcing on the phase dynamics.

The phase coupling function $\Gamma_{p}(\phi)$ is plotted in Fig. 3(a), which is smooth despite that the phase sensitivity function is discontinuous because $\Gamma_{p}(\phi)$ is a convolution of $Z_{y}$ with a smooth function. The solution to $\Gamma_{p}(\phi)=-\Delta$ gives the phase-locked state and it is linearly stable (unstable) when $\Gamma^{\prime}_{p}(\phi)<0$ [$\Gamma^{\prime}_{p}(\phi)>0$] if the noise is not considered. Figure 3(b) shows the time evolution of the relative phase $\phi$ converging to the stable phase-locked state for the cases with $\Delta=0.5$ and $\Delta=-0.5$. The results of MC simulations of the system (6) are in good agreement with the theoretical prediction by the noiseless phase model with the phase coupling function in Fig. 3(a). When the frequency difference is above the critical value, i.e., $|\Delta|>|\Delta_{c}|\approx 1.2596$, the relative phase will continue to increase or decrease with time because the noiseless system does not have stable phase-locked states, and phase drift will occur Pikovsky2001book . Figure 3(c) shows that this can also be well predicted by using the reduced phase model.

When $|\Delta|$ is slightly below the critical value, noise-induced phase slipping Pikovsky2001book ; Nakao2016 ; Stankovski2012PRL ; Berglund2014ArXiv can be observed. Figure 3(d) shows the results of MC simulations of Eq. (6). The relative phase $\phi$ converges to the stable phase-locked state, but it can occasionally cross the unstable phase-locked state due to noise and exhibits phase slips. This phenomenon can also be well reproduced by the reduced phase equation (7) as shown in Fig. 3(e), plotting the results of MC simulations for the averaged phase model. The discrepancy of Fig. 3(e) from Fig. 3(d) mainly originates from the averaging approximation. By directly performing simulations of the reduced phase equation before averaging, the discrepancy from the original model can be significantly reduced as shown in Fig. 3(f).

Two-coupled SISR oscillators.- Next, we consider two weakly coupled identical SISR oscillators described by

where $\mu$ ($0<\mu\ll 1$) represents weak coupling strength. The Gaussian white noise terms in Eq. (8) are mutually independent and satisfy $\left<\nu_{i}(t)\right>=0$ and $\left<\nu_{i}(t)\nu_{j}(\tau)\right>=\delta_{ij}\delta(t-\tau)$. The diffusive coupling $G_{y}(y_{i},y_{j})=\mu(y_{j}-y_{i})$ is introduced only between the slow variables. Similarly to Eq. (7), the reduced phase equation for each oscillator is given by $\dot{\theta}_{i}=\omega_{e}+\sqrt{D_{e}}\xi_{i}(t)+\mu Z_{y}(\theta_{i})G_{y}(\theta_{i},\theta_{j})$. By considering the phase difference $\phi=\theta_{1}-\theta_{2}$, which is a slow variable, and applying the averaging procedure, we can derive the equation for $\phi$ as Kuramoto1984

where $\Gamma_{d}(\phi)=\Gamma(\phi)-\Gamma(-\phi)$ is the antisymmetric part of the phase coupling function $\Gamma(\phi)=\frac{1}{2\pi}\int_{0}^{2\pi}Z_{y}(\phi+\psi)G_{y}(\phi+\psi,\psi)\mathrm{d}\psi$. The solution of $\Gamma_{d}(\phi)=0$ represents a synchronized state of the two SISR oscillators, which is stable (unstable) when $\Gamma^{\prime}_{d}(\phi)<0$ [$\Gamma^{\prime}_{d}(\phi)>0$]. As the two oscillators are identical and the coupling is symmetric, the in-phase ($\phi=0$) and antiphase ($\phi=\pm\pi$) synchronized states are always the solutions as shown in Fig. 4(a). It is notable that both synchronized states are stable in the parameter regime considered here (although the stability of the antiphase synchronization is weaker).

Since noise is present in our coupled system (8), the phase difference $\phi$ does not converge to a fixed value but forms a stationary distribution with peaks corresponding to the stable synchronized states Zhu2020ND . This distribution depends on both the noise intensity and the coupling strength. We performed MC simulations of the coupled system (8) with initial phase differences uniformly distributed in $[-\pi,\pi]$. As shown in Fig. 4(b)-(c), the phase difference tends to localize around the in-phase and antiphase synchronized states. Increasing the coupling strength can enhance the localization and more clearly separate the two states, which shows the competing relationship between the coupling-induced synchronization and noise-induced desynchronization. The distributions of the phase difference obtained by MC simulations of the original system can be well reproduced by the reduced phase equation (9) with the coupling function $\Gamma_{d}(\phi)$ obtained theoretically, wherein a higher peak is observed at the more stable in-phase state than at the less stable antiphase state. As expected, for smaller coupling strength, the phase-difference distribution is more accurately predicted.

Conclusion.- We have investigated the phase dynamics of the SISR oscillator exhibiting noise-induced coherent oscillations. The transition positions on each branch were accurately obtained via DMC or FPTD, and an approximate hybrid system was established by connecting the dynamics on the two slow branches by discontinuous transitions. We performed phase reduction on the hybrid system and obtained the reduced phase equation by further incorporating the effective frequency and effective noise intensity. The reduced equations were applied to the analysis of a periodically forced SISR oscillator and a pair of mutually coupled identical SISR oscillators. The good agreement between the predicted dynamics and the results of the original model proved the accuracy and efficiency of our reduction method. The analysis in this Letter can be readily extended to more complex situations such as nonidentical coupled oscillations and networks. Also, more accurate results would be obtained by considering higher-order approximations Rosenblum2019Chaos ; Leon2019PRE . Moreover, despite that the considered SISR oscillator has only one-dimensional slow dynamics, the present approach can also be extended to systems with higher-dimensional slow dynamics as long as the oscillation is coherent. More details will be reported in our future works.

Acknowledgments.- We thank the anonymous reviewers for valuable and insightful comments. J.Z. acknowledges support from JSPS KAKENHI JP20F40017, Natural Science Foundation of Jiangsu Province of China (Grant No. BK20190435), and the Fundamental Research Funds for the Central Universities (Grant No. 30920021112). Y.K. thanks JSPS KAKENHI JP20J13778 and JP22K14274 for financial support. H.N. thanks JSPS KAKENHI JP22K11919, JP22H00516, JPJSBP120202201, and JST CREST JP-MJCR1913 for financial support.

Supplemental Material
Phase dynamics of noise-induced coherent oscillations in excitable systems
J. Zhu, Y. Kato and H. Nakao