Planar S-systems: Global stability and the center problem

S-systems are simple examples of power-law dynamical systems (polynomial systems with real exponents). For planar S-systems, we study global stability of the unique positive equilibrium and solve the center problem. Further, we construct a planar S-system with two limit cycles.

Andronov-Hopf bifurcations of planar S-systems are discussed in [7], and the first focal value is used to construct a stable limit cycle in [17]. In fact, the first mathematical model of glycolytic oscillations by Selkov [14] is a planar S-system. In previous work, we studied planar S-systems arising from a chemical reaction network (the Lotka reactions) with power-law kinetics. A global stability analysis is given in [1], and the center problem is solved in [2]. (For an interpretation of S-systems as dynamical systems arising from chemical reaction networks with generalized mass-action kinetics [9], we refer the reader to Appendix A.) In this paper, we provide a global stability analysis of arbitrary planar S-systems (Section 3). In particular, we characterize the real exponents for which the unique equilibrium of a planar S-system is globally stable for all positive coefficients. Further, we determine the parameters for which the unique equilibrium is a center (Section 4). In particular, we characterize global centers (Subsection 4.5), and finally we construct a planar S-system with two limit cycles bifurcating from a center (Subsection 4.6). It remains open whether there exist planar S-systems with more than two limit cycles, and we discuss a "fewnomial version" of Hilbert's 16th problem asking for an upper bound on the number of limit cycles for planar power-law systems in terms of the number of monomials. For an illustration of our analysis, we provide figures in Appendix B.
In the following section, we introduce planar S-systems, bring them into exponential form, and discuss the resulting symmetries.
By a nonlinear transformation, we obtain a planar system with the origin as the unique equilibrium. In this exponential form, nullclines are straight lines and symmetries in the exponents can be exploited.

Exponential form
Given the ODE (2), we perform the nonlinear transformations x 1 = e γ1u and x 2 = e γ2v and obtainu = e a1u+b1v − e a2u+b2v , v = e a3u+b3v − e a4u+b4v , defined on R 2 , where The ODE (3) admits the equilibrium (0, 0), and the Jacobian matrix at (0, 0) is given by We abbreviate the ODE (3) by its 8 parameters, more specifically, by the scheme For any a, b ∈ R, the ODE abbreviated by the parameter scheme is obtained by multiplying the vector field in the ODE (3) with e −au−bv and is hence orbitally equivalent to (3). Thus the number of parameters could be reduced from 8 to 6.
Applying a uniform scaling transformation (u, v) → (cu, cv) with c > 0 and rescaling time accordingly is equivalent to dividing all parameters by c. Hence, the parameter space could be reduced to a 5-dimensional compact manifold.

Symmetry operations
In the proofs of our main results (Sections 3 and 4), we exploit symmetries in the parameters, in order to avoid tedious case distinctions.
In fact, the family of ODEs (3) is invariant under the symmetry group of the square (the dihedral group D 4 ) which consists of the following eight elements (rotations and reflections in R 2 ): The question arises how these symmetry operations transform the ODE (3). We start with r 1 , the rotation by 90 So r 1 transforms the ODE (3), abbreviated by the parameter scheme (6), into the ODE abbreviated by The other operations transform the parameter scheme (6) as follows: Note that the symmetry operations r 0 , r 2 , s 0 , s 2 keep the roles of a i and b i (as coefficients of u and v, respectively), whereas the other four operations interchange them. Only the subgroup consisting of r 0 , s 1 keeps the signs of both a i and b i .
3 Global stability Ultimately, we are interested in stability properties of the unique positive equilibrium of the ODE (1). Let G = (g ij ) ∈ R 2×2 and H = (h ij ) ∈ R 2×2 . A short calculation shows that the condition det(G − H) = 0 ensures that the ODE (1) admits a unique positive equilibrium for all given values of the positive parameters α 1 , α 2 , β 1 , β 2 . Then (1, 1) is the unique positive equilibrium of the ODE (2), and (0, 0) is the unique equilibrium of the ODE (3) with (4). On the other hand, if det(G − H) = 0, then the ODE (1) admits either no equilibrium or infinitely many equilibria, depending on the specific values of α 1 , α 2 , β 1 , β 2 .
We call an equilibrium (in R 2 + ) of the ODEs (1) or (2) or an equilibrium (in R 2 ) of the ODE (3) with (4) globally asymptotically stable if it is Lyapunov stable and from each initial condition the solution converges to the equilibrium.
Below, we characterize the parameters G and H (the real exponents) for which the resulting ODEs admit a unique equilibrium that is (globally) asymptotically stable for all other parameters (the positive coefficients). To begin with, we state the obvious relation between the stability properties of the ODEs (1), (2), and (3) with (4).
In our main results, we consider the ODE (3) with (4) and write the Jacobian matrix (5) with (4) as Clearly, sign J = sign(G−H) and sign(det J) = sign(det(G−H)). In particular, det J = 0 if and only if det(G − H) = 0. First, we characterize asymptotic stability.
(ii) det J > 0 and sign J equals one of the sign matrices In particular, these conditions are independent of γ 1 , γ 2 > 0.
Proof. Statement (ii) implies det J > 0 and tr J < 0, and statement (i) follows by a theorem of Lyapunov. By the same theorem (together with the assumption det J = 0), statement (i) implies det J > 0 and tr J ≤ 0 for all γ 1 , γ 2 > 0. The trace condition is equivalent to both diagonal entries of J being non-positive. However, both diagonal entries being zero makes the origin a center, as can be seen using the integrating factor e −a1u−b4v , cf. case S in Subsection 4.1. The signs of the off-diagonal entries follow from det J > 0.
In our main result, we characterize global stability.
The proof of Theorem 3 requires two auxiliary results, Lemmas 4 and 5. There we study the ODE (3) without the substitutions (4). Recall that its Jacobian matrix is given by (5), that is, Lemma 4 on the non-existence of periodic solutions will also be useful in Subsection 4.1, where we look for first integrals of the ODE (3). Lemma 5 on the boundedness of solutions will also be useful in Subsection 4.5, where we solve the global center problem.
Proof. We start by proving (a). In order to prove the boundedness of the solutions in the case a 1 < a 2 , b 3 < b 4 , and det J > 0, we consider all possible signs of a 3 −a 4 and b 1 −b 2 and the corresponding nullcline geometries. For phase portraits in the nine cases, see Figure 1. In two cases (top left and bottom right), solutions may spiral around the origin. Since the divergence of a scaled version of the right-hand side of the ODE (3) is negative (see Lemma 4), they can spiral inwards only (anti-clockwise and clockwise, respectively). In the other seven cases, two of the four regions bounded by the nullclines are forward invariant, hence solutions are ultimately monotonic in both coordinates and converge to the origin.
The symmetry operations (of the square) introduced in Subsection 2.2 preserve det J > 0 and the boundedness of solutions. Hence, statements (c), (d), and (e) follow from (b) by applying the operations s 0 or s 2 , s 1 or s 3 , and r 1 or r 3 , respectively, and it suffices to prove (b).
To prove (b1), first note that a 3 ≤ a 2 follows from a 1 ≤ a 4 by applying the operation r 2 . Since a 2 < a 1 follows from the definition of the sign matrix, it suffices to prove that a 1 ≤ a 4 is necessary for the boundedness. Assume a 1 > a 4 and a 1 > a 2 . A short calculation shows that the set (3) if γ < 0, |γ| is small enough, and u 0 > 0 is large enough. All the solutions starting in this forward invariant set are monotonic in both coordinates and unbounded. For an illustration, see the top panel in Figure 2.
We now show the necessity of a 3 ≤ a 2 = a 1 ≤ a 4 for the boundedness in (b2). The same argument as in the proof of (b1) shows that it suffices to prove that a 1 ≤ a 4 is necessary for the boundedness. Assume a 1 > a 4 and a 1 = a 2 and consider the auxiliary ODEu which can be solved by separation of variables. For v > 0, the curve given by is an orbit of the ODE (15) with u → +∞, v → 0 for t → ∞. All solutions of the ODE (3) that start above this curve are monotonic in both coordinates and unbounded. For an illustration, see the bottom panel in Figure 2.
It is left to show the sufficiency of a 3 ≤ a 2 = a 1 ≤ a 4 for the boundedness in (b2). One can use a Lyapunov function V : Assuming a 3 < a 4 and b 1 > b 2 (recall the assumption on sign J), the boundedness of the sublevel sets of V is equivalent to a 3 ≤ a 1 ≤ a 4 and b 2 ≤ b 4 ≤ b 1 , see Figure 3 for the illustration of the level sets of V . Thus, if in addition to a 3 ≤ a 1 ≤ a 4 , the inequalities b 2 ≤ b 4 ≤ b 1 also hold, the boundedness of the solutions of the ODE (3) follows. In case the inequalities b 2 ≤ b 4 ≤ b 1 do not hold, we also need to take into account the sign structure of the vector field in order to conclude the boundedness of the solutions.
is bounded and forward invariant for all c and for all sufficiently large d.
is bounded and forward invariant for all c and for all sufficiently negative d. For an illustration of the constructed sets, see Figure 4.
Finally, we prove our main result.
Proof of Theorem 3. We have to show that among the systems fulfilling condition (ii) in Proposition 2 exactly those do not admit periodic or unbounded solutions that meet condition (ii) in the present theorem.
In fact, all systems fulfilling condition (ii) in Proposition 2 are covered by Lemma 4 and hence do not admit a periodic solution. Now, statements (a), (b2), (c2), (d2), (e2) in Lemma 5 characterize those systems that do not admit an unbounded solution.

The center problem
An equilibrium is a center if all nearby orbits are closed.
Our aim is to characterize all parameters a 1 , a 2 , a 3 , Additionally, we characterize all the parameters for which the origin is a global center. Finally, we construct a system with two limit cycles.
Let J be the Jacobian matrix of the ODE (3) at the origin. For the origin to be a center, it is a prerequisite that tr J = 0 and det J > 0. (If det J = 0, then the origin lies on a curve of equilibria.) Hence, we assume these conditions throughout Subsections 4.1 and 4.2.

First integrals
We look for first integrals (constants of motion) for the ODE (3) and try an integrating factor of the form e −au−bv . As we have seen in the proof of Lemma 4, the divergence is proportional to First, we consider a 1 = a 2 and b 3 = b 4 . (case S) Setting a = a 1 and b = b 4 , all four terms vanish, and the system is integrable. In fact, the ODE (3) is orbitally equivalent tȯ a system with separated variables. This case has codimension 2 in the parameter space.
Next, we consider The divergence (17) simplifies to Setting a = a 2 and b = b 4 , the last two terms vanish, and the first term is zero due to tr J = 0. This case has codimension 3 in the parameter space.
The following cases can be treated in the same way (and have codimension 3 in the parameter space): It is easy to see that case S is invariant under all symmetry operations.
On the other hand, cases I1-I4 can be obtained from each other by symmetry operations. Below, we apply all symmetry operations to case I1: The corresponding first integrals are displayed in Table 1.

Reversible systems
Let R : R 2 → R 2 be a reflection along a line. A vector field F : R 2 → R 2 (and the resulting dynamical system) is called reversible w. : an equilibrium of a reversible system which has purely imaginary eigenvalues and lies on the symmetry line of R is a center.
The above definition can be generalized and the fact still holds: A vector field  (14) to (11) is equivalent to the original scheme, This holds if and only if there exist a, b ∈ R such that or, equivalently, The ODE (3) is reversible w.r.t. s 3 , the reflection along the line u = −v, if the system transformed by s 3 followed by time reversal is orbitally equivalent to the original system. That is, if applying (14) to (13) is equivalent to the original scheme, This holds if and only if there exist a, b ∈ R such that that is, or, equivalently, The two families of reversible systems given by (18) and (19), respectively, have codimension 3 in the parameter space. The other two reflections, s 0 and s 2 (across the u-and v-axis), also lead to reversible systems, however, they are already covered by case S.
Cases R1 and R2 can be obtained from each other by symmetry operations. Below, we apply all symmetry operations to case R1: R1 r0,r2,s1,s3 r1,r3,s0,s2 / / R2 Finally, we remark that neither for R1 nor for R2 we were able to find a first integral. However, for systems that are in the intersection of R1 and R2, the functions 1 + e r(u+v) [e qu + e qv ] − r q and e −a1u−b4v (e qu + e qv ) − q+r q serve as first integral and integrating factor, respectively, where q = a 4 − a 1 and r = a 3 − a 1 .

Main result
In Table 2, we display the seven cases of centers we identified in Subsections 4.1 and 4.2. Indeed these are all possible centers of the ODE (3).
Theorem 6. Let J be the Jacobian matrix of the ODE (3) at the origin, that is, The following statements are equivalent: 1. The origin is a center of the ODE (3).
2. The eigenvalues of the Jacobian matrix at the origin are purely imaginary, that is, tr J = 0 and det J > 0, and the first two focal values vanish.
Proof. 1 ⇒ 2: If J has a zero eigenvalue, that is, det J = 0, then the origin lies on a curve of equilibria and cannot be a center. Hence, the eigenvalues of J are purely imaginary, and all focal values vanish.
After performing the substitutions the ODE (20) is orbitally equivalent to (3). Using we compute det J and the first two focal values, L 1 and L 2 . We find and note that det J > 0 implies a 3 = a 4 and b 2 = 0. Further, using the Maple program in [6], we find Expressions for L 2 (in case L 1 = 0) will be given below.
We show that all parameters (a 2 , b 2 , a 3 , b 3 , a 4 , b 4 ) ∈ R 6 in the ODE (20) for which tr J = L 1 = L 2 = 0 and det J > 0 belong to one of the seven cases in Table 3.
To begin with, L 1 = 0 implies either Case (a), where b 3 = b 4 and a 2 = 0 (due to tr J = 0), corresponds to case S in Table 3.
In case (b), where D = 0 (and b 3 , b 4 = 0 due to b 2 = 0), we find using the Maple program in [6]. Now, L 2 = 0 implies that at least one of six factors is zero: -As shown above, the first subcase b 3 − b 4 = 0 is covered by case S in Table 3.
-The subcase a 4 + b 4 = 0 implies D = −a 4 b 3 = b 3 b 4 and hence b 2 = a 3 − a 4 . Adding a 2 = b 3 −b 4 (due to tr J = 0) yields a 2 +b 2 = a 3 +b 3 , and the situation is covered by case R1. , and hence which need not be considered further.
In case (c2), where a 4 + b 4 = 0 and b 3 = 0, we find Now, L 2 = 0 implies that at least one of five factors is zero: -The first subcase a 3 = 0 (and b 3 = 0) is covered by case I1 in Table 3.
-As mentioned above, the subcase a 3 − a 4 = 0 implies det J ≤ 0 which need not be considered further.
To obtain the case distinction for the ODE (3), we perform the substitutions (21) in Table 3. The result is displayed in Table 2.

Global centers
We say that the origin is a global center of the ODE (3) if all orbits are closed and surround the origin.
Proof. For the cases S, I1, I2, I3, I4, the theorem follows immediately by investigating the level sets of the first integrals, see Table 1.
Below, we will implicitly use the easily checkable fact that condition (23) is invariant under any of the symmetry operations r 0 , r 1 , r 2 , r 3 , s 0 , s 1 , s 2 , s 3 .
It suffices to show the theorem for the case R1, because the case R2 then follows by applying any of the symmetry operations r 1 , r 3 , s 0 , s 2 . In the sequel, we consider only R1. Also, we can assume that the system under consideration is not in case S, and thus sign J is one of The 1st and the 3rd of these four cases can be transformed to each other by s 1 and s 3 . The same applies to the 2nd and the 4th. Thus, we restrict our attention to the cases Another short calculation shows that the two chains of inequalities in (23) are equivalent for R1. Note also that in the case R1 the ODE (3) can be written in the orbitally equivalent forṁ v = e a3u+a2v − e a4u+a1v .
Therefore, we have to show that (i) if a 1 > a 2 and a 3 < a 4 then the origin is a global center for the ODE (24) if and only if a 3 ≤ a 2 ≤ a 1 ≤ a 4 and (ii) if a 1 < a 2 and a 3 < a 4 then the origin is a global center for the ODE (24) if and only if a 3 ≤ a 1 ≤ a 2 ≤ a 4 .
In both of the cases (i) and (ii), the necessity of the inequality chain between a 1 , a 2 , a 3 , a 4 follows immediately from Lemma 5.
The sufficiency in the case (i) follows directly by taking into account the nullcline geometry, the sign structure of the vector field, the fact that all the orbits are symmetric w.r.t the u = v line, and the easily checkable fact thatu +v < 0 whenever u > v, whileu +v > 0 whenever u < v, see the top panel in Figure 5.
(The sign ofu +v equals the sign of v − u, because both of the differences e a1u+a4v − e a4u+a1v and e a3u+a2v − e a2u+a3v are nonpositive (respectively, nonnegative) for u > v (respectively, for u < v) andu +v can be zero only if u = v, because a 1 = a 4 and a 2 = a 3 would imply det J = 0.) In case (ii), we consideru −v instead ofu +v. We cannot determine where exactly it is positive and negative. However, it is enough that we know that it is negative (respectively, positive) whenever both of (a 1 − a 3 )u + (a 4 − a 2 )v and (a 4 − a 2 )u + (a 1 − a 3 )v are negative (respectively, positive), see the bottom panel in Figure 5. Starting from an initial point withu < 0 andv > 0, the solution will cross the u-nullcline and enter the region, whereu > 0 andv > 0.
Then the solution will reach the region, whereu −v > 0. Afterwards, it hits the v-nullcline and then the u-nullcline, after whichu < 0 and therefore the solution will reach the region, whereu −v < 0. From there, it will hit the v-nullcline again.
Finally, we remark that the center is clockwise (respectively, anticlockwise) if and only if a 3 < a 4 and b 1 > b 2 (respectively, a 3 > a 4 and b 1 < b 2 ).

Limit cycles
For the ODE (3), we are also interested in asymptotic stability when the trace of the Jacobian matrix vanishes, that is, when linearization does not give any information. In fact, using the (sign of the) first focal value computed in Subsection 4.4, we characterize super-and subcritical Hopf bifurcations resulting in a stable or unstable limit cycle, see also [7].
Proposition 8. For the ODE (3), let det J > 0 and tr J = 0 at the origin and If 1 < 0, the origin is asymptotically stable. If 1 > 0, it is repelling.
If we consider a one-parameter family of ODEs (3) along which the eigenvalues of the Jacobian matrix cross the imaginary axis with positive speed, for example, with parameter µ = tr J, then an Andronov-Hopf bifurcation occurs at µ = 0. If 1 < 0, the bifurcation is supercritical (and there exists an asymptotically stable closed orbit for small µ > 0). If 1 > 0, it is subcritical (and there exists a repelling closed orbit for small µ < 0).
Further, we are interested in a degenerate Hopf or Bautin bifurcation resulting in two limit cycles, see [5,Section 8.3]. Indeed, using the first two focal values computed in Subsection 4.4, we construct an S-system with two limit cycles.
In particular, we consider case (c4) in Subsection 4.4: we set a 1 = b 1 = b 4 = 0, a 3 = b 3 = a 2 and hence tr J = L 1 = 0 and choose a 2 , b 2 , a 4 such that L 2 < 0 (and det J > 0) with L 2 given by Equation (22), for example, a 2 = −1, b 2 = −2, a 4 = 4. By slightly decreasing b 3 and a 2 (thereby keeping b 3 = a 2 and tr J = 0), we obtain L 1 > 0, and the resulting system has a stable limit cycle. Finally, by slightly increasing a 2 such that tr J < 0, we create a small unstable limit cycle via a subcritical Hopf bifurcation.
It remains open, whether the ODE (3) admits more than two limit cycles. In fact, one could formulate a "fewnomial version" of the second part of Hilbert's 16th problem for planar power-law systems defined on the positive quadrant: Khovanskii [4] gives an explicit upper bound on the number of nondegenerate positive solutions of n generalized polynomial equations in n variables in terms of the number of distinct monomials; see also [15]. Similarly, we can ask for an upper bound on the number of limit cycles of planar power-law systems (with finitely many equilibria) in terms of the number of monomials.
In analogy to the cyclicity problem (the local version of Hilbert's 16th problem), we can also ask for an upper bound on the number of limit cycles that can bifurcate from a center, when we fix the number of monomials and their signs and perturb the positive coefficients and real exponents. Our example shows that in the simplest case with two binomials this upper bound is at least two. For a computational algebra approach to this question for planar polynomial systems with real or complex coefficients and integer exponents of small degree, see [11].
containing two connected components with two vertices and two edges each, To each vertex, one assigns a stoichiometric complex (either the zero complex 0 or one of the molecular species X 1 and X 2 ), in particular, one specifies the reversible reactions representing the production and consumption of X 1 and X 2 .
For mass-action systems, the deficiency (a nonnegative integer) plays a crucial role in the analysis of the dynamical behaviour. For example, if the deficiency is zero, then periodic solutions are not possible. For the generalized mass-action system (25), the stoichiometric deficiency [9] is given by δ = 4 − 2 − 2 = 0, since there are 4 vertices and 2 connected components in the graph and the stoichiometric subspace has dimension 2. In contrast to mass-action systems with deficiency zero, this system gives rise to rich dynamical behaviour.
Analogously, every n-dimensional S-system can be specified as a generalized mass-action system in terms of [9] with deficiency zero. In fact, every generalized mass-action (GMA) system in terms of biochemical systems theory (BST) can be specified as a generalized mass-action system in terms of [9]. More specifically, every power-law dynamical system arises from a generalized chemical reaction network, that is, a digraph without self-loops and two functions assigning to each vertex a stoichiometric complex and to each source vertex a kinetic-order complex. Thereby, complexes need not be different, as in the case of the zero complex 0 in the generalized mass-action system (25).     : Illustration of the proof of Theorem 7, case R1, to show the sufficiency of a 3 ≤ a 2 < a 1 ≤ a 4 (and a 3 ≤ a 1 < a 2 ≤ a 4 , respectively) for the origin being a global center of the ODE (3). Both panels display the nullcline geometry, the sign structure of the vector field, the line of reflection, and the signs ofu +v andu −v.