Two Nested Limit Cycles in Two-Species Reactions

: We search for limit cycles in the dynamical model of two-species chemical reactions that contain seven reaction rate coefﬁcients 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 veriﬁcations 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.


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 twoand 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].
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.
with the induced kinetic differential equatioṅ x = k 1 x − k 2 xy + k 4 x 2ẏ = 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 a supplementary material.

A Method to Find Limit Cycles
Here-following [30]-we summarize briefly how we can find two limit cycles in the planar differential systemu = P • (u, v),v = Q • (u, v), where P and Q are functions (here: polynomials) dependent on some parameters, by a local investigation of a singular point.
1. 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 aṡ Let us denote the Jacobian matrix of this system at the origin by J and let λ 1 and λ 2 be the eigenvalues. Next, the parameters are chosen in such a way that trace(J) = 0.
2. Then we look for a polynomial Φ(x, y) = 6 ∑ k+s=2 φ ks x k y s such that ∂Φ ∂xẋ and the quadratic part of Φ(x, y), 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. 3. Keeping trace(J) = 0 and Φ 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 Φ 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.
4. Finally, the parameters are perturbed in such a way that trace(J) = 0 and Re(λ 1,2 ) < 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 Re(λ 1,2 ) > 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.

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: where x(t) ≥ 0 and y(t) ≥ 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.

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ẋ = 0,ẏ = 0 for d 1 and y 0 , we get that 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 iṡ 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).
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 then the formula for g 2 is ).
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 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 ≈ −0.0278896 and also the condition trace(J) = 0 holds. From A 11 > 0 and det(A) > 0 it follows that Φ 2 is a positive definite quadratic form, so the function Φ(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 trace(J) = 0, g 1 = 0, g 2 > 0, but we failed to find such values.

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.
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 = 1229 5700 ≈ 0.215614,

The Appearance of the Second Limit Cycle
Next, we perturb the parameter c 1 such that trace(J) becomes negative. If c 1 = 22 100 , then we get that 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)     y (c) Figure 5. Model 1 with a stable outer and an unstable inner limit cycle. The trajectory is going outward, approaching the large limit cycle as t tends to +∞ and the small limit cycle as t tends to −∞. The distance of the initial point from the singular point is d = 0.042.  6. Model 1 with a stable outer and an unstable inner limit cycle. The trajectory is going outward, approaching the large limit cycle as t tends to +∞ and the small limit cycle as t tends to −∞. The distance of the initial point from the singular point is d = 0.038.  Figure 7. Model 1 with a stable outer and an unstable inner limit cycle. The trajectory is going inward, approaching the singular point as t tends to +∞ and the small limit cycle as t tends to −∞. The distance of the initial point from the singular point is d = 0.01.

Model 2
We investigate the following dynamical system: where x(t) ≥ 0 and y(t) ≥ 0 denote the concentrations of two chemical species and 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  Figure 8. The reaction inducing the system (19).
Setting the values f 1 and f 2 as 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 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 ≈ 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 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 e 2 = 278 1000 (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 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 c 1 = 233 1000 (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 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 Figures 9-11.  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 = 79 100 , d 1 = 26331 21100 ≈ 1.24791, d 2 = 4 10 , e 1 = 3 10 , e 2 = 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 where J denotes the Jacobian matrix at A which is a stable focus. The trajectories can be seen in Figures 12-14, where the initial point (x 0 , y 0 + d) is denoted in red and A(x 0 , y 0 ) is denoted in orange.  Figure 13. Model 2 with a stable outer and an unstable inner limit cycle. The trajectory is going outward, approaching the outer limit cycle as t → ∞ and the inner limit cycle as t → −∞ (d = 0.05).  Figure 14. Model 2 with a stable outer and an unstable inner limit cycle. The trajectory is going inward, approaching the point A as t → ∞ and the small limit cycle as t → −∞ (d = 0.01).

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.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.