Two Nested Limit Cycles in Two-Species Reactions

Abstract

We search for limit cycles in the dynamical model of two-species chemical reactions that contain seven reaction rate coefficients as parameters and at least one third-order reaction step, that is, the induced kinetic differential equation of the reaction is a planar cubic differential system. Symbolic calculations were carried out using the Mathematica computer algebra system, and it was also used for the numerical verifications to show the following facts: the kinetic differential equations of these reactions each have two limit cycles surrounding the stationary point of focus type in the positive quadrant. In the case of Model 1, the outer limit cycle is stable and the inner one is unstable, which appears in a supercritical Hopf bifurcation. Moreover, the oscillations in a neighborhood of the outer limit cycle are slow-fast oscillations. In the case of Model 2, the outer limit cycle is unstable and the inner one is stable. With another set of parameters, the outer limit cycle can be made stable and the inner one unstable.

1. Introduction

This paper is a part of a series of works [1,2,3,4,5] on the existence or absence of limit cycles in two- and three-species chemical reactions endowed with mass action kinetics. A very recent paper [6] by Valenzuela et al. is similar to ours, but they use the normal form theory, which is computationally not very efficient. (A few of our papers on models with non mass action—rational—kinetics are [7,8].)

Let us review the history of the topics shortly. Although similar models were known even in the XIXth century [9], the first model with some chemical relevance and showing oscillatory behavior was the Lotka–Volterra reaction [10,11,12].

$$\mathrm{X}\to 2\mathrm{X},\mathrm{X}+\mathrm{Y}\to 2\mathrm{Y},\mathrm{Y}\to 0.$$

Frank-Kamenetsky [13] used it to describe the oxidation of hydrocarbons, or—in [14]—as a model of cold flames. The induced kinetic differential equation of this reaction shows conservative oscillations: its stationary point is a center around which a family of periodic trajectories—parametrized by the initial conditions—appears. People interested in applications were looking for a natural model with (stable) limit cycles, because this corresponds to the experimental observation that the trajectories tend to a periodic one on the long run. According to the general view, the first reaction (called Brusselator) with an induced kinetic differential equation having a limit cycle was constructed by Prigogine and Lefever in [15]. However, Frank-Kamenetsky and Salnikov in an extremely well-written early paper [16] constructed the reaction.

$$\mathrm{X}\underset{}{\overset{k_1}{\to }}2\mathrm{X},\mathrm{X}+\mathrm{Y}\underset{}{\overset{k_2}{\to }}2\mathrm{Y},2\mathrm{X}\underset{}{\overset{k_4}{\to }}3\mathrm{X},\mathrm{Y}\underset{k_5}{\overset{k_3}{\leftrightarrow }}0$$

with the induced kinetic differential equation

$$\dot{x}=k_{1}x-k_{2}xy+k_4{x}^2\dot{y}=k_{2}xy-k_{3}y+k_5$$

having a limit cycle and only containing second-order terms (second-order reaction steps).

In the sixties and seventies of the twentieth century, the oscillatory Belousov–Zhabotinsky reaction was the topic of active (mainly experimental) research [17,18]. Ding-Hsü [19] was the first to use explicitly Hopf’s theorem to prove the existence of a limit cycle in the induced kinetic differential equation of a reaction, namely that of the Oregonator, leading to a three-variable differential equation with a second degree right hand side. It is worth citing him to show how happy he was to find this tool for this purpose: “thereby to publicize Hopf’s theorem. This theorem is not as well known and available as it should be”.

As to the more theoretical investigations, one should mention the Póta–Hanusse–Tyson–Light theorem, actually proven by Póta [20] stating that two-species bimolecular systems cannot have limit cycles (see an alternative proof in [21]). Schnakenberg [22] and following him Császár et al. [23] constructed classes of reactions which showed limit cycles—numerically.

A possible classification of problems in reaction kinetics (but also in other fields of applied sciences) might be this [24] [Ch. 11]. Models are formulated and their properties are investigated by the methods of the qualitative theory, or by numerical and other methods: this approach is the direct approach. A more interesting direction for the chemist (biologist, etc.) is the inverse approach: given some measurements of qualitative properties of a phenomenon, we look for models with this behavior. A special case of this is parameter estimation, when the model structure is known and it is only the numerical values of the model parameters we are looking for. All the research mentioned up to this point solved direct problems: given a reaction or its kinetic differential equation, some properties of it are found. (Note that whereas the route from reactions to differential equations is straightforward [24,25], it is far from being true in the opposite direction, see, e.g., [26].) Our present paper also belongs to this category. However, Escher [27,28] formulated and solved a series of inverse problems: given a dynamic behavior (e.g., existence of a limit cycle) he constructed reactions to show the given phenomenon.

In this paper, we present two reactions among two species and not higher than third-order reaction steps, then show that both models can have two limit cycles with appropriate values of the parameters (reaction rate coefficients by meaning). Although we present illustrating figures coming from numerical calculations, we emphasize that our proofs are rigorous and use no approximations. The novelty of our treatment is that we use focus quantities to find limit cycles symbolically in a kinetic differential equation. Although the framework of our work is simple, the computations are cumbersome both for the human and the computer, and also needs ingenuity.

The structure of our paper is as follows. Section 2 describes a method to find limit cycles. Next, in Section 3 we discuss our Model 1 in detail. For Model 2, the same investigations will be carried out in less detail in Section 4. Figures illustrate the rigorous results throughout. Finally, we mention that our Mathematica notebooks are available online [29] and also as Supplementary Materials.

2. A Method to Find Limit Cycles

Here—following [30]—we summarize briefly how we can find two limit cycles in the planar differential system $\dot{u}=P\circ (u,v),\dot{v}=Q\circ (u,v)$, where P and Q are functions (here: polynomials) dependent on some parameters, by a local investigation of a singular point.

• As a first step, the singular point $(u_0,v_0)$ of the system is shifted into the origin using the substitution $x=u-u_0,y=v-v_0,$ so in the new coordinates the system is written as

  $$\dot{x}=\tilde{P}\circ (x,y),\dot{y}=\tilde{Q}\circ (x,y).$$

  (1)
  Let us denote the Jacobian matrix of this system at the origin by J and let ${\lambda }_1$ and ${\lambda }_2$ be the eigenvalues. Next, the parameters are chosen in such a way that $trace\left(J\right)=0.$
• Then we look for a polynomial ${\displaystyle \Phi (x,y)=\sum _{k+s=2}^6{\varphi }_{ks}{x}^k{y}^s}$ such that

  $$\frac{\partial \Phi }{\partial x}\dot{x}+\frac{\partial \Phi }{\partial y}\dot{y}=g_1{(x^2+y^2)}^2+g_2{(x^2+y^2)}^3+h.o.t$$

  (2)

  and the quadratic part of $\Phi (x,y)$,

  $${\Phi }_0(x,y)={\varphi }_{20}{x}^2+{\varphi }_{11}xy+{\varphi }_{02}{y}^2$$

  is a positive definite quadratic form. The coefficients $g_1$ and $g_2$ in (2) depend on parameters of system (1) and are called focus (or Lyapunov) quantities.
• Keeping $trace\left(J\right)=0$ and ${\Phi }_0$ positive definite we look for values of parameters of system (1) to set the values of $g_1$ and $g_2$ in the following way.

  (a)

  First, if  $g_1=0$ and $g_2<0$ then, since $\Phi$ is a positive definite Lyapunov function, the origin is a stable focus.

  (b)

  If we now take a small perturbation, so that $g_1$ becomes positive (while $g_2$ remains negative), then an unstable focus arises at the origin and a stable limit cycle appears around the singular point.

• Finally, the parameters are perturbed in such a way that $trace\left(J\right)\ne 0$ and $\mathrm{Re}\left({\lambda }_{1,2}\right)<0$. In this case, the origin becomes stable, and if the perturbation is sufficiently small, then the outer stable limit cycle is preserved (but can be shifted) and an unstable limit cycle appears between the origin and the outer stable limit cycle as a result of a supercritical Hopf bifurcation. Since we cannot say in advance what perturbations are “sufficiently small”, the existence of two limit cycles in a specific perturbed system should be also verified numerically.
  Similarly, if in step 3a we look for parameters such that $g_1=0$ and $g_2>0$ at the beginning and achieve that $\mathrm{Re}\left({\lambda }_{1,2}\right)>0$, $g_1<0$, $g_2>0$ in the end, then the outer limit cycle will be unstable and the inner one will be stable.

3. Model 1

We investigate the induced kinetic differential equation of the reaction in Figure 1.

Assuming mass action type kinetics, i.e., the dynamical system:

$$\begin{array}{cc}\hfill \dot{x}=& x^{2}y+xy-c_1{x}^2-d_{1}x+e_{1}y+f_1,\hfill \\ \hfill \dot{y}=& -x^{2}y-xy+c_1{x}^2+d_{2}x-e_{2}y+f_2.\hfill \end{array}$$

(3)

where $x\left(t\right)\ge 0$ and $y\left(t\right)\ge 0$ denote the concentrations of two chemical species and $c_1,d_1,d_2,e_1,e_2,f_1,f_2$ are the reaction rate coefficients, all supposed to be positive. The meaning of reaction rate coefficients can be found in standard textbooks on chemical kinetics or in [24]. The reader who is not an insider will understand the signs if (s)he glances a look at Figure 1.

3.1. Symbolic Preparations

Let us denote the singular point of system (3) in the first quadrant by $A(x_0,y_0).$ To simplify the calculations, we consider the case when $x_0=1.$ Solving the system $\dot{x}=0,\dot{y}=0$ for $d_1$ and $y_0,$ we get that

$$\begin{array}{cc}\hfill d_1=& \frac{d_2(2+e_1)+c_1(e_1-e_2)+(2+e_2)f_1+(2+e_1)f_2}{2+e_2},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill y_0=& \frac{c_1+d_2+f_2}{2+e_2}.\hfill \end{array}$$

(5)

Let us emphasize again that the assumption $x_0=1$ implies that the reaction rate coefficients are not independent, (4) should hold among them.

Now, shifting the singular point $A(x_0,y_0)$ into the origin with the transformation $x_1=x-x_0,y_1=y-y_0,$ we get that the transformed system is

$$\begin{array}{cc}\hfill {\dot{x}}_1=& -\frac{1}{2+e_2}(c_1{x}_1-d_2{x}_1+c_1{e}_1{x}_1+d_2{e}_1{x}_1+c_1{e}_2{x}_1+2f_1{x}_1\hfill \\ \hfill & +e_2{f}_1{x}_1-f_2{x}_1+e_1{f}_2{x}_1+c_1{x}_1^2-d_2{x}_1^2+c_1{e}_2{x}_1^2-f_2{x}_1^2-4y_1\hfill \\ \hfill & -2e_1{y}_1-2e_2{y}_1-e_1{e}_2{y}_1-6x_1{y}_1-3e_2{x}_1{y}_1-2x_1^2{y}_1-e_2{x}_1^2{y}_1),\hfill \\ \hfill {\dot{y}}_1=& -\frac{1}{2+e_2}(-c_1{x}_1+d_2{x}_1-2c_1{e}_2{x}_1-d_2{e}_2{x}_1+3f_2{x}_1-c_1{x}_1^2\hfill \\ \hfill & +d_2{x}_1^2-c_1{e}_2{x}_1^2+f_2{x}_1^2+4y_1+4e_2{y}_1+e_2^2{y}_1+6x_1{y}_1+3e_2{x}_1{y}_1\hfill \\ \hfill & +2x_1^2{y}_1+e_2{x}_1^2{y}_1).\hfill \end{array}$$

(6)

Let J denote the Jacobian matrix of system (6) at the origin. The necessary condition for a Hopf bifurcation at the origin (and also at the point $A(x_0,y_0)$) is that the matrix J has pure imaginary eigenvalues (and the necessary condition for this is that the trace of J is zero).

Calculating the trace of $J,$ we get that

$$trace\left(J\right)=-\frac{4+c_1-d_2+c_1{e}_1+d_2{e}_1+4e_2+c_1{e}_2+e_2^2+2f_1+e_2{f}_1-f_2+e_1{f}_2}{2+e_2}$$

(7)

Solving $\mathrm{trace}\left(J\right)=0$ for $c_1$ gives

$$c_1=\frac{d_2-d_2{e}_1-(2+e_2)(2+e_2+f_1)+f_2-e_1{f}_2}{1+e_1+e_2}.$$

(8)

In this case, the eigenvalues of J are of the form ${\lambda }_{1,2}=\pm i\sqrt{\beta }$, where

$$\begin{array}{cc}\hfill \beta =& \frac{1}{1+e_1+e_2}(-2e_1+2d_2{e}_1+d_2{e}_1^2+2e_2-2d_2{e}_2+e_1{e}_2-d_2{e}_1{e}_2-e_2^2+e_1{e}_2^2-e_2^3\hfill \\ & +2f_1+e_1{f}_1+4e_2{f}_1+2e_1{e}_2{f}_1+2f_2+5e_1{f}_2+2e_1^2{f}_2).\hfill \end{array}$$

(9)

The eigenvalues will be purely imaginary if $\beta >0$. Then system (6) with the substitution for $c_1$ in (8) (and using the variables x and y instead of $x_1$ and $y_1$) will be the following:

$$\begin{array}{cc}\hfill \dot{x}=& \frac{1}{1+e_1+e_2}(2x+2e_{1}x+3e_{2}x+e_1{e}_{2}x+e_2^{2}x+2x^2+d_2{e}_1{x}^2\hfill \\ & +3e_2{x}^2+e_2^2{x}^2+f_1{x}^2+e_2{f}_1{x}^2+e_1{f}_2{x}^2+2y+3e_{1}y+e_1^{2}y\hfill \\ & +2e_{2}y+e_1{e}_{2}y+3xy+3e_{1}xy+3e_{2}xy+x^{2}y+e_1{x}^{2}y+e_2{x}^{2}y),\hfill \\ \hfill \dot{y}=& -\frac{1}{1+e_1+e_2}(2x+d_2{e}_{1}x+5e_{2}x-d_2{e}_{2}x+2e_2^{2}x+f_{1}x+2e_2{f}_{1}x\hfill \\ & +f_{2}x+2e_1{f}_{2}x+2x^2+d_2{e}_1{x}^2+3e_2{x}^2+e_2^2{x}^2+f_1{x}^2+e_2{f}_1{x}^2\hfill \\ & +e_1{f}_2{x}^2+2y+2e_{1}y+3e_{2}y+e_1{e}_{2}y+e_2^{2}y+3xy+3e_{1}xy+3e_{2}xy+x^{2}y+e_1{x}^{2}y+e_2{x}^{2}y).\hfill \end{array}$$

(10)

Now, in order to apply the procedure described in the previous section, we look for a polynomial

$$\Phi (x,y)=\sum _{k+s=2}^6{\varphi }_{ks}{x}^k{y}^s$$

(11)

such that

$$\frac{\partial \Phi }{\partial x}\dot{x}+\frac{\partial \Phi }{\partial y}\dot{y}=g_1{(x^2+y^2)}^2+g_2{(x^2+y^2)}^3+h.o.t.$$

(12)

Calculating the quadratic part of (11) we get that

$${\Phi }_2=\frac{1}{2}{\varphi }_{11}\left(\frac{(d_2(e_1-e_2)+(1+2e_2)(2+e_2+f_1)+f_2+2e_1{f}_2)}{(2+e_2)(1+e_1+e_2)}{x}^2+2xy+\frac{2+e_1}{2+e_2}{y}^2\right),$$

where ${\varphi }_{11}$ can be chosen arbitrarily, so we set ${\varphi }_{11}=1$. Now, let

$$A:=\left(\begin{array}{cc}\frac{d_2(e_1-e_2)+(1+2e_2)(2+e_2+f_1)+f_2+2e_1{f}_2}{2(2+e_2)(1+e_1+e_2)}& \frac{1}{2}\\ \frac{1}{2}& \frac{2+e_1}{2(2+e_2)}\end{array}\right)$$

The determinants of the upper-left main minors of A are

$$\begin{array}{cc}\hfill A_{11}=& \frac{d_2(e_1-e_2)+(1+2e_2)(2+e_2+f_1)+f_2+2e_1{f}_2}{(1+e_1-e_2)(2(2+e_2)(1+e_1+e_2)},\hfill \\ \hfill det\left(A\right)=& -\left(\right(2e_1-2d_2{e}_1-d_2{e}_1^2-2e_2+2d_2{e}_2-e_1{e}_2+d_2{e}_1{e}_2\hfill \\ & +e_2^2-e_1{e}_2^2+e_2^3-2f_1-e_1{f}_1-4e_2{f}_1-2e_1{e}_2{f}_1\hfill \\ & -2f_2-5e_1{f}_2-2e_1^2{f}_2)/\left(4{(2+e_2)}^2(1+e_1+e_2)\right)).\hfill \end{array}$$

(13)

At this point, the positive definiteness of A together with the previous conditions cannot be decided, so in order to simplify the further calculations, we fix the values of two more reaction rate coefficients $f_1$ and $f_2$ as

$$f_1=1,f_2=2$$

(14)

and proceed with calculating $g_1$. This gives

$$\begin{array}{cc}\hfill g_1=& (-36-105e_1-30d_2{e}_1-120e_1^2-50d_2{e}_1^2-6d_2^2{e}_1^2-56e_1^3-28d_2{e}_1^3\hfill \\ \hfill & -7d_2^2{e}_1^3-6e_1^4-6d_2{e}_1^4-2d_2^2{e}_1^4-75e_2+12d_2{e}_2-159e_1{e}_2-8d_2{e}_1{e}_2\hfill \\ \hfill & +6d_2^2{e}_1{e}_2-113e_1^2{e}_2-9d_2{e}_1^2{e}_2+7d_2^2{e}_1^2{e}_2-18e_1^3{e}_2-d_2{e}_1^3{e}_2+2d_2^2{e}_1^3{e}_2\hfill \\ \hfill & -36e_2^2+13d_2{e}_2^2-65e_1{e}_2^2+18d_2{e}_1{e}_2^2-30e_1^2{e}_2^2-6e_1^3{e}_2^2-4d_2{e}_1^3{e}_2^2\hfill \\ \hfill & +9e_2^3+d_2{e}_2^3-6e_1{e}_2^3+7d_2{e}_1{e}_2^3-e_1^2{e}_2^3+4d_2{e}_1^2{e}_2^3+6e_2^4+7e_1{e}_2^4\hfill \\ \hfill & -2e_1^2{e}_2^4+2e_1{e}_2^5)/((2+e_2)(123+234e_1+34d_2{e}_1+137e_1^2+30d_2{e}_1^2\hfill \\ \hfill & +3d_2^2{e}_1^2+26e_1^3+2d_2{e}_1^3+3e_1^4+330e_2-34d_2{e}_2+348e_1{e}_2+16d_2{e}_1{e}_2\hfill \\ \hfill & -6d{2}^2{e}_1{e}_2+68e_1^2{e}_2+6e_1^3{e}_2+307e_2^2-46d_2{e}_2^2+3d_2^2{e}_2^2+126e_1{e}_2^2\hfill \\ \hfill & +10d_2{e}_1{e}_2^2+11e_1^2{e}_2^2+116e_2^3-12d_2{e}_2^3+12e_1{e}_2^3+16e_2^4\left)\right).\hfill \end{array}$$

(15)

We calculated the value of $g_2$ in (12) as a function of $d_2,e_1,e_2$, however, the expression is very complicated, it is contained in [29]. To simplify it, again we fix the values of two further reaction rate coefficients, $e_1$ and $d_2$ as

$$e_1=1/2,d_2=20,$$

(16)

then the formula for $g_2$ is

$$\begin{array}{cc}\hfill g_2=& -\left(\right(2(-126512084933352295+478399692835658985e_2\hfill \\ \hfill & -762458375238510816e_2^2+730329622432012315e_2^3\hfill \\ \hfill & -466875108129002500e_2^4+180057388265535584e_2^5\hfill \\ \hfill & -42412278548678749e_2^6+7565318705371668e_2^7-1899913140260206e_2^8\hfill \\ \hfill & +200501244715092e_2^9+44607817747872e_2^{10}+18408318367872e_2^{11}\hfill \\ \hfill & +8352995841088e_2^{12}+1190085589888e_2^{13}+348608893184e_2^{14}\hfill \\ \hfill & +33760183296e_2^{15}+6787676672e_2^{16}+473953280e_2^{17}+64299008e_2^{18}\hfill \\ \hfill & +2555904e_2^{19}+196608e_2^{20}\left)\right)/(3(2+e_2)(83-42e_2+8e_2^2)(-73+85e_2\hfill \\ \hfill & (+e_2^2+2e_2^3){)}^2(109-25e_2+32e_2^2+4e_2^3)(17163-19172e_2+12044e_2^2\hfill \\ \hfill & -1888e_2^3+256e_2^4)(23933-29628e_2+23764e_2^2-2976e_2^3+512e_2^4)\left)\right).\hfill \end{array}$$

(17)

Using the requirement that negative cross-effect cannot be present in a kinetic differential equation [24] [Theorem 6.27] and conditions (4), (5), (8), (9), (13)–(17) and the Reduce command of Mathematica [31], the numerical solution of the semi-algebraic system

$$g_1=0\wedge g_2<0\wedge A_{11}>0\wedge \det \left(A\right)>0\wedge y_0>0\wedge \beta >0\wedge c_1>0\wedge d_1>0\wedge e_2>0$$

for $e_2$ is

$$e_2\approx -0.771291.$$

(18)

Similarly, we found that $g_2>0$ is not possible for $g_2$ defined by (17). For the value of $e_2$ in (18), the numerical value of $g_2$ is

$$g_2\approx -0.0278896$$

and also the condition $\mathrm{trace}\left(J\right)=0$ holds. From $A_{11}>0$ and $\det \left(A\right)>0$ it follows that ${\Phi }_2$ is a positive definite quadratic form, so the function $\Phi (x,y)$ defined in (11) is a positive definite Lyapunov function in a sufficiently small neighbourhood of the origin. This together with conditions $g_1=0$ and $g_2<0$ means that its Lie derivative given in (12) is negative definite, so the origin is a stable focus.

From (15), we can see that $g_1$ is a polynomial in $d_2,e_1,e_2$ and $d_2=20,e_1=1/2$ and $e_2$ defined by (18) is a simple root of the equation $g_1=0$. Thus, we can choose a small perturbation of any of the parameters $d_2,e_1,e_2$ such that $g_1$ becomes negative. Hence, after such a perturbation a stable limit cycle bifurcates from the origin.

Now, from (7) and (8) it is obvious that we can take an arbitrarily small perturbation of any of the parameters $d_2$, $e_i$, $f_i$ ($i=1,2$) such that an unstable limit cycle appears from the origin in a supercritical Hopf bifurcation. Thus, the resulting perturbed system has two limit cycles.

Remark 1.

We mention that we have performed some computational experiments changing the parameters $f_1$, $f_2$, $d_2,e_1,e_2$ trying to find values of the parameters for which $\mathrm{trace}\left(J\right)=0,$$g_1=0$, $g_2>0$, but we failed to find such values.

3.2. Numerical Results for Model 1

In this section we present the numerical study confirming the existence of two limit cycles in system (3) and illustrate the results with figures created with the Wolfram Language.

3.2.1. The Appearance of the First Limit Cycle

For the values of the parameters defined by (4), (8), (14), (16) and (18) we first perturb the parameter $e_2$ such that $g_1$ becomes positive. If

$$e_2=\frac{78}{100}$$

then we get that

$$g_1\approx 0.015511>0,g_2\approx -0.666999<0,\mathrm{trace}\left(J\right)=0.$$

It means that the origin becomes unstable and a stable limit cycle appears around the singular point.

In this case, the parameter settings for system (3) are the following: $c_1={\displaystyle \frac{1229}{5700}}\approx 0.215614$, $d_1={\displaystyle \frac{59173}{2850}}\approx 20.7625$, $d_2=20$, $e_1={\displaystyle \frac{1}{2}}$, $e_2={\displaystyle \frac{78}{100}}$, $f_1=1$, $f_2=2$. The singular point of the system in the first quadrant is $A(x_0,y_0)=(1,7.99123)$ and the eigenvalues of the Jacobian matrix are ${\lambda }_{1,2}\approx \pm 1.06195i$.

The trajectories can be seen in Figure 2 and Figure 3, where the initial point $(x_0,y_0+d)$ is denoted in red and $A(x_0,y_0)$ is denoted in orange. Picture (a) shows the behavior of the trajectories in the phase plain, and pictures (b) and (c) show the solutions as a function of time. Where the shape of the trajectories makes it possible, we also denoted the direction of the trajectories with an arrow.

3.2.2. The Appearance of the Second Limit Cycle

Next, we perturb the parameter $c_1$ such that $\mathrm{trace}\left(J\right)$ becomes negative. If $c_1={\displaystyle \frac{22}{100}}$, then we get that

$$g_1\approx 0.015511>0,g_2\approx -0.666999<0,\mathrm{trace}\left(J\right)\approx -0.00359712<0.$$

Since the trace of J is now negative, it means that the eigenvalues of J have negative real parts. It means that the origin becomes stable and an unstable limit cycle appears between the origin and the outer stable limit cycle as a result of a supercritical Hopf bifurcation. The perturbation for $c_1$ must be sufficiently small to ensure that when the inner limit cycle appears, the outer one is still preserved.

In this case, the parameter settings for system (3) are the following: $c_1={\displaystyle \frac{22}{100}}$, $d_1={\displaystyle \frac{72148}{3475}}\approx 20.762$, $d_2=20$, $e_1={\displaystyle \frac{1}{2}}$, $e_2={\displaystyle \frac{78}{100}}$, $f_1=1$, $f_2=2$. The singular point of the system in the first quadrant is $A(x_0,y_0)=(1,7.99281)$ and the eigenvalues of the Jacobian matrix are ${\lambda }_{12}\approx -0.00179856\pm 1.0619i$.

The trajectories can be seen in Figure 4, Figure 5, Figure 6 and Figure 7, where the initial point $(x_0,y_0+d)$ is denoted in red and $A(x_0,y_0)$ is denoted in orange. Picture (a) shows the behavior of the trajectories in the phase plain, and pictures (b) and (c) show the solutions as a function of time.

4. Model 2

We investigate the following dynamical system:

$$\begin{array}{cc}\hfill \dot{x}=& \phantom{-}{x}^{2}y-xy-c_1{x}^2-d_{1}x+e_{1}y+f_1,\hfill \\ \hfill \dot{y}=& -x^{2}y+xy+c_1{x}^2+d_{2}x-e_{2}y+f_2.\hfill \end{array}$$

(19)

where $x\left(t\right)\ge 0$ and $y\left(t\right)\ge 0$ denote the concentrations of two chemical species and

$$c_1,d_1,d_2,e_1,e_2,f_1,f_2\ge 0.$$

(20)

are the reaction rate coefficients. Actually, the system is the induced kinetic differential equation of the reaction in Figure 8 assuming mass action type kinetics.

Compared to Model 1, here the signs of the terms $xy$ are changed. After repeating the steps (4)–(13) with this new system, we get that

$$\begin{array}{ccc}\hfill d_1& =& {\displaystyle \frac{c_1(e_1-e_2)+e_2{f}_1+e_1(d_2+f_2)}{e_2}}\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill c_1& =& {\displaystyle \frac{d_2-d_2{e}_1-e_2(e_2+f_1)+f_2-e_1{f}_2}{-1+e_1+e_2}}.\hfill \end{array}$$

(22)

Setting the values $f_1$ and $f_2$ as

$$f_1={\displaystyle \frac{1}{2}},f_2={\displaystyle \frac{1}{10}}$$

(23)

and proceeding with the calculations, we obtain that both $g_1$ and $g_2$ are functions of $d_2,e_1,e_2$. Investigating the possible values of $e_1$ and $d_2$ numerically, we find that when for example $d_2=1/100$ and $0.128<e_1<0.213$ then $g_1=0,g_2>0$ is possible; and when $d_2=1/100$ and $0.264<e_1<0.5$ or $0.213<e_1<0.228$ then $g_1=0,g_2<0$ is possible.

Here, we want to achieve that the outer limit cycle is unstable (that is, $g_2>0$), so we choose these parameters as

$$d_2={\displaystyle \frac{1}{100}},e_1={\displaystyle \frac{18}{100}}.$$

(24)

With this, we find that $g_1=0,g_2<0$ and $g_1=g_2=0$ cannot be the case and $g_1=0,g_2>0$ is only possible when

$$e_2\approx 0.29582$$

(25)

where $e_2$ is a root of a fifth-degree polynomial. With the parameter settings (21)–(25) we get that the system has three singular points in the first quadrant: $A(1,1.30837)$, $B(0.809127,2.04744)$ and $C(0.457223,3.41004)$ and

$$g_1\approx 0g_2\approx 0.0276312>0\mathrm{trace}\left(J\right)\approx 0$$

where J denotes the Jacobian matrix at A. In this case the point A is an unstable focus.

As a next step, keeping $c_1,d_1,d_2,e_1,f_1,f_2$ the same as in (21)–(24), we perturb $e_2$ as

$$\begin{array}{c}\hfill e_2={\displaystyle \frac{278}{1000}}\end{array}$$

(26)

and obtain that the system has three singular points in the first quadrant: $A(1,1.23247)$, $B(0.7895,2.26181)$ and $C(0.414968,4.09327)$ and

$$\begin{array}{ccc}\hfill g_1& \approx & -0.0090896<0\hfill \\ \hfill g_2& \approx & 0.0400065>0\hfill \\ \hfill \mathrm{trace}\left(J\right)& \approx & 0\hfill \end{array}$$

where J denotes the Jacobian matrix at A. In this case, the point A becomes stable and an unstable limit cycle appears around A.

Finally, perturbing $c_1$ as

$$\begin{array}{c}\hfill c_1={\displaystyle \frac{233}{1000}}\end{array}$$

(27)

we get that the system has three singular points in the first quadrant: $A(1,1.23381)$, $B(0.789796,2.26142)$ and $C(0.414926,4.09403)$ and

$$g_1\approx -0.0090896<0g_2\approx 0.0400065>0\mathrm{trace}\left(J\right)\approx 0.000726619>0$$

where J denotes the Jacobian matrix at A. In this case, the point A becomes unstable and a stable limit cycle appears between A and the outer unstable limit cycle. We have to check numerically that the perturbations are small enough, that is, the outer cycle is preserved when the inner cycle appears. This is shown in Figure 9, Figure 10 and Figure 11.

To sum up the final step, the parameter settings are $c_1={\displaystyle \frac{233}{1000}}$, $d_1={\displaystyle \frac{67983}{139000}}\approx 0.489086$, $d_2={\displaystyle \frac{1}{100}}$, $e_1={\displaystyle \frac{18}{100}}$, $e_2={\displaystyle \frac{278}{1000}}$, $f_1={\displaystyle \frac{1}{2}}$, $f_2={\displaystyle \frac{1}{10}}$. The system has three singular points in the first quadrant: $A(1,1.23381)$ (unstable focus, orange), $B(0.789796,2.26142)$ (saddle, green), $C(0.414926,4.09403)$ (real sink, blue). The initial point is denoted by red and the distance of the initial point from the singular point by d.

Finally, we would like to remark that in Model 2 it is also possible that the outer limit cycle is stable and inner one is unstable. To achieve this, we repeated the procedure described above with the following parameter settings: $c_1={\displaystyle \frac{79}{100}}$, $d_1={\displaystyle \frac{26331}{21100}}\approx 1.24791$, $d_2={\displaystyle \frac{4}{10}}$, $e_1={\displaystyle \frac{3}{10}}$, $e_2={\displaystyle \frac{633}{1000}}$, $f_1=1$, $f_2=1$. We obtained that the system has one singular point in the first quadrant: $A(x_0,y_0)=(1,3.45972)$ and

$$g_1\approx 0.00642203>0g_2\approx -0.301207<0\mathrm{trace}\left(J\right)\approx -0.00119905<0$$

where J denotes the Jacobian matrix at A which is a stable focus. The trajectories can be seen in Figure 12, Figure 13 and Figure 14, where the initial point $(x_0,y_0+d)$ is denoted in red and $A(x_0,y_0)$ is denoted in orange.

5. Discussion

We investigated two reaction networks, and have found two nested limit cycles in both of them. We did the calculations symbolically, as far as it did not become impossible. Recent methods in the qualitative theory of differential equations helped to do this.

Either with the development of computers or theory, we hope that larger parts of such investigations will be possible to be carried out symbolically, thus making possible the mapping of the whole parameter space. Furthermore, we hope to understand better the role of the common $x^{2}y$ term.