Control of birhythmicity: A self-feedback approach

TL;DR

This paper introduces a self-feedback control method to eliminate birhythmicity in diverse oscillatory systems, validated through analytical and numerical analysis across biological and engineering models.

Contribution

A novel self-feedback control scheme is proposed and validated for controlling birhythmicity in various natural and artificial systems.

Findings

Control scheme effectively eliminates birhythmicity.

Can induce transitions between bistable states.

Applicable to biochemical and engineering systems.

Abstract

Birhythmicity occurs in many natural and artificial systems. In this paper we propose a self-feedback scheme to control birhythmicity. To establish the efficacy and generality of the proposed control scheme, we apply it on three birhythmic oscillators from diverse fields of natural science, namely, an energy harvesting system, the p53-Mdm2 network for protein genesis (the OAK model) and a glycolysis model (modified Decroly-Goldbeter model). Using the harmonic decomposition technique and energy balance method we derive the analytical conditions for the control of birhythmicity. A detailed numerical bifurcation analysis in the parameter space establishes that the control scheme is capable of eliminating birhythmicity and it can also induce transitions between different forms of bistability. As the proposed control scheme is quite general, it can be applied for control of several real systems, particularly in biochemical and engineering systems.

Equations

$$\displaystyle\frac{d^{2}y}{dt^{2}}+\mu F\bigg{(}\frac{dy}{dt}\bigg{)}+\Omega_{0}^{2}y$$

$$\frac{dv}{dt}+\gamma v$$

$$F\bigg{(}\frac{dy}{dt}\bigg{)}=-\frac{dy}{dt}+\frac{1}{3}\bigg{(}\frac{dy}{dt}\bigg{)}^{3}-\frac{\alpha}{5}\bigg{(}\frac{dy}{dt}\bigg{)}^{5}+\frac{\beta}{7}\bigg{(}\frac{dy}{dt}\bigg{)}^{7}.$$

$$\ddot{x}+\mu F(\dot{x})+x=G(x,\dot{x},t).$$

$$\displaystyle\frac{d^{2}y}{dt^{2}}+\mu F\bigg{(}\frac{dy}{dt}\bigg{)}+\Omega_{0}^{2}y$$

$$\frac{dv}{dt}+\gamma v$$

$$\ddot{y}+\mu F(\dot{y})+{\Omega }_0^{2}y+d\dot{y}+\frac{{\eta }_0}{\gamma }({\eta }_1\dot{y}+\dot{v})=0.$$

$$\ddot{y}+\mu F(\dot{y})+{\Omega }_0^{2}y+d\dot{y}=0.$$

$$y(t)=A\cos (\omega t),$$

$$\begin{split}(&-\omega^{2}+\Omega_{0}^{2})A\cos(\omega t)\\ =&\mu\omega\bigg{(}-1+\frac{1}{4}\omega^{2}A^{2}-\frac{\alpha}{8}\omega^{4}A^{4}+\frac{5\beta}{64}\omega^{6}A^{6}\bigg{)}A\sin(\omega t)\\ &+d\omega A\sin(\omega t)\\ &+\mu\omega^{3}\bigg{(}-\frac{1}{12}+\frac{\alpha}{16}\omega^{2}A^{2}-\frac{3\beta}{64}\omega^{4}A^{4}\bigg{)}A^{3}\sin(3\omega t)\\ &+\mu\omega^{5}\bigg{(}-\frac{\alpha}{80}+\frac{\beta}{64}\omega^{2}A^{2}\bigg{)}A^{5}\sin(5\omega t)\\ &-\frac{\beta}{448}\mu\omega^{7}A^{7}\sin(7\omega t).\end{split}$$

$$\begin{split}(&-\omega^{2}+\Omega_{0}^{2})A\cos(\omega t)\\ =&\mu\omega\bigg{(}-1+\frac{1}{4}\omega^{2}A^{2}-\frac{\alpha}{8}\omega^{4}A^{4}+\frac{5\beta}{64}\omega^{6}A^{6}\bigg{)}A\sin(\omega t)\\ &+d\omega A\sin(\omega t)+\mathcal{H},\end{split}$$

$${\omega }^2-{\Omega }_0^2=0,$$

$$\mu\bigg{(}1-\frac{1}{4}\omega^{2}A^{2}+\frac{\alpha}{8}\omega^{4}A^{4}-\frac{5\beta}{64}\omega^{6}A^{6}\bigg{)}-d=0.$$

$$\mu\bigg{(}1-\frac{1}{4}\Omega_{0}^{2}A^{2}+\frac{\alpha}{8}\Omega_{0}^{4}A^{4}-\frac{5\beta}{64}\Omega_{0}^{6}A^{6}\bigg{)}-d=0.$$

$$y(t)=A\cos (t+\varphi ),$$

$$y(t)\backsimeq A\cos (t),$$

$$\Delta E$$

$$f(A^{2})=\mu\bigg{(}1-\frac{1}{4}A^{2}+\frac{\alpha}{8}A^{4}-\frac{5\beta}{64}A^{6}\bigg{)}-d=0.$$

$$\frac{d}{dA}\bigg{(}\Delta E(A)\bigg{)}\bigg{\rvert}_{\mbox{limit cycle}}<0.$$

$$\lambda_{1,2}=\frac{1}{2}\bigg{[}(\mu-d)\pm\sqrt{(d-\mu)^{2}-4\Omega_{0}}\bigg{]},$$

$$d_{\text{mbox}HB}=\mu ,$$

$$\frac{dP}{dt}$$

$$\frac{dMc}{dt}$$

$$\frac{dMn}{dt}$$

$$\frac{d\alpha }{dt}$$

$$\frac{d\gamma }{dt}$$

$$\varphi (\alpha ,\gamma )=\frac{\alpha (1+\alpha )(1+\gamma {)}^2}{L+(1+\alpha {)}^2(1+\gamma {)}^2}.$$

$$\ddot{y}+\mu F(\dot{y})+y+d\dot{y}=0.$$

$$y(\tau )=y_0(\tau )+\mu y_1(\tau )+{\mu }^2{y}_2(\tau )+\cdots ,$$

$$\omega =1+\mu {\omega }_1+{\mu }^2{\omega }_2+\cdots ,$$

$$d=\mu d_1+{\mu }^2{d}_2+\cdots$$

$${\ddot{y}}_0+y_0=0.$$

$${\ddot{y}}_1+y_1=-2{\omega }_1{\ddot{y}}_0-F({\dot{y}}_0)-d_1{y}_0.$$

Full text

Multistability appears in diverse forms and their study is an exciting topic of research in science and engineering. A particular form of multistability is bistability: it shows many variant, such as, the coexistance of two stable steady states, one stable steady state and one stable limit cycle, two stable limit cycles, or two chaotic attractors. Birhythmicity is the phenomenon of coexistence of two stable limit cycles separated by an unstable limit cycle with different amplitudes and frequencies. In many physical systems birhythmicity is undesirable as in energy harvesting systems but in most biological systems, e.g., enzymatic oscillations, it is desirable. Therefore, control of birhythmicity is of utmost importance. Although the control of multistability is a well studied topic, the control of birhythmicity has not been explored to that extent. In this paper we propose a control scheme that can effectively control and, whenever required, can eliminate birhythmicity. We theoretically explore and numerically establish the technique of control of birhythmicity and transitions to any desired attractor. A number of engineering and biological systems are investigated with the proposed control scheme to establish the efficacy and generality of the scheme. The main essence of this control scheme lies in the fact that it is easily realizable and offers an efficient mean to control birhythmicity.

I Introduction

Birhythmicity is a variant of multistability Pisarchik and Feudel (2014), which arises in many natural and artificial systems in the field of physicsKwuimy and Nataraj (2015), biologyGoldbeter (1996, 2002); Kar and Ray (2004) and chemistryMohamed and R (1983). The coexistence of two stable limit cycles of different amplitudes (and frequencies) separated by an unstable limit cycle is the signature of birhythmicity. Birhythmicity may appear in chaotic oscillators also, e.g., the coexistence of two chaotic attractors has been reported in the literatureLozi and Ushiki (1991) and studied in detail in Ref. Pisarchik et al., 2006. The appearance of birhythmicity plays a crucial role in living systems as it helps to maintain different modes of oscillations that organize various biochemical processes in response to variations in their environment Mohamed and R (1983). That is why most of the biochemical oscillators are birhythmic. A few prominent examples are: glycolytic oscillator and enzymatic reactions Goldbeter (1996, 2002); Kar and Ray (2003, 2004), intracellular $\mbox{Ca}^{2+}$ oscillations Pisarchik et al. (2006), birhythmic oscillations due to the complex regulatory properties of allosteric enzymes, namely phosphofructokinase (PFK), which is activated by ADP and feedback due to this ADP to ATPKar and Ray (2004), birhytmicity in the p53-Mdm2 networkAbou-Jaoudé, Ouattara, and Kaufman (2009); Vogelstein, Lane, and Levine (2000), oscillatory generation of cyclic AMP during the aggregation of slime mold Dictyostelium discodeumMartiel and Goldbeter (1987) or circadian oscillation in PER and TIM proteins in DrosophilaLeloup and Goldbeter (1999). Unlike in living systems, birhythmicity is often undesirable in physical Yamapi and Filatrella (2014) and engineering systems Kwuimy and Nataraj (2015). For example, birhythmicity limits the efficiency of an energy harvesting system Kwuimy and Nataraj (2015): In an energy harvesting system the wind-induced vibration shows birhythmic oscillations, thus, depending upon the initial vibrational energy of the wind the system may oscillate in a small amplitude limit cycle and results in less mechanical deformation that, in turn, yields less electrical energy. Therefore, for an energy harvesting system the oscillation with a large amplitude limit cycle is always desirable to have larger energy production. Clearly in some situations birhythmicity is undesirable while in others it is a necessity. This marks the importance of control of birhythmicity.

A recent extensive review work on control of multistability by Pisarchik and Feudel (2014) suggests that although several control mechanisms have been reported for the control of bistable systems containing oscillations and stable steady statePisarchik and Goswami (2000); Pisarchik (2001), the control of birhythmicity is a less explored topic. Only a few works are reported on the control of bithythmicity. Ghosh et al. (2011) showed that time delay feedback control, which was originally proposed by Pyragas (1992) to control chaos, is able to control birhythmicity as well in a modified birhythmic van der Pol oscillator. However, owing to the presence of time delay in their control scheme a detailed bifurcation analysis for the controlled system is a difficult task. Also, the implementation of a delayed signal to be fed is challenging. Sevilla-Escoboza et al. (2015) showed that application of a harmonic modulation and the presence of a positive feedback along with a proper choice of the parameters can transform a multistable system with coexisting periodic and chaotic attractors to a monostable one. Recently, in Ref. Biswas, Banerjee, and Kurths, 2016 we have proposed a technique to control birhythmicity by using a conjugate self-feedback method. This technique was verified using a variant of the van der Pol oscillator with birhythmic oscillations and it was shown that the conjugate self-feedback in that oscillator is capable of removing birhtymicity by inducing monorhythmic oscillation. However, in our scheme in Ref. Biswas, Banerjee, and Kurths, 2016 one requires the access of two variables: the variable of interest and its conjugate counterpart.

In this paper we propose a more general and experimentally feasible control technique that employs only one accessible variable. We establish the effectiveness of this control technique using three real systems from diverse field of physics and biology, namely (i) an energy harvesting system Kwuimy and Nataraj (2015), (ii) the p53-Mdm2 network popularly known as the OAK model Abou-Jaoudé, Ouattara, and Kaufman (2009) and (iii) a variant of glycolytic oscillators Goldbeter (1996); Kar and Ray (2004). The control and taming of birhythmic oscillations in these three oscillators of different origin also establish the generality of our control scheme. To establish the efficacy of our scheme, we carry out an extensive theoretical analysis using the harmonic decomposition technique and the energy balance method. Also, we employ a rigorous numerical bifurcation analysis to identify the parametric zone of occurrence of bi- and mono-rhythmic oscillations and their exact genesis.

The paper is organized in the following manner: The next section describes the details of control of birhythmicity in an energy harvesting model. We carry out a detailed analysis for the onset of birhythmicity. In sections III and IV we consider the control of birhythmicity in the p53-Mdm2 network (OAK model) and glycolytic oscillator, respectively. Finally, section V concludes the outcomes of the study.

II Energy harvesting system

II.1 The model

At first we describe the original model of an energy harvesting system. Energy harvesting systems generate electrical energy from ambient energy arising from sources like structural vibration, wind flow, physiological and chemical reactions, etc. Abdelkefi, Yan, and Hajj (2013); *Kwon10; *Kwuimy14pla; *Kwuimy12; *Litak12; *Tekam14; *Tekam15. Kwuimy and Nataraj (2015) considered an energy harvesting model which is implemented with an arrangement of a cantilever attached to piezoelectric patches under the action of transverse wind flow. The physical model consists of an electrical circuit with a load resistance and a flexible beam of distributed piezoelectric patches. The dimensionless form of the original model is given by the following set of equations:

[TABLE]

Here $y$ is the dimensionless transversal beam deflection function and $v$ is the dimensionless form of the voltage generated by the piezoelectric element. Also, $\mu$, $\Omega_{0}$, $\gamma$, $\eta_{0}$ and $\eta_{1}$ are all positive parameters and the nonlinear function F\bigg{(}\frac{dy}{dt}\bigg{)} is given by

[TABLE]

The system is birhythmic for the following parametersKwuimy and Nataraj (2015): $\mu=0.1$, $\alpha=0.144$, $\beta=0.005$, $\Omega_{0}=1$, $\eta_{0}=0.1$, $\eta_{1}=0.25$ and $\gamma=0.2$. It shows three distinct limit cycles (LCs) (two stable and one unstable) depending upon two sets of initial conditions (IC), namely, \mathcal{I}_{1}\equiv\big{(}y(0),dy(0)/dt,v(0)\big{)}\in(0.1,0,0.3) (small amplitude stable LC) and \mathcal{I}_{2}\equiv\big{(}y(0),dy(0)/dt,v(0)\big{)}\in(7,0,0.3) (large amplitude stable LC). The unstable LC determines the basin boundary of these stable LCs. The time series and phase plane plots for these sets of initial conditions are shown in Fig. 1(a) and Fig. 1(b), respectively. From Fig. 1(b) we see that the system shows two limit cycles separated by an unstable LC indicating birhythmicity.

II.2 The self-feedback control scheme

We consider a dynamical system modeled by the following equation:

[TABLE]

Where $x\in\mathbb{R}$, $F(\dot{x})$ is a nonlinear function, and $\mu$ is a system parameter that determines the intrinsic dynamics of the considered system. Here $G(x,\dot{x},t)$ is the control term. In our present case it is given by $G(x,\dot{x},t)=-d\dot{x}$, where $d$ controls the strength of the self-feedback term. Note that to implement the system one needs the access to a single (scalar) variable. This type of self-feedback has been used and implemented earlier in the Fabry-Perot laser diode system Su et al. (2014), biological system to control complex motor task Sigrist et al. (2011) and learning Chiviacowsky and Wulf (2002). Another variant of the self-feedback control is the paradigmatic Pyragas control technique Pyragas (1992) where an additional time-delayed version of the variable is used. In this paper our goal is to study the effect of self-feedback term on birhythmicity; we concentrate on the fact that how the self-feedback parameter $d$ affects the birhythmic oscillation. Since birhythmicity involves global bifurcations therefore we use continuation based rigorous bifurcation analysis along with theoretical analysis to track the complete system behavior.

II.3 Control of the energy harvesting system with self-feedback

Next we apply the self-feedback control scheme to (1), which now reads

[TABLE]

The term $-d\frac{dy}{dt}$ in (4a) represents the self-feedback proportional to the time rate of change of the transversal beam deflection with $d$ as the self-feedback strength. $d$ also determines the nature of the self-feedback: $d>0$ represents a positive feedback, and $d<0$ represents a negative one and $d=0$ implies no feedback.

In order to analyze the controlled system, we reduce Eq. (4) to the following single equation

[TABLE]

According to Refs. Kwuimy and Nataraj, 2015; Kwuimy and Litak, 2014, the variable $v$ does not affect the dynamical property of the fixed point of the system and hence they showed that one can put $\eta_{0}=\eta_{1}=0$. Then Eq. (5) becomes

[TABLE]

Next, we will constitute the amplitude equation of Eq. 6 to predict the kind of bifurcation structures associated with the system. Generally, for this, the most common technique in the literature is to apply the harmonic decomposition technique (see for example Refs.Yamapi, Nbendjo, and Kadji, 2007; Ghosh et al., 2011; Biswas, Banerjee, and Kurths, 2016). Although the harmonic decomposition technique is not an asymptotic method, however it can predict the amplitude equation in a simple yet effective way. In this context one more technique, which is a much more suitable technique for weakly nonlinear systems, is the Poincaré-Lindstedt technique which is discussed in Appendix-A; it is shown that both the analyses give equivalent amplitude equations and match well with the numerical results.

According to the harmonic decomposition technique we assume the approximate solution of Eq. (6) as

[TABLE]

where $A$ is the amplitude and $\omega$ is the frequency of the oscillator with feedback. Substitution of this in Eq. (6) yields

[TABLE]

According to Ref. Jordan and Smith, 1999 we can treat the higher harmonic terms as forcing terms and ignore them. Thus Eq. (8) reduces to

[TABLE]

where $\mathcal{H}$ denotes higher harmonic terms.

From Eq. (9) we get the following frequency and amplitude equations:

[TABLE]

and

[TABLE]

From Eq. (11) we infer that it is not a purely amplitude equation as it contains the frequency $\omega$ in it. So we substitute the value of $\omega$ from Eq. (10) in (11) and get

[TABLE]

It may be noted that for $d=0$ Eq. (12) represents the amplitude equation of the uncontrolled energy harvesting system. It is interesting to note that the amplitude of the system does not depend on $\mu$ unless $d\neq 0$. The frequency in the harmonic limit becomes $\omega=1$. Also the frequency equation viz. Eq. (10) states that the frequency of the system does not depend on the feedback strength $d$ (unlike Ref. Biswas, Banerjee, and Kurths, 2016), thus, leaving the original frequency of the system intact. The three roots of Eq. (12) corresponds to the amplitudes of the three limit cycles (two stable and one unstable).

To test the stability of the system, we apply the energy balance method as suggested in Ref Ghosh et al., 2011. For $\mu=0$ and $d=0$, the harmonic solution of Eq. (6) may be given as Yamapi, Nbendjo, and Kadji (2007)

[TABLE]

where $\phi$ is the initial phase and may be considered $\phi=0$ for convenience. The phase plane of this is a circle with period $T=2\pi$. In the presence of a self-feedback we can approximate

[TABLE]

Let us consider \big{(}-\mu F(\dot{y})-d(\dot{y})\big{)} as the external forcing term to calculate the change in energy $\Delta E$ during one period, i.e., $0\leq t\leq T$, with $T=2\pi$. The change in energy is given by

[TABLE]

$\Delta E=0$ for a periodic solution (limit cycle). Hence from the above integral we get using the condition of Eq. (14)

[TABLE]

Note that, Eq. (16) is equivalent to Eq. (11) for $\omega=1$. The number of limit cycles can be obtained by solving Eq. (16) by normalizing the frequency to unity. The number of positive roots gives the number of LCs. The stability of the limit cycle is determined by the slope of the curve of Eq. (16) at the zero crossing point. The negativity of the slope determines the stability of the LC. The condition of stable limit cycle thus can be written as

[TABLE]

We solve the amplitude equation (16) (which is the same as Eq. (12) for $\Omega_{0}=1$) by graphical method, i.e., we plot the polynomial $f(A^{2})$ with $A^{2}$. The zero crossing points of $f(A^{2})$ are the solutions of the equation. Here we consider the following parameter values: $\mu=0.1$, $\alpha=0.144$, $\beta=0.005$. These are the values of the parameters for which the original system (4) exhibits birhythmicity in the absence of self-feedback. Now we vary the control parameter $d$ to get different solutions. The number of solutions of Eq. (16) determines the number of limit cycles. The solutions for different values of the control parameter $d$ is shown in Fig. 2. We vary $d$ from positive high to negative high (curves in Fig. 2 from lower o upper). From the figure we see that there is no zero crossing for $d=0.2$ (yellow line with hollow circles). Thus there is no limit cycle for this value of the control parameter. In other words at $d=0.2$ the system is in a stable steady state (SSS). Decrease of $d$ causes $f(A^{2})$ to shift upwards and eventually to cross the zero line, thus giving rise to a stable limit cycle. This case is shown for $d=0.05$ in the figure (green line with solid circles). At this value of the parameter there is only one stable limit cycle with low amplitude. The negative slope of the curve at the zero crossing indicates that the the LC is stable. Next, at $d=0$ (solid line) there are three zero crossing points. Two of them have negative slopes and one has positive slope at the zero crossing points. Thus, there exists birhythmicity with two stable LCs (with negative slopes) and one unstable LC (with positive slopes) in between them. At $d=-0.2$ (purple line with hollow squares) there is only one zero crossing of $f(A^{2})$ with a negative slope at the zero crossing exhibiting the presence of only one stable LC with larger amplitude.

In the absence of control the system undergoes only a global bifurcation, namely, a saddle-node bifurcation of limit cycle (SNLC) and a codimension-2 cusp bifurcation in the $\alpha-\beta$ parameter space. However, the presence of the control term causes the system to experience a local bifurcation, namely, a Hopf bifurcation (HB). The eigenvalues of the Jacobian of Eq. (6) around the stable points $(y,\dot{y})=(0,0)$ are given by

[TABLE]

and the condition of Hopf bifurcation reads

[TABLE]

where $d_{\mbox{HB}}$ is the critical value of $d$ for a Hopf bifurcation to occur.

To investigate the detailed bifurcation scenario in the system we use the continuation package XPPAUTErmentrout (2002) in the $d-\mu$ parameter space. The two-parameter bifurcation diagram in the $d-\mu$ plane is shown in Fig. 3(a).

The parameter values are $\alpha=0.144$, $\beta=0.005$, $\Omega_{0}=1$, $\eta_{0}=0.1$, $\eta_{1}=0.25$ and $\gamma=0.2$: this set of parameter values confirm birhythmic oscillations in the uncontrolled system. From Fig. 3(a) we observe that the $d-\mu$ plane is divided by SNLC curves and a HB curve. The birhythmic regime exists between two SNLC curves [cyan (lighter gray) zone]. In this zone three LCs coexist, one with larger amplitude, another with lower amplitude and the third one is unstable. The SNLC curves also govern the transition between the birhythmic zone to the monorhythmic one. The monorhythmic zones are shown by the yellow (light gray) region. The HB line governs the transition between LC and the stable steady state (SSS) [blue (dark) zone]. The HB line exactly matches with the analytically obtained result of Eq. (19). It may be noted that for $d=0$ (i.e., the uncontrolled case) the system is in a birhythmic regime for all values of $\mu$ for the chosen set of other parameter values.

Next, we take $\mu=0.1$ (horizontal broken line in Fig. 3(a)) and vary $d$. The resulting one-parameter bifurcation diagram is shown in Fig. 3(b). The system is in the birhythmic zone for $d=0$. An increase in $d$ brings the system to a monorhythmic zone through an SLNC for $d>d_{U}$. Here the limit cycle is a small-amplitude limit cycle. Further increase in $d$ causes the small amplitude LC to loose stability through an inverse Hopf bifurcation and the system rests in a stable steady state. Therefore, increasing the positive value of the control parameter brings the birhythmic oscillator to a monorhythmic one with small-amplitude oscillation and eventually to a stable steady state. On the other hand, decreasing value of $d$ in negative direction causes the system to experience an SNLC at $d=-d_{L}$ and monorhythmic oscillation with larger amplitude emerges. To summarize, the proper choice of the control parameter $d$ may cause the system to induce monorhythmicity either with a small amplitude LC (for $d_{U}<d<d_{HB}$) or a large amplitude LC (for $d<-d_{L}$). It is to be noted that, there is a hysteresis around $d=0$ with a width of $\Delta d=d_{U}-d_{L}$, which is shown by the cyan (lighter grey) zone in Fig. 3(b). In this zone the state of the system depends on the choice of the initial conditions. It may be noted that the hysteresis zone increases with increasing the value of $\mu$. We plot the hysteresis width $\Delta d$ with $\mu$ in Fig. 4 for the aforementioned parameters, which shows that the hysteresis width increases with an increase in $\mu$. Therefore, to eliminate the birhythmicity a higher feedback strength ($d$) is necessary for higher values of $\mu$.

To demonstrate the effectiveness of the control scheme, we plot the time series and the phase plane diagrams for different values of the control parameter $d$ (Fig. 5). Here we use two sets of initial conditions, namely, $\mathcal{I}_{1}=(0.1,0,0.3)$, which is around the origin, targeting to have the smaller amplitude LC and the other $\mathcal{I}_{2}=(7,0,0.3)$, away from origin, targeting the larger amplitude LC.

The solid line indicates the result for $\mathcal{I}_{1}$ and the line with hollow circles that for $\mathcal{I}_{2}$. Figure 5(a) shows time series and Fig. 5(b) shows the phase plane plot for $d=-0.2$ (i.e., $d<-d_{L}$). We observe that irrespective of the initial conditions the system always shows a large amplitude LC. Figure 5(c-d) demonstrate the scenario for $d=-0.01$, i.e., $-d_{L}<d<d_{U}$, for which the system is in the birhythmic region. The occurrence of LCs with two different amplitudes confirms the presence of birhythmicity in the system. For $d>d_{U}$ there exists only small amplitude LC. This is shown for $d=0.05$ in Fig. 5(e-f). Finally, the oscillation in the system ceases to a stable fixed point for $d>d_{HB}$: This is demonstrated for $d=0.2$ in Fig. 5(g-h).

III The p53-Mdm2 network: OAK model

The control of proliferation of abnormal cells by protein in mammals is modeled through the p53-Mdm2 network, which is called the OAK model originally proposed by by Abou-Jaoudé et al.Abou-Jaoudé, Ouattara, and Kaufman (2009) The OAK model describes the interaction between p53, cytoplasmic Mdm2 and nuclear Mdm2 Ciliberto, Novak, and Tyson (2005). A detailed description of the model may be found in Ref. Abou-Jaoudé, Ouattara, and Kaufman, 2009; W. A-Jaoudé, 2011. Nuclear Mdm2 accelerates the degradation of p53 by ubiquitination and by blocking its functional activity. p53 enhances the transcription of gene MDM2 and thus regulates cytoplasmic Mdm2 level. The translocation of Mdm2 from the cytoplasm to the nucleus is inhibited by p53. Although the actual model consists of a 4-dimensional differential equations But in Ref. W. A-Jaoudé, 2011 the model has been reduced to a 3-dimensional one. We apply the self-feedback control term $-dMc$ and rewrite the equations as follows

[TABLE]

where $P$, $Mc$ and $Mn$ represent the concentrations of p53, cytoplasmic Mdm2 and nuclear Mdm2, respectively. Here the term $dMc$ in the second equation is the proposed self-feedback. The schematic of the control scheme is represented in Fig. 7.

The self-feedback occurs in $Mc$: a portion of $Mc$ is feed to itself (red line). The negative $d$ invokes the positive interaction and is shown by the sharp arrowhead in Fig. 7(I). The positive $d$ represents the negative interaction, i.e., suppression of $Mc$ which is shown in Fig. 7(II).

For $d=0$ we get the original system which shows birhythmicity for the parameter set: $k_{p}=5$, $K_{p}=0.2$, $d_{p}=2.5$, $k_{Mc}=0.1$, $k_{Mc}^{\prime}=1.2$, $K_{Mc}=0.4$, $k_{in}=0.45$, $k_{in}^{\prime}=0.4$, $K_{Mn}=0.1$, $d_{Mc}=0.6$, $V_{r}=10$, $d_{Mn}=1.9$ and $n=6$. The small amplitude LC arising from this system has the initial condition $\mathcal{I}_{1}\equiv(P(0),Mc(0),Mn(0))\in(0.6,0.3,0.4)$ and the large amplitude LC is a consequence of the initial condition $\mathcal{I}_{2}\equiv(P(0),Mc(0),Mn(0))\in(3,0.3,0.2)$. The time series and phase plane plots for the original system (i.e., $d=0$ in Eq.(III)) is shown in Fig. 6(a) and (b), respectively. To control the birhythmicity, we have applied the self-feedback scheme by using $d\neq 0$. We analyze the bifurcation scenario appearing in this system. The two-parameter bifurcation diagram in the $d-d_{Mn}$ parameter space is shown in Fig. 8(a). The yellow (light gray) zone shows the monorhythmic zone. The cyan (lighter gray) zone gives the birhythmic zone and the blue (dark) zone is the zone of a stable steady state (SSS).

In this case also we find that the space is divided by SNLC bifurcations and a supercritical Hopf bifurcation. The global SNLC bifurcation distinguishes between the birhythmic and monorhythmic states. The stable limit cycle looses stability through an inverse Hopf bifurcation and a stable steady state emerges.

To understand the scenario in more detail, we draw the one parameter bifurcation diagram sweeping the control parameter $d$ by fixing $d_{Mn}=1.9$ as shown in Fig. 8(b). The variation of $d$ is considered along the broken yellow horizontal line in Fig. 8(a). For $d=0$, the system is in the birhythmic zone for the preferred set of parameter values. Increase in $d$ brings the system to monorhythmic oscillations through a global SNLC for $d>d_{U}$. This monorhythmic oscillation is of small amplitude. Further increase in $d$ causes the stable LC to loose its stability and a stable steady state emerges through a Hopf bifurcation at $d=d_{HB}$. Decrease in $d$ below zero causes the system to enter the monorhythmic zone again but with large amplitude oscillation for $d<-d_{L}$. From the one-parameter bifurcation diagram we see that there is a hysteresis zone of width $\Delta d=(d_{U}-d_{L})$, which is governed by the SNLC curves.

Finally, we draw the time series and phase plane plots of the system for different values of the control parameter $d$ in Fig. 9. The solid line shows the case of initial condition $\mathcal{I}_{1}$, and the line with hollow circles shows the same for $\mathcal{I}_{2}$, respectively.

Fig. 9(a) gives the time series and (b) shows the phase plane plots for $d=-0.35$. Both LCs are of the same large amplitude irrespective of initial conditions indicating monorhythmic oscillation. For $d=0.04$, the system is in the bithythmic zone shown in Fig. 9(c-d). Increasing $d$ brings the system to the monorhythmic zone with small amplitude oscillation for $d>d_{U}$. The situation for $d=0.1$ is shown in Fig. 9(e-f). Finally the system enters the stable steady state for $d>D_{HB}$: Fig. 9(g-h) demonstrates the scenario for $d=0.35$. Thus, it is worth noting that the birhythmicity may be controlled and eliminated with proper choice of the control parameter $d$.

IV The Glycolysis model: modified Decroly-Goldbeter model

Enzymatic oscillations with periodicity (of several minutes) have two-fold interest in biology, first, in metabolic pathways and second, as general models for biological rhythms. These examples include glycolytic oscillations in yeast and muscle and the periodic synthesis of cAMP during the aggregation of the slime mold Dictyostelium discoideum Martiel and Goldbeter (1987); Gorbunova and Spitzer (2002). Kar and Ray (2004) consider a product-activated enzyme model which is a modified version of the well known Decroly-Goldbeter modelGoldbeter (1996). According to them, the allosteric enzymes consist of multiple identical subunits. These subunits undergo conformational transition between more reactive (R) and less reactive (T) states. Here the substrate ($(S)$) injection rate $\nu$ is constant. The product $P$ is resulted from the bindings of $S$ with $R$ and $T$ states of the enzyme. Then the product is removed with a rate proportional to its concentration, which results in a positive feedback and activates the $T$ to $R$ transitions. When the product $P$ gives a positive feedback to substrate $S$ birhythmicity results. The system dynamics is described by the following equations along with the self-feedback mechanism

[TABLE]

with,

[TABLE]

Here $\alpha$ is the normalized substrate concentration and $\gamma$ is the normalized product concentration. The self-feedback control is provided by the $-d\gamma$ term in the second subequation of (21). The schematic of the pathways is shown in Fig. 10 Goldbeter (1996); Kar and Ray (2004); the self-feedback term is shown using red lines in Fig. 10. According to Ref. Goldbeter, 1996 the reaction product ($P$) leaves the system at a rate proportional to its concentration. Therefore, in our case a positive self-feedback represented by Fig. 10 (I) means an accumulation of the reaction product ($P$), whereas, a negative feedback represented by Fig. 10 (II) means an extraction of the reaction product ($P$). Here we make use of the kinetic assumption that for a step catalysed by a Michaelian enzyme is not saturated by its substrate Goldbeter (1996). Earlier, theoretical Goldbeter (1996); Kar and Ray (2004) and experimental observations Danø, Sørensen, and Hynne (1999) of the glycolytic oscillations were made by changing the substrate injection rate $\nu$. However, in the present case we keep $\nu$ constant and study the effect of the control on the reaction product ($P$).

For $d=0$ the system is reduced to the original one and is capable of showing birhythmicity for the following set of parameters: $\nu=0.255$, $q=1.0$, $K_{s}=0.06$, $L=3.6\times 10^{6}$, $\sigma=10$, $\sigma_{i}=1.3$, $n=4$ and $K=10.0$. The system exhibits large amplitude LC for the initial conditions $\mathcal{I}_{1}\equiv(\alpha(0),\gamma(0))\in(100,5)$ and small amplitude LC for $\mathcal{I}_{2}\equiv(\alpha(0),\gamma(0))\in(80,5)$, respectively. The time series and phase plane plots for the original system (i.e., $d=0$ in Eq.(21)) are shown in Fig. 11(a) and (b), respectively.

The control is active for $d\neq 0$. We investigate the possible bifurcation scenario appearing in the system. In Fig. 12(a) we present the two-parameter bifurcation diagram in the $d-\sigma_{i}$ parameter space. It is divided in birhythmic, monorhythmic and stable steady state zones by SNLC bifurcation and Hopf bifurcation curves. The SNLC bifurcation governs the transition between birhythmicity and monorhythmicity and the (inverse) Hopf bifurcation brings the system to a stable steady state. The yellow (light gray) zone in Fig. 12(a) shows the monorhythmic oscillatory zone. The cyan (lighter gray) zone shows the birhythmic regime and the blue (dark) zone the stable steady state. Here another interesting bistable state exists, namely, the coexistence of one large amplitude LC and one stable steady state between the Hopf and SNLC curves (the purple zone).

To have a better understanding of this, we further draw the one-parameter bifurcation diagram with $d$ as the parameter (sweeping $d$ along the broken horizontal line in Fig. 12(a)). The diagram is shown in Fig. 12(b). The birhythmic zone is presented by the shaded region in the figure.

For $d<d_{\mbox{HB}}$ the system is in stable steady states. At $d=d_{\mbox{HB}}$ a supercritical Hopf bifurcation occurs and oscillation of small amplitude emerges. For $d_{\mbox{HB}}<d<d_{L}$ the system is in a monorhythmic zone with small amplitude LC. The birhythmic zone lies between $d_{L}<d<d_{U}$. For $d>d_{U}$ the system shows monorhythmicity with large amplitude oscillation. Note that the variation of $d$ causes the system to undergo a transition, which is inverse to the previous two models discussed in earlier subsections. It is interesting to note that, there is a hysteresis governed by the SNLC bifurcations (cyan zone) with the hysteresis width $\Delta d=(d_{U}-d_{L})$. The SNLC curves meet each other at $d=d_{C}$ indicating the presence of a codimension-2 cusp type of bifurcation.

Finally, we demonstrate the time series of the system for different values of the control parameter $d$ with two different initial conditions, namely, $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ in Fig. 13. The solid line in the figure indicates the result for $\mathcal{I}_{1}$ and the line with hollow circles that for $\mathcal{I}_{2}$. The two initial conditions are chosen in such a way that the system shows small amplitude oscillation for $\mathcal{I}_{1}$ and large amplitude oscillation for $\mathcal{I}_{2}$ when $d=0$ (cf. Fig. 11).

Fig. 13(a,b) shows time series and phase plane plots, respectively, for $d=-0.012$ $(d<d_{HB})$. The system rests in a stable steady point. A supercritical Hopf bifurcation occurs at $d=d_{HB}$ and the system revives from its dormant state and LC of small amplitude emerges. The situation for $d=-0.001$ is shown in Fig. 13(c-d). In the range $d_{L}<d<d_{U}$ the system exhibits a birhythmic nature. Fig. 13(e-f) shows the case for $d=0.001$. Finally, the system creates a large amplitude LC for $d>d_{U}$. The scenario for $d=0.012$ is shown in Fig. 13(g-h). Thus we can conclude that, like the two previous models, the control scheme can effectively eliminate the birhythmic behavior and may lead to a preferred monorhythmic state.

V Conclusions

In conclusion, we have studied the effect of self-feedback mechanism on birhythmic oscillations. We have shown that the self-feedback through a proper variable is able to eliminate birhythmic oscillations and select monorhythmic oscillation of either large amplitude or small amplitude, as desired. We have successfully applied the self-feedback mechanism to three realistic models from different branches of natural science, e.g., an energy harvesting system and biochemical models such as the OAK model and glycolysis model. A rigorous analysis using harmonic decomposition and energy balance method have established the efficacy of the self-feedback mechanism. Further, we have explored the possible bifurcation scenarios to get a deep understanding of the genesis of monorhythmic oscillation that comes out of birhythmic oscillation. In comparison with our previously proposed control scheme in Ref. Biswas, Banerjee, and Kurths, 2016 the proposed self-feedback control technique uses self-feedback of only a single variable, hence the physical implementation of this control scheme is comparatively easy. We believe that our proposed coupling scheme is general enough to be applied effectively to control birhythmicity in several physical and biochemical processes.

Acknowledgements.

Authors thankfully acknowledge the insightful suggestions by the anonymous referees. DB acknowledges CSIR, New Delhi, India. TB acknowledges Science and Engineering Research Board (Department of Science and Technology, India) [grant No. SB/FTP/PS-05/2013]. DB acknowledges Haradhan Kundu, Department of Mathematics, University of Burdwan, for his useful suggestions regarding computations.

Appendix A Stability analysis using Poincaré-Lindstedt method

As we discussed in the main text, the harmonic decomposition technique is not an asymptotic method and in this context the Poincaré-Lindstedt technique is a much more suitable technique for weakly nonlinear systems. Since we are interested to find out the amplitude equation of Eq. (6), therefore in the following we show that both the techniques give similar amplitude equations.

Let us consider the system equation Eq. (6). With $\Omega_{0}=1$ it becomes

[TABLE]

With $F(\dot{y})$ given by Eq. (2). We consider the solution to be of the form

[TABLE]

where $\tau\equiv\omega t$ and $\omega$ is a known frequency and $y_{i}(\tau)$, ($i=1,2,\dots$) is the periodic function of periodic $2\pi$. We also write

[TABLE]

where $\omega_{i}$ are unknown ad needed to determine. Further, we decompose the control parameter $d$ as

[TABLE]

Substituting Eqs. (24), (25) and (26) into Eq. (23) and equating the coefficients of different powers of $\mu$, we get

Coeff. of $\mu^{0}$:

[TABLE]

Coeff. of $\mu^{1}$:

[TABLE]

Coeff. of $\mu^{2}$:

[TABLE]

The solution of Eq. (27) is given by

[TABLE]

where $A$ is the complex amplitude to be determined. Substitution of Eq. (30) in Eq. (29) invokes the solvability condition as

[TABLE]

Equating the real and imaginary parts of the above equation

[TABLE]

and

[TABLE]

Therefore, with $A$ replaced by $\frac{A}{2}$, the amplitude equation of Eq. (33) is similar to that obtained in Eq. (12) (using the harmonic decomposition technique). Apart from amplitude equation one can show $\omega_{2}\neq 0$ (note that $\omega_{1}=0$) by applying the similar steps carried above but for the coefficient of $\mu^{2}$: However, this analysis does not affect our amplitude equations.