Oscillations in Planar Deficiency-One Mass-Action Systems

Whereas the positive equilibrium of a planar mass-action system with deficiency zero is always globally stable, for deficiency-one networks there are many different scenarios, mainly involving oscillatory behaviour. We present several examples, with centers or multiple limit cycles.


Introduction
In this paper we study mass-action systems in dimension two with a unique positive equilibrium and deficiency one.The deficiency is a non-negative integer associated to any chemical reaction network, and is explained in Section 2. If the deficiency is zero, the Deficiency-Zero Theorem gives a rather complete picture: existence and uniqueness of a positive equilibrium is nicely characterized through the underlying directed graph of the chemical network.And if it exists, it is globally asymptotically stable (at least in dimension two).In contrast, the Deficiency-One Theorem is a purely static statement: there is at most one positive equilibrium, and if it exists, it is regular.However, nothing is said about the dynamic behaviour.In the present paper we take a step towards filling this gap, at least in dimension two.
We show that Andronov-Hopf bifurcations may occur, even generalized ones, that produce more than one limit cycle near the equilibrium, and also degenerate ones that produce centers.Note that the uniqueness and regularity of the positive equilibrium rules out the classical fixed point bifurcations such as saddle-node and pitchfork bifurcations.
We start with a brief summary of reaction network theory in Section 2. Then in Section 3 we study the simplest weakly reversible network with deficiency one, an irreversible cycle along a quadrangle.This is always permanent, but we present examples with up to three limit cycles.On the other hand, we also give a sufficient condition that guarantees global stability of the positive equilibrium for all rate constants.
In Section 4 we consider irreversible chains of three reactions.Even though this is not weakly reversible, a positive equilibrium may exist, and could be globally stable, or may be surrounded by up to three limit cycles.Furthermore, the equilibrium may be surrounded by a continuum of closed orbits and a homoclinic orbit.
In Section 5 we study three separate reactions, and produce two types of centers, and an example with four limit cycles.
Finally, in Section 6 we give a simple example where existence of the positive equilibrium depends on the rate constants.

Planar mass-action systems
In this section we briefly introduce mass-action systems and related notions that are necessary for our exposition.We restrict to the case of two species.For more details about mass-action systems, consult e.g.[7], [9].The symbol R + denotes the set of positive real numbers.
Definition 1 A planar Euclidean embedded graph (or a planar reaction network ) is a directed graph (V, E), where V is a nonempty finite subset of R 2 .
Denote by (a 1 , b 1 ), (a 2 , b 2 ), . . ., (a m , b m ) the elements of V , and by X and Y the two species.Accordingly, we often refer to (a i , b i ) as a i X + b i Y.We assume throughout that the reaction vectors (a j − a i , b j − b i ) ∈ R 2 ((i, j) ∈ E) span R 2 .The concentrations of the species X and Y at time τ are denoted by x(τ ) and y(τ ), respectively.
Definition 2 A planar mass-action system is a triple (V, E, κ), where (V, E) is a reaction network and κ : E → R + is the collection of the rate constants.
We remark that the translation of a network by (α, β) ∈ R 2 (i.e., taking (a i + α, b i + β) instead of (a i , b i ) for i = 1, 2, . . ., m) amounts to multiplying the differential equation (1) by the monomial x α y β , an operation that does not have any effect on the main qualitative properties.Thus, any behaviour shown in this paper can also be realized with a i , b i ≥ 0 for all i = 1, 2, . . ., m, a setting that is more standard in the literature.
In some cases, a network property alone has consequences on the qualitative behaviour of the differential equation (1).For instance, weak reversibility implies permanence [5,Theorem 4.6].We now define these terms.
Definition 3 A planar mass-action system (V, E, κ) is weakly reversible if every edge in E is part of a directed cycle.
Theorem 1 Weakly reversible planar mass-action systems are permanent.
We now recall two classical theorems on the number of positive equilibria for mass-action systems with low deficiency.The deficiency of a planar reaction network (V, E) is the non-negative integer δ = m − − 2, where m = |V | and is the number of connected components of the directed graph (V, E).
(i) There is no periodic solution that lies entirely in R 2 + .(ii) If the underlying network is weakly reversible then there exists a unique positive equilibrium.Furthermore, it is asymptotically stable.(iii) If the underlying network is not weakly reversible then there is no positive equilibrium.
Notice that the combination of Theorems 1 and 2 yields that the unique positive equilibrium of a weakly reversible deficiency-zero planar mass-action system is in fact globally asymptotically stable.
For stating the second classical result, we need one more term.For a directed graph (V, E), denote by t the number of its absorbing strong components.
Theorem 3 (Deficiency-One Theorem [8]) Assume that the deficiency of a planar mass-action system is one.Further, assume that = t = 1.Then the following statements hold.
(i) If the underlying network is weakly reversible then there exists a unique positive equilibrium.(ii) If the underlying network is not weakly reversible then the number of positive equilibria is either 0 or 1. (iii) The determinant of the Jacobian matrix at a positive equilibrium is nonzero.
We now highlight the main differences between the conclusions of the above two theorems.For a planar mass-action system that falls under the assumptions of the Deficiency-One Theorem, (A) in case the underlying network is not weakly reversible, (a) a positive equilibrium can nevertheless exist, (b) whether there exists a positive equilibrium might depend on the specific values of the rate constants, (c) even if there exists a unique positive equilibrium, there could be unbounded solutions as well as solutions that approach the boundary of R 2 + , (B) regardless of weak reversibility, (a) the unique positive equilibrium could be unstable (however, the Jacobian matrix there is guaranteed to be nonsingular), (b) there is no information about the existence of periodic solutions.
Points (a) and (b) in (A) above are studied in detail in [1] and [2], respectively.In this paper, we touch these questions only briefly: we show a network in Section 6 for which the existence of a positive equilibrium is dependent on the specific choice of the rate constants.
Investigation of points (a) and (b) in (B) above is the main motivation for the present paper.The only (published) example so far of a reaction network that satisfies the Deficiency-One Theorem but with an unstable positive equilibrium (and presumably a limit cycle) seems to be the three-species network which is due to Feinberg [8, (4.12)].In this paper we show that such examples are abundant already for two species.Note that Feinberg's example is a bimolecular one.It is shown in [14] that the only bimolecular two species system with periodic solutions is the Lotka reaction [13].
Sections 3 and 4 are devoted to studying the cycle of four irreversible reactions and the chain of three irreversible reactions, respectively.In Section 5 we examine a generalization of the latter: three irreversible reactions that do not necessarily form a chain.The networks in Sections 3, 4, and 6 all satisfy the assumptions of the Deficiency-One Theorem.
Finally, we remark that the Deficiency-One Theorem has a version for the case of more than one connected component (i.e., ≥ 2).However, this does not cover the deficiency-one networks in Section 5. Nevertheless, the conclusions (ii) and (iii) in Theorem 3 hold.

Quadrangle
In this section we study the mass-action system and its associated differential equation under the non-degeneracy assumption that (a 1 , b 1 ), (a 2 , b 2 ), (a 3 , b 3 ), (a 4 , b 4 ) are distinct and do not lie on a line.Note that the mass-action system (2) is weakly reversible and its deficiency is δ = 4 − 1 − 2 = 1.By the Deficiency-One Theorem [8, Theorem 4.2], there exists a unique positive equilibrium for the differential equation (3).Moreover, the determinant of the Jacobian matrix at the equilibrium does not vanish [8,Theorem 4.3].Since additionally the system is permanent [5,Theorem 4.6], the index of the equilibrium is +1 [10,Theorem 19.3].Hence, the determinant is positive, and consequently, the unique positive equilibrium is asymptotically stable (respectively, unstable) if the trace is negative (respectively, positive).In case the trace is zero, the eigenvalues are purely imaginary and some further work is required to decide stability of the equilibrium.
In Section 3.1, we present a system for which the unique positive equilibrium is unstable and a stable limit cycle exists.In Section 3.2, we prove that even three limit cycles are possible for the differential equation (3).Finally, in Section 3.3, we describe a subclass of the quadrangle networks (2) that are globally stable for all rate constants.

Unstable equilibrium and a stable limit cycle
Let us consider the mass-action system (2) with Thus, the network and its associated differential equation take the form A short calculation shows that the unique positive equilibrium is given by and the trace of the Jacobian matrix at the equilibrium is positive if and only if .
By picking rate constants that make the trace positive, one gets a system, where the positive equilibrium is repelling, and, by combining permanence and the Poincaré-Bendixson Theorem, there must exist a stable limit cycle.

Three limit cycles
Let us consider the mass-action system (2) with Thus, the network and its associated differential equation take the form Our goal is to show that there exist rate constants κ 1 , κ 2 , κ 3 , κ 4 such that the above differential equation has three limit cycles.Linear scaling of the differential equation by the equilibrium (x, y), followed by a multiplication by x yields where As a result of the scaling, the positive equilibrium is moved to (1, 1), and the focal value computations become somewhat more convenient.Note that From this, we obtain that for some γ > 0. After dividing by κ 4 , the differential equation ( 4) thus becomes where K > 0 and γ > 0. One finds that the trace of the Jacobian matrix at the equilibrium (1, 1) vanishes for γ = 16 + 1 K .Under this, the first focal value is which is zero for K = K 0 ≈ 0.06862, negative for 0 < K < K 0 , and positive for K > K 0 .Assuming K = K 0 , one finds that the second focal value, L 2 , is approximately 0.01293, a positive number.Take now K = K 0 and γ = 16 + 1 K0 .Since the first nonzero focal value is positive, the equilibrium (1, 1) is repelling.First, perturb K to a slightly smaller value, and simultaneously perturb γ in order to maintain the relation γ = 16 + 1 K .Then L 1 < 0, and thus the equilibrium (1, 1) becomes asymptotically stable, and an unstable limit cycle Γ 1 is created.Next perturb γ to a slightly larger value.Then the trace becomes positive, and thus the equilibrium (1, 1) becomes unstable again, and a stable limit cycle Γ 0 is created.Finally, by the permanence of the system, the Poincaré-Bendixson Theorem guarantees that a stable limit cycle surrounds Γ 1 .Therefore, we have shown that there exist K > 0 and γ > 0 such that the differential equation (5) has at least three limit cycles.
We conclude this subsection by a remark.By keeping b 4 > 2 a parameter (instead of fixing its value to 5), one could find parameter values for which L 1 = 0, L 2 = 0, L 3 < 0 holds (with b 4 ≈ 4.757 and K ≈ 0.0909).Then one can bifurcate three small limit cycles from the equilibrium.

Global stability of the equilibrium
As we have seen in Sections 3.1 and 3.2, the unique positive equilibrium of the differential equation ( 3) could be unstable for some rate constants.However, under a certain condition on the relative position of the four points (a 1 , b 1 ), (a 2 , b 2 ), (a 3 , b 3 ), (a 4 , b 4 ), one can conclude global asymptotic stability of the unique positive equilibrium for all rate constants.
The differential equation ( 3) is permanent and has a unique positive equilibrium.Furthermore, the determinant of the Jacobian matrix is positive there.Hence, by the Poincaré-Bendixson Theorem, global asymptotic stability of the equilibrium is equivalent to the non-existence of a periodic solution.One can preclude the existence of a periodic solution by the Bendixson-Dulac test: if there exists a function h : R 2 + → R + such that div(hf, hg) < 0 then the differential equation ẋ = f (x, y), ẏ = g(x, y) cannot have a periodic solution that lies entirely in R 2 + .With f (x, y) and g(x, y) denoting the r.h.s. of the equations for ẋ and ẏ in (3), respectively, and taking h(x, y) = x −α y −β , one finds where a 5 = a 1 and b 5 = b 1 by convention.Ignoring the degenerate case On the other hand, if a 1 < a 4 < a 2 < a 3 then no matter how one fixes α, at least one of (α − a 2 )(a 2 − a 3 ) and (α − a 4 )(a 4 − a 1 ) is positive.Notice that we covered all configurations with a 1 = min(a 1 , a 2 , a 3 , a 4 ).All the other cases are treated similarly.Also, it works analogously with the b j 's and β.Then there is no periodic solution and the unique positive equilibrium is globally asymptotically stable.
Proof By the above discussion, one can find α and β such that after multiplying by h(x, y) = x −α y −β , the r.h.s. of the differential equation ( 3) has negative divergence everywhere.Then, by the Bendixson-Dulac test, there is no periodic solution and therefore the unique positive equilibrium is globally asymptotically stable.
In other words, if there exists a periodic solution then at least one of a i < a i+3 < a i+1 < a i+2 and b j < b j+3 < b j+1 < b j+2 in Proposition 1 holds.The index-free way to express a i < a i+3 < a i+1 < a i+2 and b j < b j+3 < b j+1 < b j+2 is to say that the projection of the quadrangle to a horizontal line and a vertical line, respectively, take the form , where some arrows are bent in order to avoid overlapping.
Finally, since the mass-action systems in Sections 3.1 and 3.2 have a periodic solution for some rate constants, at least one of a i < a i+3 < a i+1 < a i+2 and b j < b j+3 < b j+1 < b j+2 must hold.Indeed, in each subsection b j < b j+3 < b j+1 < b j+2 holds with j = 2.

Chain of three reactions
In this section we study the mass-action system and its associated differential equation under the non-degeneracy assumption that By the Deficiency-One Theorem, the number of positive equilibria for the differential equation ( 7) is either 0 or 1.Our first goal is to understand when is it 0 and when is it 1.Crucial for this is the relative position of the four points where ∆(ijk) = det(P j − P i , P k − P i ) is twice the signed area of the triangle P i P j P k .The quantity ∆(ijk) is thus positive (respectively, negative) if the sequence P i , P j , P k , P i of points are positively (respectively, negatively) oriented.The quantity ∆(ijk) is zero if the three points P i , P j , P k lie on a line.Note also that Denote by f (x, y) and g(x, y) the r.h.s. of the equations for ẋ and ẏ in (7), respectively.By taking one obtains after a short calculation that the equilibrium equations take the form Thus, if there exists a positive equilibrium, all have the same sign.If all of them are zero then P 1 , P 2 , P 3 , P 4 lie on a line, contradicting the non-degeneracy assumption (8).If the common sign is nonzero then in particular h 4 = −(h 1 + h 2 + h 3 ) = 0, so P 1 , P 2 , P 3 do not lie on a line, and thus the obtained binomial equation ( 9) has exactly one positive solution for each choice of the rate constants.Let us stress that the existence of a positive equilibrium does not depend on the specific choice of the rate constants.Next, we discuss the geometric meaning of Assume that 0, the sequence P 2 , P 3 , P 1 , P 2 is oriented counterclockwise, too.Thus, P 1 and P 4 lie on the same side of the line through P 2 and P 3 (the green open half-plane in the left panel in Figure 1 shows where P 4 can be located for h 1 < 0 to hold).Since additionally h 1 + h 2 < 0 holds, P 1 and P 4 lie on the same side of the line that is through P 3 and is parallel to the line through P 1 and P 2 (the red open half-plane in the left panel in Figure 1 shows where P 4 can be located for h 1 + h 2 < 0 to hold).This latter follows from the fact that h 1 + h 2 = det(P 2 − P 1 , P 4 − P 3 ).In other words, the sum of the angles P 1 P 2 P 3 and P 2 P 3 P 4 is smaller than 180 • (see the two red arcs in right panel in Figure 1).The case h 1 > 0, h 1 + h 2 > 0, h 1 + h 2 + h 3 > 0 is treated similarly, and we obtain the following result.
Proposition 2 Consider the differential equation (7).Then the following four statements are equivalent.
In particular, the existence of a positive equilibrium is independent of the values of κ 1 , κ 2 , κ 3 .
We remark that the equivalence of (a), (b), and (c) in Proposition 2 also follows from [2, Corollaries 4.6 and 4.7], where the existence of a positive equilibrium is discussed for general deficiency-one mass-action systems.Now that we understand when the mass-action system (6) has a positive equilibrium, our next goal is to find parameter values for which the equilibrium is surrounded by three limit cycles (Section 4.1) or by a continuum of closed orbits (Section 4.2).We prepare for these by moving the equilibrium to (1, 1).
•(a1, b 1 ) •(a2, b 2 ) •(a3, b 3 ) half-plane with green vertical lines: location of P 4 for h 1 < 0 to hold half-plane with red horizontal lines: Fig. 1 For a positive equilibrium to exist, the point P 4 is located in the sector that is the intersection of the green and red open half-spaces (left panel).Equivalently, the sum of the two angles indicated is less than 180 • (right panel).
As a result of the scaling, the positive equilibrium is moved to (1, 1).Further, it follows by (9) that for some λ = 0, which is positive (respectively, negative) if h 1 , h 1 + h 2 , h 1 + h 2 + h 3 are all positive (respectively, negative).Denote by J the Jacobian matrix of (10) at the equilibrium (1, 1).A short calculation shows that and thus, det J > 0.
Next, perturb r to a slightly smaller value, and simultaneously perturb K in order to maintain the relation K = 2 r−q(2q+1) .Then L 1 < 0, and thus the equilibrium (1, 1) becomes asymptotically stable, and an unstable limit cycle Γ 1 is created.
Finally, perturb K to a slightly larger value.Then tr J > 0, and thus the equilibrium (1, 1) becomes unstable, and a stable limit cycle Γ 0 is created.
We remark (without proving) that the mass-action systems of this subsection are permanent for all q > 0 and r > 0. In particular, the ones with at least three limit cycles are permanent.

Reversible center
Let us consider now the mass-action system ( 6) with where pq < 0 and p+q = 0. We will prove that the unique positive equilibrium of this mass-action system is a center, provided the rate constants κ 1 , κ 2 , κ 3 are set appropriately.By taking λ = − 1 p 2 −q 2 in (11), we have κ 1 = p−q p , κ 2 = − q p−q , κ 3 = 1, which are indeed all positive under the assumptions on p and q.Setting K = − p q , the associated scaled differential equation (10) then takes the form ẋ = (p − q) + qx p y q − px q y p , ẏ = (q − p) + px p y q − qx q y p .( Proposition 4 The equilibrium (1, 1) is a center of the differential equation ( 14), provided pq < 0 and p + q = 0 hold.
Proof Note that the Jacobian matrix at (1, 1) equals to (p 2 −q 2 ) 0 −1 1 0 .Thus, the eigenvalues are purely imaginary.Since the differential equation ( 14) is of the form ẋ = f (x, y), the system is reversible w.r.t. the line x = y and (1, 1) is indeed a center.
We depicted the typical phase portraits in Figure 2. The one for p + q > 0 suggests that the closed orbits are surrounded by a homoclinic orbit at the origin.In Proposition 5 we show that this is indeed the case.
Fig. 2 The mass-action systems (7) with the substitution (13) (left column), and the phase portraits of the corresponding scaled differential equation ( 14) (right column).The top row is for p + q > 0, while the bottom row is for p + q < 0. Notice that the union of closed orbits is bounded for p + q > 0, and unbounded for p + q < 0.

Three reactions
In this section we study the mass-action system and its associated differential equation under the non-degeneracy assumptions that (a 1 , b 1 ), (a 2 , b 2 ), (a 3 , b 3 ) do not lie on a line, none of (c 1 , d 1 ), (c 2 , d 2 ), (c 3 , d 3 ) equals to (0, 0), and Our first goal is to understand the number of positive equilibria.We find that (x, y) ∈ R 2 + an equilibrium if and only if Notice that, by the non-degeneracy assumptions (17), the three numbers all be zero.Thus, taking also into account that (a 1 , b 1 ), (a 2 , b 2 ), (a 3 , b 3 ) do not lie on a line, the existence of a positive equilibrium is equivalent to Furthermore, once there exists a positive equilibrium, it is unique.Note also that whether there exists a positive equilibrium is independent of the choice of the rate constants κ 1 , κ 2 , κ 3 .Now that we understand when the mass-action system (15) has a positive equilibrium, our next goal is to find parameter values for which the equilibrium is surrounded by four limit cycles (Section 5.1) or by a continuum of closed orbits (Sections 5.2 and 5.3).We remark that the center problem is solved in the special case when one of the reactions is vertical and another one is horizontal [3].Further, the existence of two limit cycles is also discussed there.
We prepare for the rest of this section by moving the equilibrium to (1, 1).Linear scaling of the differential equation ( 16) by the equilibrium (x, y), followed by a multiplication by x yields where As a result of the scaling, the positive equilibrium is moved to (1,1).Further, it follows by (18) that for some λ = 0, which is positive (respectively, negative) if the common sign in ( 19) is positive (respectively, negative).Denote by J the Jacobian matrix of (20) at the equilibrium (1, 1).A short calculation shows that

Four limit cycles
In this subsection we discuss why we strongly conjecture that there exist parameter values for which the differential equation (16) has at least 4 limit cycles.
Let us consider now the mass-action system (15) with , 3 , and it turns out the third focal value, L 3 , is negative at the former point and positive at the latter one.Thus, there exist a, b, d such that L 1 = L 2 = L 3 = 0 (namely, we numerically find that this happens at a = a ≈ 1.01282, b = b ≈ 0.65463, d = d ≈ 3.28862).Since, again numerically, we see that L 4 is negative for a, b, d, we conjecture that there exist parameter values for which the unique positive equilibrium of the differential equation ( 16) is asymptotically stable, and is surrounded by four limit cycles Γ 0 , Γ 1 , Γ 2 , Γ 3 , which are unstable, stable, unstable, stable, respectively.Since the formulas for L 2 , L 3 , and L 4 get complicated, we cannot handle them analytically.This is why we leave the existence of four limit cycles a conjecture.
In any of the above cases, Proposition 6 concludes the proof.
x p y q − 2x q y p ẏ = −1 + 2x p y q − x q y p Fig. 3 Some reaction networks that all fall under Corollary 2, along with the differential equation (25).Among the graphs, in the top row we have p > 0 and −p < q < p, while in the bottom row we have p < 0 and p < q < −p.
We depicted in Figure 3 some reaction networks that all fall under Corollary 2. Now fix p, q, c 1 , c 2 , c 3 , d 1 , d 2 , d 3 such that all the assumptions of Corollary 2 are fulfilled, and consider the mass-action system (16) with (24).How to choose κ 1 , κ 2 , κ 3 in order that the unique positive equilibrium is a center?First note that there exists an (x, y) ∈ R 2 + for which (21) with (24) holds if and only if By (26), κ 1 is arbitrary and κ3 .
Thus, the answer to the above question is that one has to choose κ 1 , κ 2 , κ 3 such that κ 1 is arbitrary and κ Finally, we remark that the condition d1 c1 = d2 c2 d3 c3 in Corollary 2 sheds light on why we had to set (a 4 , b 4 ) = q − p, p + q 2 p in Section 4.2.With this choice, the absolute values of the slopes of the three reactions are q p , 1, q 2 p 2 , the first one being the geometric mean of the latter two.

Liénard center
Let us consider now the mass-action system (15) with Its associated scaled differential equation is then This concludes the proof.positive, and hence λ > 0. Thus, K λ is also positive.Further, under the assumptions of this corollary, it is straightforward to verify the condition (28).Proposition 7 then concludes the proof.
We depicted in Figure 4 some reaction networks that all fall under Corollary 3.

Zigzag
We conclude with an example of a reaction network, where the existence of a positive equilibrium does depend on the choice of the rate constants, a phenomenon that appeared for neither of the networks in Sections 3 to 5.
The mass-action system under investigation in this section, along with its associated differential equation takes the form There is a unique positive equilibrium at 2 − κ for κ < 2 and no positive equilibrium for κ ≥ 2. The determinant of the Jacobian matrix at the equilibrium is positive, while its trace is 5κ − 9, which becomes positive for κ > 9  5 .The Andronov-Hopf bifurcation at κ = 9 5 is subcritical, since L 1 = 5π 13 > 0. The x-axis is invariant, consists of equilibria, and attracts nearby points from R 2 + if κ > 1.For κ > 9  5 it seems that all orbits except the positive equilibrium converge to the x-axis.

2 c 2 2 ẋFig. 4 Corollary 3
Fig. 4 Some reaction networks that all fall under Corollary 3, along with the differential equation (27), and the corresponding phase portraits.