Chemical Systems with Limit Cycles

The dynamics of a chemical reaction network (CRN) is often modeled under the assumption of mass action kinetics by a system of ordinary differential equations (ODEs) with polynomial right-hand sides that describe the time evolution of concentrations of chemical species involved. Given an arbitrarily large integer \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K \in {\mathbb N}$$\end{document}K∈N, we show that there exists a CRN such that its ODE model has at least K stable limit cycles. Such a CRN can be constructed with reactions of at most second-order provided that the number of chemical species grows linearly with K. Bounds on the minimal number of chemical species and the minimal number of chemical reactions are presented for CRNs with K stable limit cycles and at most second order or seventh-order kinetics. We also show that CRNs with only two chemical species can have K stable limit cycles, when the order of chemical reactions grows linearly with K.


Introduction
Chemical reaction networks (CRNs) are often modeled by reaction rate equations, which are systems of first-order, autonomous, ordinary differential equations (ODEs) describing the time evolution of the concentrations of chemical species involved. Considering CRNs which are subject to the law of mass action, their reaction rate equations have polynomials on their right-hand sides (Yu and Craciun 2018;Craciun et al. 2020). The mathematical investigation of ODEs with polynomial right-hand sides has a long history and includes a number of challenging open mathematical problems, for example, Hilbert's 16 th Problem (Ilyashenko 2002), which asks questions about the number and position of limit cycles of the planar ODE system of the form where f (x, y) and g(x, y) are real polynomials of degree at most n. Denoting H (n) the maximum number of limit cycles for the system (1-2), neither the value of H (n) (for n ≥ 2) nor any upper bound on H (n) has yet been found (Ilyashenko 1991). Since a quadratic system with 4 limit cycles has been constructed (Shi 1980), we know that H (2) ≥ 4. Similarly, H (3) ≥ 13, because cubic systems with at least 13 limit cycles have been found (Li et al. 2009;Yang et al. 2010).
Considering CRNs with two chemical species undergoing chemical reactions of at most n-th order, their reaction rate equations are given in the form (1-2), where f (x, y) and g(x, y) are real polynomials of degree at most n. In particular, if we denote by C(n) the maximum number of stable limit cycles in such reaction rate equations, then we have C(n) ≤ H (n). Considering CRNs with two chemical species undergoing chemical reactions of at most second order, it has been previously shown (Póta 1983;Schuman and Tóth 2003) that their reaction rate equations cannot have any limit cycles (i.e., C(2) = 0), while general ODEs with quadratic right-hand sides can have multiple limit cycles, with H (2) ≥ 4. In particular, we observe that finding CRNs with two chemical species which have, under mass action kinetics, multiple stable limit cycles, is even more challenging than finding planar polynomial ODEs with multiple limit cycles. Considering cubic systems, we have H (3) ≥ 13, but most of the chemical systems (with at most third-order reactions) in the literature often have at most one limit cycle (Field and Noyes 1974;Schnakenberg 1979;Plesa et al. 2016). A chemical system with two stable limit cycles has been constructed (Plesa et al. 2017), giving C(3) ≥ 2, but this is still much less than 13 limit cycles which can be found in some ODE systems with cubic right-hand sides in the literature (Li et al. 2009;Yang et al. 2010). To obtain multiple stable limit cycles in chemical systems, we have to consider higher-order chemical reactions or systems with more than two chemical species Hofbauer 2021, 2022).
Considering CRNs with N chemical species undergoing chemical reactions of at most n-th order, their reaction rate equations are given as the following system of ODEs where x = (x 1 , x 2 , . . . , x N ) ∈ R N is the vector of concentrations of N chemical species and its right-hand side f : R N → R N is a vector of real polynomials of degree at most n. In this paper, we prove the following first main result: Theorem 1 Let K be an arbitrary positive integer. Then there exists a CRN with N (K ) chemical species which are subject to M(K ) chemical reactions of at most seventh order such that (i) Reaction rate equations (3) have at least K stable limit cycles, (ii) We have N (K ) ≤ K + 2 and M(K ) ≤ 29 K .
Theorem 1 provides a stronger result than finding K limit cycles in a polynomial ODE system of the form (3), because not every polynomial ODE system corresponds to a CRN and, therefore, the set of reaction rate equations is a proper subset of ODEs with polynomial right-hand sides. To make the existence of K limit cycles possible while restricting to polynomials of degree at most n = 7, we allow for more than two chemical species, replacing the ODE system (1-2) by a more general ODE system (3) with N (K ) equations. In particular, the next important question is how small the CRN can be so that it has K limit cycles. Our answer is partially given in part (ii) of Theorem 1 where we provide upper bounds on the number of chemical species involved and the number of chemical reactions (of at most seventh order). Another important parameter to consider is the maximum order of the chemical reactions involved, i.e., the degree n of the polynomials on the right-hand side of the ODE system (3). Since systems of at most second-order reactions (the case n = 2) are of special interest in the theory of CRNs and applications (Wilhelm 2000), we state our second main result as: Theorem 2 Let K be an arbitrary positive integer. Then there exists a CRN with N (K ) chemical species which are subject to M(K ) chemical reactions of at most second order such that (i) Reaction rate equations (3) have at least K stable limit cycles, (ii) We have N (K ) ≤ 7K + 14 and M(K ) ≤ 42 K + 24.
By restricting to second-order (bimolecular) reactions, we obtain CRNs with more realistic second-order kinetics, but our construction increases the number of species and chemical reactions involved, as it can be seen by comparing parts (ii) of Theorems 1 and 2. The precise definitions of CRNs, mass action kinetics, reaction rate equations and limit cycles in N -dimensional systems are given in Sect. 2. In both Theorems 1 and 2, we restrict our considerations to systems described by polynomial ODEs where the degree of polynomials is bounded by a constant independent of K , i.e., we consider polynomials of the degree at most n = 7 (in Theorem 1) or at most n = 2 (in Theorem 2), and we increase the number of chemical species, N (K ), as K increases, to get K stable limit cycles. Another approach is to restrict our considerations to chemical systems with only N = 2 chemical species. In Sect. 8, we construct two-species CRNs with K stable limit cycles which include chemical reactions of at most n(K )-th order, where n(K ) = 6K − 2. This establishes our third main result: Theorem 3 Let C(n) be the maximum number of stable limit cycles of reaction rate equations (1-2) corresponding to CRNs with two chemical species undergoing chemical reactions of at most n-th order. Then we have where the floor function . denotes the integer part of a positive real number.
To prove Theorems 1, 2 and 3, we first construct a planar system given by Eqs.
(1-2), where f and g are continuous non-polynomial functions chosen in such a way that the ODE system (1-2) has K stable limit cycles in the positive quadrant [0, ∞) × [0, ∞). Such a planar non-polynomial ODE system is constructed in Sect. 3. In Sect. 4, we then increase the number of chemical species from 2 to N (K ) to transform the nonpolynomial ODE system to a polynomial one. In Sect. 5, we modify this construction by using an x-factorable transformation to arrive at reaction rate equations corresponding to a CRN (Samardzija et al. 1989). Theorem 1 is then proven in Sect. 6 by showing that the reaction rate equations have at least K stable limit cycles. This is followed by our proofs of Theorems 2 and 3 in Sects. 7 and 8, respectively.

Notation and Mathematical Terminology
Definition 1 A chemical reaction network (CRN) is defined as a collection (S, C, R) consisting of chemical species S, reaction complexes C and chemical reactions R.
We denote by N the number of chemical species and by M the number of chemical reactions, i.e., |S| = N and |R| = M. Each chemical reaction is of the form where X i , i = 1, 2, . . . , N , are chemical species, and linear combinations N i=1 ν i, j X i and N i=1 ν i, j X i of species with non-negative integers ν i, j and ν i, j are reaction complexes.
Definition 2 Let (S, C, R) be a CRN with N chemical species and M chemical reactions. Let x i (t) be the concentration of chemical species X i ∈ S, i = 1, 2, . . . , N . The time evolution of x i (t) is, under the assumption of the mass action kinetics, described by the reaction rate equations, which are written as a system of N ODEs in the form where k j , j = 1, 2, . . . , M, is a positive constant called the reaction rate for the j-th reaction.
To illustrate Definitions 1 and 2, we present a simple example system, which is known to have a stable limit cycle for certain parameter values (Schnakenberg 1979;Erban and Chapman 2020).
Example Consider a chemical reaction network with two chemical species X 1 and X 2 which are subject to the following four chemical reactions Then N = |S| = 2, M = |R| = 4 and the set of reaction complexes is given as C = {2X 1 + X 2 , 3X 1 , X 1 , X 2 , ∅}. Denote ν and ν as matrices with elements ν i, j and ν i, j . Then, both are N × M = 2 × 4 matrices given by The reaction rate equations (6) are given as the following system of two ODEs which describes the time evolution of the concentrations x 1 (t) and x 2 (t) of chemical species X 1 and X 2 , respectively. In general, the reaction rate equations (6) in Definition 2 are ODEs of the form (3), where the right-hand side f : R N → R N is a vector of real polynomials. However, not every polynomial ODE system can be obtained as the reaction rate equations of a CRN as we formalize in Lemma 1, where we provide a necessary and sufficient condition (9) when a polynomial ODE system can be written as reaction rate equations of a CRN. The condition is that any polynomial right-hand side terms not proportional to x i in the equation for chemical species X i are non-negative. The necessity of the condition (9) is shown based on the fact that any reaction rate equations of a CRN satisfy this condition. The sufficiency of the condition is provided by showing that each term satisfying the condition in a polynomial ODE system can correspond to a chemical reaction.
Lemma 1 Consider a system of N ODEs with polynomial right-hand sides describing the time evolution of x i (t), i = 1, 2, . . . , N , in the form where α i, j are real constants and ν i, j are non-negative integers, for i = 1, 2, . . . , N and j = 1, 2, . . . , M. Then the polynomial ODE system (8) can be written as reaction rate equations (6)  Proof Reaction rate equations (6) are of the form (8). The non-negativity condition (9) follows from ν i, j = 0 and the non-negativity of both k j and ν i, j in Eq. (6). Conversely, if an ODE is of the form (8) and α i, j > 0, then we can choose ν i, j = ν i, j + 1 in Eq. (6) and put the reaction rate as k j = α i, j . On the other hand, if α i, j < 0, then the condition (9) implies that ν i, j ≥ 1, because ν i, j are non-negative integers. Therefore, we can put ν i, j = ν i, j − 1 and k j = −α i, j > 0.
In this paper, we prove the existence of limit cycles in chemical systems in Sects. 6, 7 and 8 by proving the existence of limit cycles in systems of ODEs (8) with polynomial right-hand sides satisfying the condition (9). Then the approach used in the proof of Lemma 1 can be used to construct the corresponding CRN. However, the construction of a CRN corresponding to reaction rate equations is not unique. For example, consider a term of the form −x 3 1 on the right-hand side of Eq. (8). Using the construction in the proof of Lemma 1, it would correspond to the chemical reaction 3X −→ 2X with the rate constant equal to 1, but the same term can also correspond to the chemical reaction 3X −→ X with the rate constant equal to 1/2. We conclude this section by a formal definition of a stable limit cycle.
Definition 3 Consider a system of N ODEs given by (3), where their right-hand side f : In Definition 3, the constant T is the period of the limit cycle and the property (iii) states that the limit cycle attracts all trajectories which start sufficiently close to it. The distance in the property (iii) of Definition 3 is the Euclidean distance defined by dist{z, x c } = min

Planar ODE Systems with Arbitrary Number of Limit Cycles
In this section, we construct a planar ODE system of the form (1-2) with K limit cycles in the positive quadrant. It is constructed as a function of 2K parameters denoted by a 1 , a 2 , . . . , a K and b 1 , b 2 , . . . , b K , as An illustrative dynamics of the ODE system (10-11) is shown in Fig. 1 Fig. 1 a Twenty illustrative trajectories of the ODE system (10-11) for K = 4 and the parameter choices a 1 = b 1 = a 2 = b 3 = 2 and a 3 = b 2 = a 4 = b 4 = 6. As t → ∞, all presented trajectories approach one of the four limit cycles, which are plotted as the black dashed lines. b Twenty illustrative trajectories of the ODE system (10-11) for K = 4 and the parameter choices a 1 = b 1 = a 2 = b 3 = 2 and a 3 = b 2 = a 4 = b 4 = 4. As t → ∞, some trajectories converge to the stable limit cycle denoted by the black dashed line, while some trajectories, which started inside the limit cycle, converge to the stable fixed point denoted as the red dot (Color figure online) where the ODE system has four limit cycles, which is highlighted in Fig. 1a by plotting some representative trajectories. The existence of K stable limit cycles for the ODE system (10-11) can also be proven analytically, as it is done in Lemma 2.
Remark To construct the ODE system (10-11), one can first consider a planar ODE system written as dr /dt = r (1 − r 2 ) in the polar coordinate system. It has one stable limit cycle (a unit circle) and can be rewritten in the Cartesian coordinate system as To obtain the ODE system (10-11), we translate the right-hand sides of ODEs (12-13) by (x, y) → (x − a k , y − b k ) so that the corresponding limit cycle is centered at the point (a k , b k ). Then we add the K terms (corresponding to K limit cycles) together with suitable weights. Note that the denominator of the k-th term in (10-11) is approximately one when (x, y) is close to (a k , b k ), while the k-th term in (10-11) is approximately zero when (x, y) is far away from the point (a k , b k ).
In Fig. 1a, we have presented an example with K = 4 and parameter choices In particular, the distance between points (a i , b i ), i = 1, 2, 3, 4, is at least four. If we decrease this distance, then the ODE system (10-11) can have less limit cycles. This is highlighted in Fig. 1b, where we present an example with K = 4 and parameter In Fig. 1b, we observe that there is only one limit cycle, denoted as the black dashed line. This limit cycle is stable and a number of illustrative trajectories converge to this limit cycle as t → ∞. However, there is also a stable equilibrium point at (3, 3), which attracts some of the trajectories. In particular, we can only expect that the ODE system (10-11) will have K stable limit cycles provided that points (a i , b i ) are sufficiently separated. This is formally proven in Lemma 2, by defining K disjoint open sets (15), where each is an annulus with a center (a i , b i ) for i = 1, 2, . . . , K . First, we show that each annulus does not contain any equilibrium points. Then, we show that a directional vector of the ODE system (10-11) on the boundary of each annulus points inside the domain. Applying the Poincaré-Bendixson theorem (Strogatz 2015), we conclude that each annulus contains at least one stable limit cycle; thus, the ODE system (10-11) has at least K stable limit cycles.
Proof We define the sets Then the condition (14) implies that i.e., the sets i are pairwise disjoint sets. We will show that each of them contains at least one stable limit cycle. The boundary of consists of two parts: outer and inner circles defined by and respectively, that is, Define the following functions for k = 1, 2, . . . , K : Then, the ODE system (10-11) can be rewritten as where First, we will show that i for i = 1, 2, . . . , K does not contain any equilibrium points. Let us consider any point (x * , y * ) ∈ i . Substituting in the terms for k = i in (22), we obtain g(x * , y * ) = r sin θ {1 − r 2 } + r cos θ 1 + r 6 cos 6 θ + r 6 sin 6 θ + The first terms in (24) and (25) can be rewritten as whereθ = α with tan α = r 2 − 1 and π/2 < α < 3π/2 in the case of (24) and θ = α − π/2 in the case of (25). Since we have for any θ and α, at least one of the two terms expressed in the form (26) is greater than which has a minimum √ 2/9 when 1/2 < r 2 < 2. Therefore, at least one of the absolute values of the i-th components, (24) and (25) at any point (x * , y * ) ∈ i is greater than or equal to √ 2/9. Without loss of generality, we assume . We want to show that the second term in (24) (i.e., the sum) has a smaller magnitude than the first term Using |z i | ≤ c and (28), we estimate (27) as Since (x * , y * ) ∈ i and (a k , b k ) ∈ k where k = i, our assumption (14) implies that c 2 ≥ 2. Thus, (29) becomes Therefore, the magnitude of the second term in (24) is bounded by 4(K − 1)/c 3 . Since Using the assumption (14), which implies the sufficient condition (31). Therefore, (x * , y * ) is not an equilibrium point.
Next, consider an arbitrary point (x b , y b ) ∈ ∂ i1 . Let us calculate the scalar product of vectors Using (22), we obtain this scalar product as The first two terms in (34) become which has a magnitude greater than 2/9 by using (28) we can estimate the third and fourth terms in (34), namely, we have where Then the sum of the third and fourth terms in (34) is bounded by 8 √ 2(K − 1)/d 3 . Therefore, the sufficient condition that the scalar product in (34) is negative is Using the assumption (14), which implies the sufficient condition (36). Therefore, the vector always points inside the domain i for each boundary point ( Similarly, for an arbitrary point (x b , y b ) ∈ ∂ i2 , we can show that the scalar product of vectors in (33) is always positive due to that the sum of the first two terms in (34) is equal to which is greater than 2/9 by using (28) with c 2 = 1/2, and the sum of the third and fourth terms in (34) is bounded by 8(K − 1)/(d 3 √ 2). Therefore, the sufficient condition that the scalar product in (34) is positive is which is a weaker condition than the condition (36), i.e., it is again satisfied because of our assumption (14). This implies that the scalar product in (34) is positive. Thus, the directional vector always points inside the domain i on all parts of the boundary ∂ i . Therefore, applying Poincaré-Bendixson theorem (Strogatz 2015), we conclude that each i contains at least one stable limit cycle. Since i , i = 1, 2, . . . , K , are pairwise disjoint, this implies that the ODE system (10-11) has at least K stable limit cycles.

ODE Systems with Polynomial Right-Hand Sides and Arbitrary Number of Limit Cycles
Considering an auxiliary variable we can formally convert the non-polynomial ODE system (10-11) to a system of (K + 2) ODEs with polynomial right-hand sides (Kerner 1981). We obtain for i = 1, 2, . . . , K . The dynamics of the original ODE system (10-11) with the initial condition (x(0), y(0)) = (x 0 , y 0 ) is the same as the dynamics of the extended ODE system (39-41), when we initialize the additional variables by However, when we use a general initial condition, the trajectory of the extended ODE system (39-41) may become unbounded and it may not converge to a limit cycle. To illustrate this behavior, let us consider the initial condition where c > 0 is a constant. If c = 1, then the initial condition (43) reduces to (42). In particular, Fig. 1a shows an illustrative behavior of both the extended ODE system (39-41) for c = 1 and the planar ODE system (10-11), assuming that we use a sufficiently accurate numerical method to solve ODEs (39-41) and plot the projection of the calculated trajectory to the (x, y)-plane. Changing c = 1 to c = 0.5, we plot the dynamics of the extended ODE system in Fig. 2a, where the black dots denote the end points of the calculated trajectories at the final time (t = 100). We observe that only the trajectories which started 'inside a limit cycle' (i.e., their projections to the (x, y)-plane are initially inside a black dashed circle) seem to converge to it, while the other trajectories do not seem to approach the 'limit cycles'. This is indeed the case even if we continue these trajectories for times t > 100. In fact, depending on the accuracy of the numerical method used, all trajectories eventually stop somewhere in the phase plane, because u i (t) → 0 as t → ∞.
On the other hand, considering the extended ODE system (39-41) with the initial condition (43) for c > 1, some additional variables u i (t) tend to infinity as t → ∞, and we again do not observe sustained oscillations in our numerical experiments (results not shown). In particular, the formal conversion of the non-polynomial ODE system (10-11) into the polynomial system (39-41) does not preserve the dynamics well. Therefore, we augment our polynomial ODE system (10-11) in a different way. We introduce K new variables v i , i = 1, 2, . . . , K , and formulate the extended ODE system as the following (K + 2) equations: where ε > 0 is a constant. The first two ODEs (44-45) are the same as ODEs (39-40) with v k taking place of u k . The difference is in the dynamics of the additional variables, i.e., in Eq. (46) which removes the non-polynomial factor (38) in a different way. Rather than defining new variable u i in the form (38) and deriving ODEs which have equivalent dynamics to the ODE system (10-11) for very special initial condition (42), we have written the ODE (46) in such a way that it formally recovers the non-polynomial factor (38) in the limit ε → 0, which will be used in our proof of Lemma 3, where we consider small values of ε. However, even for larger values of ε, the ODE system (44-46) has multiple limit cycles for general initial conditions, as it is illustrated for ε = 1 and K = 4 in Fig. 2b, where all plotted trajectories finish on a limit cycle (see the final calculated positions, at time t = 100, plotted as black dots). Next, we prove that the extended ODE system (44-46) has K limit cycles in the sense of Definition 3 for general values of K . Since (44-46) is a system of (K + 2) ODEs, we cannot directly apply the Poincaré-Bendixson theorem as we did for the planar system in the proof of Lemma 2. While one possible approach to proving the existence of limit cycles is to work with generalizations of the Poincaré-Bendixson theorem to higher-dimensional ODEs (Hirsch 1982;Li and Muldowney 1996;Sanchez 2010), we will base our proof of Lemma 3 on the application of Tikhonov's theorem (Tikhonov 1952;Klonowski 1983) and the result of Lemma 2. In particular, we show that the extended system (44-46) is a polynomial system which has K limit cycles for sufficiently small values of ε provided that the points (a i , b i ) are suffi-ciently separated, as we have previously established in Lemma 2 for the planar ODE system (10-11).
Proof Let us consider ε = 0. Then the right-hand side of the ODE (46) is equal to zero. This equation can be solved for Substituting v i = q i (x, y) into (44-45), we obtain that the reduced problem in the sense of Tikhonov's theorem (Tikhonov 1952;Klonowski 1983) is equal to where functions f (·, ·) and g(·, ·) are defined in (10) and (11). This means that the reduced system (48-49) corresponding to the fast-slow extended ODE system (44-46) is the same as our original non-polynomial ODE system (10-11). Therefore, using Lemma 2, we know that the reduced system (48-49) has (at least) K stable limit cycles in the sense of Definition 3, i.e., there exist K solutions of the reduced system (48-49) which are periodic with period T i > 0 for i = 1, 2, . . . , K . Moreover, there exist ε i > 0, i = 1, 2, . . . , K , such that any solution (x(t), y(t)) of the reduced systems (48-49) approaches the limit cycle (x c,i (t), y c,i (t)) as t → ∞ provided that the initial condition (x(0), y(0)) satisfies Next, we define pairwise disjoint sets i ⊂ R K +2 for i = 1, 2, . . . , K by where functions q j (·, ·) are defined by (47). Let us define Let ε ∈ (0, ε 0 ) be chosen arbitrarily. To show that the extended fast-slow polynomial ODE system (44-46) has (at least) K stable limit cycles, it is sufficient to show that each set i contains one stable limit cycle. We will do this by applying Tikhonov's theorem (Tikhonov 1952;Klonowski 1983). Considering the ODEs (46) for i = 1, 2, . . . , K , where x > 0 and y > 0 are taken as parameters, we obtain the adjoined system as a K -dimensional system of ODEs with an isolated stable equilibrium point [q 1 (x, y), q 2 (x, y), . . . , q K (x, y)], where q i (·, ·) is defined in (47). This equilibrium point attracts the solutions of the adjoined system for any initial condition. Therefore, the ODE system (44-46) has a limit cycle in i . Moreover, this limit cycle attracts any solution x(t), y(t), v 1 (t), v 2 (t), . . . , v K (t) of the ODE system (44-46) with the initial condition satisfying

Chemical Systems with Arbitrary Many Limit Cycles
To construct a CRN with K limit cycles, we first construct a system of ODEs with polynomial right-hand sides which satisfy the condition (9) in Lemma 1, i.e., it will be a system of reaction rate equations, which correspond to a CRN. Once we have such reaction rate equations, there are infinitely many CRNs described by them, so we conclude this section by specifying some illustrative CRNs corresponding to the derived reaction rate equations. Our starting point is the polynomial ODE system (44-46), which has K limit cycles provided that the conditions of Lemma 3 are satisfied. The reaction rate equations are constructed by applying the so-called x-factorable transformation (Plesa et al. 2016) to the right-hand sides of equations (44) and (45). We do not modify the right-hand sides of ODEs (46), because they already satisfy the conditions of Definition 2. We obtain the ODE system: The illustrative dynamics of the ODE system (53-55) is presented in Fig. 3a, where we use the same parameters as we use in Fig. 2b for the ODE system (44-46). We observe that the presented trajectories converge to one of the four limit cycles as in Fig. 2b. The shape of the limit cycles is slightly modified by using the x-factorable transformation, but the limit cycles are still there as we formally prove in Sect. 6. The x-factorable transformations modify the dynamics on the x-axis and y-axis. In Fig. 3a, we present illustrative trajectories which all start with the positive values of x(0) and y(0), while in Fig. 2b, some of the illustrative trajectories have zero initial values of x(0) and y(0). To get a comparable plot, we use the same initial conditions in both Fig. 2b and Fig. 3a, with the only exception that all initial conditions with x(0) = 0 (resp. y(0) = 0) in Fig. 2(b) are replaced by x(0) = 1/2 (resp. y(0) = 1/2) in Fig. 3a. We note that if we start a trajectory of the ODE system (53-55) on the x-axis or the y-axis, then it stays on the axis.
In Fig. 3b, we present illustrative dynamics of the ODE system (53-55) for K = 9, showing that each computed trajectory converges to one of the 9 limit cycles denoted by black dashed lines. To illustrate that this behavior does not require special choices of initial conditions, we used different initial conditions for x(0) and y(0) together with the initial conditions for variables However, a similar figure can be obtained if we replace (56) with the initial condition (43), or if we initialize all values of v i , i = 1, 2, . . . , K as zero (results not shown).
A CRN corresponding to reaction rate equations (53-55) can be obtained (by applying the construction in the proof of Lemma 1) as the CRN with K +2 chemical species, i.e., using the notation in Definition 1, we have To specify the reaction complexes and chemical reactions, we expand the right-hand side of reaction rate equations (53-55). First, we rewrite ODEs (55) as where k i, j , i = 1, 2, . . . , K , j = 1, 2, . . . , 11, are positive constants given by Consequently, the right-hand side of Eq. (55) can be interpreted as the set of 14 chemical reactions for each i = 1, 2, . . . , K . We define it as Consequently, reaction rate equation (55) corresponds to 14 K chemical reactions in sets R i , i = 1, 2, . . . , K . Similarly, we rewrite ODEs (53-54) as where k i, j , i = 1, 2, . . . , K , j = 12, 13, . . . , 21, are constants given by Considering sufficiently large a i and b i (say, for a i > 1 and b i > 1), the constants (63) are positive. Moreover, since the term −v i x 2 y 2 appears in both Eqs. (61) and (62), the right-hand sides of Eqs. (53-54) can be interpreted as the set of 15 K chemical reactions. We define Then, we conclude that the reaction rate Eqs. (53-55) correspond to the CRN with N = K + 2 chemical species and 29 K chemical reactions of at most seventh order given by The CRN (S, C, R) consisting of chemical species S given by (57) and chemical reactions R given by (65) is the CRN which we will use to prove Theorem 1 in Sect. 6. The corresponding set of reaction complexes C can be inferred from the provided lists of reactions R i and R * i , i = 1, 2, . . . , K , given by (60)  The idea of the proof of Theorem 1 is similar to the one chosen in Sects. 3 and 4, where we have first proven Lemma 2 about the existence of K limit cycles in the planar ODE system (10-11) and then we have used it to prove the existence of K limit cycles in the (K + 2)-dimensional ODE system in Lemma 3. In this section, we will again start by formulating Lemma 4 for a planar ODE system, establishing that it has K stable limit cycles. This is again proven by using K disjoint sets (15). The result of Lemma 4 is then used in Lemma 5 to prove Theorem 1 by considering the (K + 2)-dimensional ODE system (53-55). The planar ODE system in Lemma 4 is derived by applying the x-factorable transformation to the planar ODE system (10-11). We obtain where we have used the notation f k (·, ·) and g k (·, ·) introduced in Eqs. (18), (19) and (22). The dynamics of the ODE system (66-67) is similar to the dynamics of the original planar ODE system (10-11) in the same way as the dynamics of the (K + 2)-dimensional extended ODE system (53-55) is similar to the dynamics of the (K + 2)-dimensional extended ODE system (44-46). We have already observed in Fig. 3(a) that the limit cycle around the point (a i , b i ) = (6, 6) of the ODE system (44-46) is relatively circular. On the other hand, the shape of the limit cycles can more significantly differ between Figs. 2b and 3a if the corresponding parameters a i and b i are not equal to each other. Motivated by this observation, we will study the case a i = b i in Lemma 4 and prove that it is possible to choose these parameters in a way that the planar ODE system (66-67) have (at least) K stable limit cycles. This result is sufficient for the proof of Theorem 1. However, we also note that the existence of limit cycles of the ODE system (66-67) is not restricted to the case a i = b i and a more general lemma could be stated and proven, as we did in Lemma 2 where the existence of K limit cycles has been proven under a relatively general condition (14). The advantage of the case a i = b i is that it simplifies the proof of Lemma 4, because we can use the approach and notations introduced in the proof of Lemma 2.
Proof Let us define regions i ⊂ R 2 , i = 1, 2, . . . , K , together with their boundary parts ∂ i1 and ∂ i2 by (15), (16) and (17). Our choice of values of a i and b i in (68) satisfies the assumption (14) in Lemma 2. Therefore, the ODE system (10-11) has with parameters given by (68) at least K stable limit cycles. Moreover, in the proof of Lemma 2, we have shown that each region i does not include any equilibrium of the planar ODE system (10-11). Any equilibrium of the ODE system (66-67) is either located on the x-axis or y-axis, or it is also an equilibrium of the ODE system (10-11). However, our assumption (68) implies that no region i , i = 1, 2, . . . , K , intersects with the x-axis or y-axis. Therefore, we conclude that each i , for i = 1, 2, . . . , K , does not contain any equilibrium of the ODE system (66-67). Next, consider any point (x b , y b ) ∈ ∂ i . We will compute the scalar product of vectors by rewriting the second vector as a sum of two vectors The scalar product of vectors has already been calculated in the proof of Lemma 2 starting with Eq. (33). We obtained that it is negative for ( always points inside the domain i on all parts of the boundary ∂ i . Next, we want to show that this conclusion also holds if vec- as it is done in Eq. (70). To do this, we note that our choice of parameters (68) implies that for all i, j = 1, 2, . . . , K , which not only satisfies the assumption (14), but it can be used in Eq. (37) to make a stronger conclusion that the scalar product of vectors (71) is at most −1.45 for (x b , y b ) ∈ ∂ i1 and at least 1.45 for (x b , y b ) ∈ ∂ i2 . Thus, we only need to show that the scalar product of vectors is in absolute value less than 1.45 to conclude that the original scalar product (69) is Using the definition of g(·, ·) in (22) and the notation z 1 = x b − a i , z 2 = y b − b i introduced in the proof of Lemma 2, we have y b − x b = z 2 − z 1 and the scalar product (72) can be written as Since we have and the second term in (73) is also less than 0.05, we conclude that the scalar product (69) is negative for (x b , y b ) ∈ ∂ i1 and positive for (x b , y b ) ∈ ∂ i2 . Therefore, the vector x b f (x b , y b ), y b g(x b , y b ) always points inside the domain i on all parts of the boundary ∂ i . In particular, applying Poincaré-Bendixson theorem (Strogatz 2015), we conclude that each i contains at least one stable limit cycle. Since i , i = 1, 2, . . . , K , are pairwise disjoint, this implies that the ODE system (66-67) has at least K stable limit cycles.
The following lemma shows that the extended ODE system (53-55) has at least K stable limit cycles when ε is small enough. The detailed proof is omitted since one can use similar steps as in the proof of Lemma 3.
The existence of K limit cycles in the CRN (65) follows by application of Lemma 5.
The right-hand sides of reaction rate equations (75-83) only include quadratic terms. Therefore, there exists a CRN corresponding to the reaction rate Eqs. (75-83) which includes (at most) second-order reactions. We can obtain it by applying the construction in the proof of Lemma 1. The right-hand sides of Eqs. (75) and (76) can be interpreted as the set of 16 K chemical reactions (compare with (64) for ODEs (53-54)) The right-hand side of equations (77) can be interpreted as the set of 14 chemical reactions for each i = 1, 2, . . . , K (compare with (60) for the right-hand side of ODE (55)) Consequently, reaction rate equations (75-77) correspond to 30 K chemical reactions in sets R s, * i and R s i , i = 1, 2, . . . , K . This is already more than 29 K chemical reactions used in Theorem 1, because we did not combine two terms on the right-hand sides into one reaction as we did in the set R * i (this is further discussed in Eq. (92) in Sect. 9). Moreover, there are additional chemical reactions corresponding to the dynamics of additional chemical species in Eqs. (78-83). The right-hand sides of equations (78-81) can be interpreted as the set of 24 chemical reactions given as Finally, the right-hand sides of Eqs. (82-83) can be interpreted as the set of 12 chemical reactions for each i = 1, 2, . . . , K given by In summary, we conclude that the reaction rate equations (75-83) correspond to the CRN with N = 7K + 14 chemical species and 42 K + 24 chemical reactions given by Using Lemma 6, we deduce that the CRN (S, C, R) consisting of chemical species S given by (74) and chemical reactions R given by (89) is an example of a CRN which satisfies Theorem 2. The corresponding set of reaction complexes C can be inferred from the provided lists of reactions R s, * i , R s i , R w and R z i , for i = 1, 2, . . . , K , given by (85), (86), (87) and (88).

Proof of Theorem 3
Theorem 3 Let C(n) be the maximum number of stable limit cycles of reaction rate equations (1-2) corresponding to CRNs with two chemical species undergoing chem-ical reactions of at most n-th order. Then we have where the floor function . denotes the integer part of a positive real number.
Given an arbitrarily large integer K ∈ N, we will show that there exists a CRN with two chemical species such that its reaction rate equations have at least K stable limit cycles and the order of the chemical reactions is at most n(K ) = 6K − 2. To do that, we start with the planar ODE system (10-11) and renormalize time t to get a planar system with polynomial ODEs. Using an auxiliary function we define our new time variable τ by ds.
Then we obtain which is a planar ODE system with its right-hand side given as polynomials of degree n(K ) − 1 = 6K − 3. Since we only rescaled the time, Fig. 1(a) provides an illustrative dynamics of the ODE system (90-91) for K = 4. The illustrative trajectories are calculated in Fig. 1a by solving ODEs (10-11) in time interval t ∈ [0, 100], and we can obtain the same result by solving ODEs (90-91) numerically in time interval τ ∈ [0, 10 −9 ]. Applying x-factorable transformation to ODEs (90-91), we obtain which is a kinetic system of ODEs with polynomials of degree n(K ) = 6K − 2 and which has K stable limit cycles. Solving for K , we obtain K = (n(K ) + 2)/6, which establishes the lower bound (4) in Theorem 3.

Discussion
The main results of this paper have been formulated as Theorems 1, 2 and 3, which show that there exist CRNs with K stable limit cycles for any integer K ∈ N. The CRN presented in our proof of Theorem 1 consisted of N (K ) = K + 2 chemical species S given by (57) and M(K ) = 29 K chemical reactions R (of at most seventh order) given by (65). The number of species and chemical reactions further increases in our proof of Theorem 2, where we restrict our investigation to CRNs with (at most) second-order kinetics. On the other hand, if we restrict to CRNs with only N = 2 chemical species, then the order of the chemical reactions increases with K as n(K ) = 6K − 2 in our proof of Theorem 3. An important question is whether we can further decrease N (K ) (the number of chemical species) and M(K ) (the number of chemical reactions) in Theorems 1 and 2 and still obtain a CRN with K stable limit cycles. One possibility to decrease M(K ) is to use one chemical reaction to interpret multiple terms on the right-hand sides of ODEs (53-55). We have already done this in the reaction set R * i given by (64) with the reaction which corresponds to terms of the form −v i x 2 y 2 appearing in both equations (53) and (54). Another way to construct a CRN with reactions modeling the two terms, −v i x 2 y 2 , in the reaction rate equations (53-54), is to use one chemical reaction per one term on the right-hand side. That is, the chemical reaction (92) could be replaced by two chemical reactions without modifying the form of the reaction rate equations (53-54). In particular, if our aim is to decrease the number M(K ) of chemical reactions, we could consider to 'merge' some other reactions, which have the same reactants. For example, reaction lists (60) and (64) contain chemical reactions If these chemical reactions had the same reaction rate constants k i,7 and k i,17 , then we could replace them by one chemical reaction given by and we would obtain a CRN which has 28 K chemical reactions rather than 29 K , which we use in Theorem 1. Consequently, there is potential to decrease the size of the constructed CRN by a careful choice of our parameters or by modifying the righthand sides of reaction rate equations (53-55). However, the focus of our paper was on the existence proofs and we leave the improvement of bounds on N (K ) and M(K ) to future work. Another possible direction to investigate is to consider more detailed stochastic description of CRNs, written as continuous-time discrete space Markov chains and simulated by the Gillespie algorithm (Erban and Chapman 2020). Such simulations would help us to investigate how our parameters a i , b i , i = 1, 2, . . . , K , need to be chosen that the system not only has the limit cycles of comparable size (as we visualized in Fig. 3 in the ODE setting), but it also follows each of these limit cycles with a similar probability (comparable to 1/K ). This could also be achieved by using the noise-control algorithm (Plesa et al. 2018) for designing CRNs. This algorithm structurally modifies a given CRN under mass-action kinetics, in such a way that (i) controllable state-dependent noise is introduced into the stochastic dynamics, while (ii) the reaction rate equations are preserved. In particular, it could be used to introduce additional chemical reactions (which do not change the ODE dynamics), but lead to controllable noise-induced switching between different limit cycles.

Author Contributions Both authors have contributed equally to this work.
Funding This work was supported by the Engineering and Physical Sciences Research Council, grant number EP/V047469/1, awarded to Radek Erban. This work was also supported by the National Science Foundation, grant number DMS-1620403, and a Visiting Research Fellowship from Merton College, Oxford, awarded to Hye-Won Kang.
Data availability Not applicable. No data are associated with this article.

Declarations
Ethical approval Not applicable.

Conflict of interest Authors have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.