Coherence and stochastic resonance in a birhythmic van der Pol system

We consider the response to uncorrelated noise and harmonic excitation of a birhythmic van der Pol-type oscillator. This system, as opposed to the standard van der Pol oscillator, is characterized by two stable orbits. The noisy oscillator can be analytically mapped, with the technique of stochastic averaging, onto an ordinary bistable system with a bistable (quasi)potential. The birhythmic oscillator can also be numerically characterized through the diagnostics of coherent resonance and the signal-to-noise-ratio. The analysis shows the presence of noise-induced coherent states, influenced by the different time scales of the oscillator.


Introduction
Self oscillatory systems exhibit spontaneous periodic oscillations, i.e. a stable orbit in the phase space. More rarely, one encounters birhythmicity, or the contemporary presence of two stable orbits (limit cycles) for the same set of parameters, as in some biochemical systems [1] or circadian oscillations [2], cell populations [3,4], neuronal dynamics [5], protein dynamics [6]. To model the oscillations, one of the first, and still nowadays prototypal system, is the van der Pol oscillator [7], that has been employed for biological modelling [8] and other oscillations [9][10][11]. Quite naturally, van der Pol oscillators have been generalized with a higher order polynomial dissipation (that entails multiple stable attractors with different natural frequencies) to be employed as a paradigm for birhythmicity, for example in some enzymatic reactions [10,[12][13][14][15][16][17][18], or energy harvesting [19]. Noise, that of course is present in real processes, perturbs the periodicity of van der Pol oscillators. In the case of birhythmic systems, noise also induces an hysteretic behavior: the displayed limit cycle depends upon the initial conditions. The transition from an attractor to another is essentially governed by an Arrhenius-type law [20][21][22] and the effective energy barrier that regulates the escapes depends upon the fluctuation spectrum and the signal (or perturbation) shape. Thus, a van der Pol-type birhythmic system shares the a e-mail: ryamapi@yahoo.fr main features of bistable systems [10,23,24]: the presence of two metastable states, an exponential dependence of the lifetime upon the inverse noise intensity, and a mean first passage time across an unstable separatrix. The similarity can be ascribed to the possibility to reduce the averaged system (averaged over the period of the self oscillations) to a standard bistable system governed by a double well quasi-potential (or pseudopotential) [25,26]. This analysis could also take advantage of the nonequilibrium potential constructed on the basis of an extremum principle [27][28][29]. In fact, it can be postulated that the behavior of the escape times, or more accurately of the mean first passage time through the separatrix between the two attractors, exponentially depends upon a quantity (the quasi-potential, or the effective energy barrier) and it is inversely proportional to the noise intensity [12,18,30]. The quasi-potential approach has also been extended to correlated noise case both driven [31] and undriven [32], that also exhibits stochastic bifurcations [33]. We have made theoretical estimations of the noise thresholds for escape times in van der Pol-type birhythmic system.
We propose to add a sinusoidal drive to such effective bistable potential, i.e. we consider a forced van der Pol type birhythmic system [34] as in some biological systems that are characterized also by a forcing term [35]. When a periodic drive is added to a static, ordinary bistable system characterized by a double well, it is known that the random oscillations between the two (meta)stable solutions induced by noise and the period of the external drive can co-operate, producing states that are, in some statistical sense, coherent with the external drive. This is the property exploited in Stochastic Resonance (SR) [36], stochastic signal detection [37,38], and dynamical phase transitions [39]. The interesting phenomenon of SR may appear in a system subject to both random and periodic force [40]. In particular, stochastic resonance is the term used to describe systems in which the presence of input or internal noise provides the optimal output of that system [41,42]. In a nonlinear system, SR ordinarily occurs under three conditions: the presence of a bistable nonlinear system, an applied periodic signal, and the presence of noise. Under these conditions, the response of the system can be similar to the resonance behavior, which is why it is called stochastic resonance. In the van der Pol birhythmic case, the potential is not an actual potential, but a quasi-potential that governs the rare escapes. Moreover, the system has several time scales: the periods of the orbits, the period of the external drive, and the average escape time from the attractors, while in the ordinary analysis of SR the first time scale (the period of the orbits) is absent. It is therefore not obvious that SR can occur, and this is the objective of our research: to investigate if a birhythmic system, where the two states are periodic orbit with an intrinsic time scale, can exhibit a (stochastic) resonance in the presence of noise, as the analogous ordinary bistable systems. Also, we wish to verify if the reduction of the model equation by means of the quasi-potential to a bistable system can also describe the co-operation between the external deterministic drive and the noise effects.
The paper is organized as follows. We first summarize some known properties of van der Pol birhythmic oscillators: in Section 2 we briefly describe the physical model of the periodically driven multi limit-cycle van der Pol oscillator. Section 3 deals with the diagnostic of the coherence between noise and deterministic oscillations in the stochastic birhythmic van der Pol system, to formalize the main tools employed in this work. The original results begin in Section 4, that deals with stochastic resonance-like phenomenon on the birhythmic van der Pol system. In this Section we numerically evaluate the tools employed to quantify the degree of coherence (or anti-coherence). In Section 5, we show that a quasi-potential can be approximately derived, also including the sinusoidal drive. Moreover, we numerically show that the system actually displays SR in the coherent detection mode, at least for some set of parameters. The last section leads to the conclusions.

Model of a noisy driven birhythmic system
In this section we set the stage for the studies of the van der Pol like birhythmic system that will be analyzed through the paper. After the model equations have been presented in Section 2.1, it will be shown that the system can be approximately mapped onto a bistable potential for the radius A of the orbits in Section 2.2.

The model of noisy driven birhythmic system
A stochastic version of a modified van der Pol system that exhibits birhythmicity, is described by the following Langevin equation [12,18,43]: Here overdots denote time derivative. The term Γ (t) represents an additive Gaussian white noise with amplitude D [40,44]: The quantities α, β are positive parameters, μ is the parameter of nonlinearity which is restricted to small values, μ 1. Adding an external deterministic drive, equation (1) becomes: where E 0 and ω are the amplitude and the frequency of the periodic force, respectively, while the properties of the stochastic term are governed by equation (2).
Equation (3) can be reduced to a nonlinear bistable oscillator in a potential for E 0 = 0 [12,14,15,17,18,43]. However, here we assume that the signal E 0 sin ωt is small enough (E 0 1), that, in the absence of noise, it cannot force the particle to move from one well to the other. Moreover, we also assume that the periodic signal is slow enough to be considered adiabatic, respect to the time scale of the self excited oscillations (i.e. ω 1). The system is a nonlinear self-sustained oscillator which possesses more than one stable limit-cycle solution when Γ (t) = 0 [2,12,18], essential to describe some biological processes. The system exhibits super-harmonic resonance structure, and symmetry-breaking crisis and intermittency. When employed to model biochemical systems like enzymatic-substrate reactions, x in equation (1) is the population of enzyme molecules in the excited polar state [13,15]. The parameters α and β measure the degree of tendency of the system to a ferroelectric instability, and μ is a parameter that effectively refers to strength of nonlinear damping (we refer to [12,18] for more details).

Analytic considerations
The standard approach to a van der Pol type equations is to approximate the periodic solutions of the free-noise equation (1) with The amplitude A has been found to be independent of the coefficient μ and implicitly given by the relation [18] 5β 64 The coefficient μ enters in the expression for the frequency Ω, which is implicitly given by the relation where ω 2 is a function of the amplitude A and of the parameters α and β [12,18]. Equation (5), depending on the value of the parameters α and β, gives rise to one or three positive real roots. Thus, for some parameters value, the system exhibits three limit cycle solutions (two stable and one unstable) [12,18,45,46], with different associated frequencies, that is the essence of birhythmicity. The unstable limit cycle represents the separatrix between the basins of attractions of the two stable limit cycles. The shaded area in Figure 1 denotes the region of the parameters α and β where birhythmicity occurs [45,46]. Assuming that the noise intensity is small [18], in the quasiharmonic regime the modified van der Pol equation (1) reduces to the following effective equation that depends on a (quasi)potential U [25,26]: whereD = D/ω 2 and the effective potential U (A) is given by [31] U (A) = μ 128 (8) The barrier of the effective potential ΔU characterizes the escapes from the attractors with limit-cycles A 1 and A 3 . In fact, a quasipotential is a Lyapunov function that summarizes the low noise properties of the system with an Arrhenius-like behavior, T ∝ exp(ΔU/D) [12,18]. In a bistable system, two barriers characterize the exit times from the two stable states: ΔU 1 and ΔU 3 for the escape from the orbit of radius A 1 and A 3 , respectively. As an example, in Tables 1 and 2 we report the energy barriers [12,30] of the multi-limit-cycle associated to the frequencies and amplitudes of the model for some values of the physical parameters α and β. In the shaded region of the parameters plane (α, β) of Figure 1, where two global minima appear, the potential U is symmetric (Tab. 1, ΔU 1 ΔU 3 , the black line of Fig. 1) or asymmetric (Tab. 2, ΔU 1 ΔU 3 or ΔU 1 ΔU 3 , the light gray region of Fig. 1) with respect the unstable amplitude A 2 [12,18].
To calculate the statistics of the exit time, i.e. the escape time of the particle from a minimum of the potential, we derive an evolution equation for the Probability Density Function (PDF) of the variable amplitude A(t). The Fokker-Planck equation corresponding to the Langevin equation (7) with δ correlated Gaussian noise sources ζ 1 (t) reads [31,47] In equation (9), P (A, t) is the PDF of the stochastic process of the limit cycle amplitude A(t). The stationary solution P (A) = P (A, t → ∞) undergoes a transition from a bimodal to a unimodal (or the opposite, from unimodal to bimodal) distribution by increasing the noise intensity D (see Figs. 5-7 in Ref. [43]). In the absence of additive noise, D = 0, the system is conned in one semi-axis (0 < A < A 2 or A 2 < A, according to the initial condition). Equation (7) gives an estimate for the change in the energy U (A) over the period 2π/ω in the quasiharmonic regime. The escape process is depicted in Figure 2 as a jump over an activation barrier ΔU 1→3 (from the leftmost minimum A 1 ) and ΔU 3→1 (from the rightmost minimum A 3 ) [12,18]. The characteristics height can be controlled by the variation of the parameters α and β [30]. The noisy birhythmic van der Pol equation (1) can be characterized through the distribution P (T ) of the escape times (denoted by T i ) from the two wells of potential U . The probability for the system to hop between the potential wells is defined through noise dependent Kramers rate [47,48], which is the inverse of the average escape time T . The quantities T 1,3 [18,30] for small noise intensity D < ΔU i are given by (here the primes refer to differentiation respect to the radius A): that describes the transition of the system from attractor with limit-cycle amplitude A 1 (left side potential well) to attractor with limit-cycle amplitude A 3 (right side potential well). A similar equation holds for the reverse passage, i.e.
The quantities (10) and (11) measure the relative stabilities pertaining to the attractors with limit cycle amplitude A 1 and A 3 through the resident time given by the relation [18,30]: The characteristics of the stability properties through equation (12) in a modified van der Pol oscillator (1) are strongly influenced by both the nonlinear coefficients α, β and the noise intensity D [30]. An escape is counted if the amplitude of the system A is greater than the separatrix amplitude A 2 for left to right escape. Analogously, if the amplitude is less than the separatrix one counts the reversal passage from right to left. In the absence of noise, the system would remain confined to its initial state. In the presence of noise, transitions eventually occur, and the two-state output evasion is a random distributed sequence. However, a birhythmic system with an intrinsic time scale there is the possibility that the noise cooperates with the internal time scales. The measure of the degree of cooperation is the subject of the next section.

Diagnostics of coherent resonance and stochastic resonance
To characterize the birhythmic system one can analyze the dynamics of the radius A of the oscillations, as per equation (4). A principal question is: does the quasi-potential Parameters domain for the existence of a single limit cycle (white area) and three limit cycles (shadowed area). The black zone is the domain of asymmetric potential, whereas the grey area is the region of symmetric potential. The nonlinear parameter reads μ = 0.001. Table 1. Amplitudes, frequencies and energy barriers of the limit cycles for the symmetric potential. All data refer to the case μ = 0.001.
approximation give a bona fide bistable potential? Put it in other words, we ask if the system evolution has the same features of the dynamics of a particle subject to an ordinary bistable potential, and if the entrained oscillations give rise to the same type of resonances expected in ordinary potentials. We thus summarize the most common diagnostics of stochastic coherence and resonance. Coherence Resonance (CR) refers to the appearance of coherent oscillations in a dynamical system in the presence of noise. In particular it can refer to the appearance of coherent behaviour at an optimal noise strength, i.e. to the emergence of orderly behaviour in the system due to the noise presence (the FitzHugh -Nagumo model [49] probably being the first prototype model for CR). We propose to extend the same diagnostic used by Pikovski and Kurths [49] to the stochastic birhythmic van der Pol model to detect the occurrence of CR. Also SR is assessed in the literature with several quantities, such as the residence-time distribution density of a particle in one of the potential wells [37], the spectral power amplitude [38,50], the hysteresis loop area (HLA) [51] and the Signal-to-Noise Ratio (SNR) [16], to name few. When adding noise any of the above tools measures an improvement of the signal quality, one speculates that the noise has induced stochastic resonance, as we shall discuss in the following.

Tools to quantify coherent resonance
To quantitatively characterize the regularity of a system in the presence of noise, let us begin with the normalized Auto-Correlation Function (ACF), namely: Noise induced coherence can also be characterized by a decay rate of the ACF, i.e. by a Correlation Time (CT) [22,49], defined as where Var(A(t)) is the variance of the amplitude A(t).
Considering A(t) as a transmitted signal, the correlation time is an indicator of the regularity or the coherence of the signal. The highest value gives the best noise intensity necessary for Coherence Resonance (CR) if, as noise is increased, the CT reaches a maximum. The idea is to compute different correlation times for different noise intensities and find the optimal intensity of the noise D c . Another way to quantify CR is to compute τ cor , to check if there is a noise level above which the random term destroys the coherence. In reference [49] it was also introduced the analysis of the standard deviation of the escape times as a noise-signal ratio, R ET . Coherence resonance is then quantified in terms of the parameter R ET , named Coefficient of Variation (CV): here T is the mean, and T 2 − T 2 is the variance of the ETs from one stable orbit to the other through the unstable orbit of amplitude A 2 . A perfectly ordered system, with escapes occurring at a regular pace, entails R ET = 0. Thus, a minimum of R ET also signals the presence of CR.
The escape times can also be described by means of an effective noise intensity defined as [52,53]: that, at variance with CV, should exhibit a maximum when coherence is enhanced.

Tools to quantify stochastic resonance
The measures introduced so far characterize the intrinsic order of the escape times. In some circumstances, when the system is investigated as a detector, it is interesting to compare how the response of the system is influenced by a periodic (sinusoidal) drive [37,38]. In these cases, one wants to measure the response of the system in the presence of the input E 0 sin ωt compared to the case when the system is solely subject to a random term, that is the Signal to Noise-Ratio (SNR) [54][55][56]. To characterize SR one can numerically calculate SNR using the mean square displacement in the presence of a weak signal and the mean square displacement induced by noise. To do so, one defines the Root Mean Square (RMS) displacement: and then SNR can then be defined as the ratio of the RMS in the presence of a mixture of noise and a signal with the case of pure noise: SR appears to be counterintuitive, inasmuch it seems to imply that the signal quality does not deteriorate as the random noise is increased. In fact, for nonlinear systems with an input signal, only in special circumstances increasing the random noise can actually improve the detection of the corrupted signal [42].

Stochastic like-resonance of birhythmic van der Pol systems
In this section, we discuss the application of the concepts illustrated in Section 3 to the birhythmic van der Pol model discussed in Section 2.1. The purpose is to investigate the effect of noise on the driven birhythmic system governed by equation (1). The two-dimensional system is reduced to a one-dimensional system through the instantaneous amplitude that is the counterpart of the approximation (4). Figures 3a and 3b display the ACF, C(τ ) as per equation (13), for six different values of noise amplitude and for two different sets of parameters α and β, when the quasi-potential is symmetric or asymmetric, respectively (see Tabs. 1 and 2). In these figures the drive term is absent, E 0 = 0, equation (1).
For the set (a) the autocorrelation increases at intermediate values of noise, namely D 10 −4 , see Figure 3a 3 . Similarly, for the set (b), in Figure 3b 2 on observes an increase at D 5 × 10 −5 . It is thus noticeable that a resonance occurs also when E 0 = 0 in equation (1), for the absence of an obvious periodic drive. However, the presence of two periodic attractors clearly gives the possibility of the noise induced escapes and synchronization with the deterministic period proper of the attractors. It is therefore decisive to compare the noise intensity at which a resonance occurs with the time-scale given by the quasi-potential, as in equation (10). For the noise intensity of Figure 3a 2 , D = 10 −4 , equation (10) with the parameter values of Table 1 gives an escape time of about T (1 → 3) ≈ 10 4 . One can conclude that the dominating time scale is given by the activation energy, as derived from the quasipotential (8), and not by the period of the van der Pol oscillations (that is ≈2π, see Tabs. 1 and 2).
Also the CT, given by equation (14), can be used as a measure of the coherence [49]. If noise can cooperate with the deterministic oscillations, we expect that τ cor is maximized at non-zero noise intensity. The dependence of this quantity on the noise amplitude is presented in Figure 4. At variance with reference [49], the CT has a clear minimum at the critical noise amplitude. Panel (a) refers to the symmetrical potential, while panel (b) refers to the asymmetric case. The presence of a minimum is a clear signature of anti-coherence resonance [57]. The minimum value τ cor,min measures the strength of the anticorrelation between the residence time in the two stable states. At the minimum the critical noise amplitude D c reads, for the panel (a) D c = 1.77, while in panel (b) it reads D c = 2.07. An estimate of the escape time through equation (10) is only valid for vanishing noise [25,26], so it cannot be used for D ≈ 1. However, numerical simulations of the full equation (1) give an escape time of about T ≈ 10, i.e. it matches the oscillation period of the attractors. We conclude that CT (measured through τ cor ), is most sensitive to periodic oscillations underlying the van der Pol birhythmic system, oscillations that are neglected by the stochastic averaging. We also conclude that the ACF (C(τ ) in Fig. 3) and the CT (τ cor in Fig. 4) are sensitive to different time scales contained in the system.
In the analysis of the ET it is tempting to employ the oscillation frequencies to have an estimate of the deterministic oscillation periods (the ratio is that in the actual system, at variance with the reduced quasi-potential, something oscillates at the frequency Ω). Using the approximated relationship between the external drive frequency and the escape time for time scale matching SR [37] We therefore seek for SR at noise levels where the average escape is the inverse of the van der Pol oscillation period.  Fig. 3. The autocorrelation function C(τ ) for different noise amplitudes, see equation (13). The set (ai) (top panels) corresponds to the symmetric case, S5, see Table 1, with α = 0.1547, β = 0.006, while the set (bi) (bottom panels) corresponds to the asymmetric case, AS3, see Table 2 16.83, and 21.79. It is evident that at these noise values the noise can cooperate with the intrinsic oscillations of the van der Pol system, giving rise to a significant reduction of the variance of the ET from the inner orbit, and conversely to greatly increase the disorder in the escapes from the other well. Tables 3 and 4 provide the values D at the peaks (dips), and the corresponding escape time T , this shows the decreasing of T 1 and the increasing of T 3 .
The phenomenon is confirmed by the behavior of the effective diffusion D eff , given by equation (16), that can be considered as an alternative measure of SR. As cit  is evident in Figures 5b 1−4 , the effective diffusion closely follows the escape times behavior.
To further investigate the appearance of different timescales in signal enhancement we include the effects of the drive term E 0 at frequency ω in the next section.

Effects of a periodic signal on the birhythmic system
In this section we extend the analytic treatment of the noisy system (1) to the case of an applied periodic drive (3). To do so, we restrict the frequency of the drive to be close to the natural frequency of the van der Pol oscillator.

Analytical treatment of the stochastic driven van der Pol birhythmic system
We consider the driven modified van der Pol system (3), where we assume that the deterministic source can oscillate harmonically close to the natural frequency, i.e. ω = Ω + ν, with ν 1. To analyse equation (3), we introduce two variables y 1 (t), y 2 (t) related to the (x,ẋ) phase-space, that rotate at the frequency of the drive ω: We use the Krylov-Bogoliubov averaging method [58], as in reference [59], with the assumption that friction parameter μ is not too large (μ 1), to obtain the following basic (averaged) equations: where P 1 (y 1 , y 2 ) = νy 2 + μy 1 2 Here, ν ∼ (Ω 2 − ω 2 )/(2ω) and the excitations ξ 1 (t) and ξ 2 (t) are two independent normalized Gaussian white noises with intensity D = D ω 2 . The associated PDFs of y 1 and y 2 are where Here, J 1 and J 2 are the probability currents or stationary state probability currents, for the Fokker Planck equation (22), which are not in general constants. However, for certain conditions on P 1 and P 2 [22], in the steady state J 1 and J 2 vanish. Under these conditions, J 1 and J 2 can be derived from a potential V , as The potential V exists if the following potential condition holds [60]: For our problem one finds where Δ ∼ ν. At this stage, two cases can be considered: case Δ = 0 and case Δ = 0. The optimal situation is the case Δ = 0, in which the harmonic forcing frequency is locked to the frequency of the multi-limit cycle oscillation, i.e. ω ∼ Ω. equation (22) represents an overdamped version of two coupled anharmonic oscillator [61,62]. In the absence of detuning Table 5. Characteristics limit of driving amplitude for bistability.
(Δ = 0) and noise terms the potential of the two coupled anharmonic oscillators is Next we transform the potential (27), into the amplitude and phase variables setting y 1 = A sin(φ), y 2 = A cos(φ). This leads to: In the new variables (A, φ), the parameters μ, α and β are fixed as in Figure 2. For this choice, depending on excitation force amplitude E 0 , the potential has two or four minima and is a two or a four well potential. Stationary oscillations under deterministic excitation are obtained setting: for given values of α and β. The amplitudes are plotted in Figure 6 as a function of the excitation E 0 , for φ (φ = − π 2 and φ = π 2 ), at different frequencies of the external drive. It is observed that the bistabilities zone widens with the increase of E 0 . In Figure 6, the black lines are the stable amplitude combined with both phases, for φ = −π/2 in Figure 6b 1 , and φ = π/2 (in Fig. 6b 2 ). The gray curves represent the unstable amplitude. The limit of the amplitude excitation to permit bistability are summarized in Table 5 for the frequency of excitation ω = 1.
Below the critical value of E 0 of Table 5 the following equations for A and φ characterize the escapes: The stability analysis of equations (30a) and the of the related potential V given by equation (28) show that only the points (A i , φ = −π/2) are stable. In many noisy systems the method leads to a one-dimensional approximation to the response amplitude that significantly simplifies the solution procedure. However, the simplified equation is coupled with the phase equation and of solution of the PDF of amplitude and phase is not possible in general. Thus, to characterize the escape time, we consider equation (30a) and replace the terms containing φ by the corresponding deterministic values. This amounts to set φ = −π/2 in equation (30a), that leads tȯ (31) Therefore, the equation for the amplitude A is independent of φ. The deterministic potential corresponding to equation (31) reads (32) It has two stable states and an unstable state; under the condition resumed in Table 5, namely, it is bistable. The reduced deterministic potential V (A), as a function of amplitude a for different values of the amplitude E 0 , is plotted in Figure 7. One can see that for the outer oscillation amplitude the potential well deepens, the inner potential well becomes shallow as E 0 increases. Through equation (32) one can derive the behavior of the potential barriers ΔV 1 and ΔV 3 , that are shown in Figure 8 as a function of the periodic signal intensity E 0 . It is evident that increasing of E 0 above some limit (reported in Tab. 5), the leftmost potential well, associated to the inner orbit, disappears. Asymmetric forms are still observed according to the value of the control parameters E 0 . This gives rise to the possible switching over the potential well in the presence of noisy excitation. To reveal how the periodic drive affect the escape times process and their distribution of the escape times on right potential well and the left potential well of the system equation (31).
The Fokker-Planck equation corresponding to the reduced equation (31) can be written as Hence, the stationary solution of equation (33) is with N being a normalization constant, and the effective potential V eff (a) is given by To analyze the effect of the deterministic excitation on the escape rates (T 1,3 ) of right limit-cycle A 3 and the left stable limit-cycle A 1 , we consider again the mean exit time  of the system from different well of the effective potential V eff (A). WhenD is small in comparison with the height of the energy barrier [12,18] Through equation (35) one can derive the behavior of the effective potential barriers ΔV effi (a) [12,18]. In fact, if one approximates the full birhythmic behavior with an effective Brownian motion of the particle in a double well potential, the effective potential barriers can be considered as an energy needed to escape from an orbit to the other. The analytic predictions of the energies ΔV eff1 and ΔV eff3 as a function of noise intensity D for different deterministic intensities E 0 are shown in Figure 9.
As in reference [43], due to the shape of Figure 7, the probability distribution is in general very asymmetric; for fixed parameters α, β the probability distribution function p(A) can be localized around a single orbit. For the set of parameters (α, β) in the gray area (the case of an asymmetrical potential) in Figure 1, one notes that a sud-den qualitative change of p(a) in the presence of noise, is denoted as P-bifurcations [43]. This can be seen in Figure 10 that the distribution with two peaks evolves into one peaks as the noise intensity D changes gradually.

Numerical investigation of the response to a sinusoidal drive
In stochastic resonance, besides the analysis of escape times as indicators of coherence, it is important to reconstruct at least the qualitative properties of the system as a detector. In this framework, one interprets the deterministic term in equation (18) as a signal to be revealed analysing the response of the van der Pol oscillating system, and noise is the disturbing term that hinders the detection. A standard tool to characterize the signal enhancement due to noise is the SNR, as per equation (18). To obtain an estimate of such enhancement, we have simulated the van der Pol system at a fixed the value of external drive frequency, ω = 0.2 -a value where the analytic approximation of Section 5.1 is not valid, for it assumes ω Ω. Averaging over 200 realizations for as long as 10 6 normalized time units and in the so-called coherent detection [41,42], i.e. the case where the initial state of the detector (in our case, the vdP oscillator) is fixed and influences the signal detection. The investigation of the interplay between noise and the coherent driving input gives rise to the phenomenon of stochastic resonance. Figure 11 shows the SNR as a function of D in an area of approximately symmetric potential for E 0 = 0.1, ω = 0.2, and μ = 0.001 with different initial condition (black: IC is around the inner orbit A 1 ; gray curve: IC is around the outer orbit A 3 ). In Figure 11a, the SNR shows no peak for both ICs. It is evident that SR does not occur with the corresponding set of parameters α = 0.0675, β = 0.0009. Figures 11b-11f display the response of the system when it appears a resonant-like behavior as a function of the noise level, that is a signature of SR. A stochastic resonancelike phenomenon is only observed when the system starts on the inner orbit with a low value of the parameter of nonlinearity, μ.
In the area where an asymmetric potential is expected (see Fig. 1), the analysis shows some analogous features, as displayed in Figure 12 at E 0 = 0.1, ω = 0.2 and μ = 0.001. In Figure 12a there is no SR-like phenomenon, while in Figures 12b-12d the SNR as a function of D exhibits a maximum for an optimal noise intensity, a characteristic feature of SR phenomena. It is clear that no SR occurs when the system starts around the outer orbit A 3 .

Conclusion
We have investigated the effect of noise added to a forced birhythmic van der Pol-like system that describes some applications as, e.g., enzymatic reactions, sleep-awake cycles, and energy harvesting. The interest also arises from the interplay between the time scales: the two oscillations periods (that are the essence of birhythmicity), the drive frequency (that is externally controlled), and the noise induced average escape time (that is determined by the quasi-potential). Moreover, the system is also analytically convenient. In fact, using the method of stochastic averaging, one can reduce the modified van der Pol equation (3) to an asymmetric bistable system driven by Gaussian white noise.
In the stochastic averaging approximation the two orbits correspond to two fix points, and the stability (in the low noise limit) can be estimated through the quasipotential, also in the presence of a forcing term at a frequency close to the frequency of self-oscillations. The numerical analysis reveals several features of the driven birhythmic system: (a) The ACF shows that the dominating time scale is dictated by the activation energy (as derived by the quasi-potential) and not by the spontaneous oscillations periods, Figure 3. (b) There is a special value for the noise D c at which the ACF reaches a minimum, a signature of anticorrelation. The D c values are slightly different for the inner and outer orbit, Figure 4. (c) The escape times distribution (as measured by the coefficient of variation CV or equivalently by the effective diffusion constant D eff ) depends upon the self generated oscillations of the birhythmic van der Pol, Figure 5. One concludes that the different tools to characterize the system are sensitive to different phenomena, ACF to (rare) escapes from the stable orbits, the escape times to the oscillations of the attractors. Thus, the escapes from the quasi-potential behave as expected for an ordinary potential, but they also keep memory of the fact that the system is in fact oscillating (the stable point in the averaged system represents an orbit) and can (stochastically) resonate also at the orbit frequency.
The analytical treatment (by means of stochastic averaging) allows to estimate the behavior of the quasipotential (or the effective activation energies, V in Fig. 7) as a function of the drive amplitude (Fig. 8) and of the white noise intensity (Fig. 9). The analysis reveals that the drive amplitude can either increase or decrease the stability of the orbits. Instead, the noise amplitude always tends to stabilize the outer orbit (i.e., to increase the effective activation energy) and to weaken the stability of the inner orbit (i.e., to decrease the effective activation energy). This type of analysis can only be performed for a drive frequency that is close to the oscillation frequency. As the stochastic averaging is only applicable when the drive frequency is equal to the oscillator frequency, we have employed numerical simulations to evaluate the SNR for slow frequency (respect to the natural frequency of the van der Pol oscillator), in the so called coherent detection approach, that keeps memory of the initial conditions. The SNR associated to a small periodic drive shows markedly different properties if the initial conditions are selected on the inner or outer orbit: depending on the parameter set, stochastic resonance is observed either for the inner or for the outer orbits. This is compatible with a picture of bistable oscillations between the two stable orbits. However, it remains to be established if the incoherent detection case, where the average washes out the influence of the initial conditions, also exhibits stochastic resonance.
R.Y. undertook this work with the support of the Max Planck Society, Germany. He acknowledges the support of the Max Planck-Institut of Mathematical Sciences, Leipzig, Germany.

Author contribution statement
This paper was proposed by Prof R. Yamapi and Prof. G. Filatralla. C. Chéagé Chamgoué is the main author of the paper, while Prof. P. Woafo was the supervisor of this works.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.