Dynamics of the Higgins–Selkov and Selkov systems

. We describe the global dynamics in the Poincar´e disc of the Higgins–Selkov model where k 2 are positive parameters, and of the Selkov model where a,b are positive parameters.


Introduction and statement of the results
The Higgins-Selkov model of glycolysis is where the unknowns x and y are concentrations which are non-negative and k i for i = 0, 1, 2 are the reaction positive constants, see [8] for the biological details of this model.
We will describe the global dynamics of the differential system (1) in the Poincaré disc for all positive values of k 0 , k 1 and k 2 . For a definition of the Poincaré disc, and of its separatrices and canonical regions see subsection 3.6 of the appendix. We denote by S (respectively R) the number of separatrices (respectively canonical regions) of a phase portrait in the Poincaré disc. Thus our first main result is: Theorem 1. The Higgins-Selkov system (1), after a rescaling of its variables, can be written as (2) x ′ = 1 − xy 2 , y ′ = ay(xy − 1), with a > 0. The global phase portraits of this system for a ∈ R is topologically equivalent to the one of shaded areas correspond to the initial conditions of the orbits having a finite final evolution, so these are the initial conditions with biological meaning. In the phase portrait C the final behavior is an equilibrium point, while in the phase portrait D is a stable limit cycle. The phase portraits B and E are non-generic because they only occur for a fixed value of the parameter. Note that the final behaviors of the orbits of the phase portraits A and F take infinite values.
• For system (2) the bifurcation values of the parameter a are a = 0, a = 1 and a = a * . We note that the Higgins-Selkov system (1) only reflects some biological meaning for some initial conditions when a ∈ (0, a * ), otherwise the orbits that do not go to infinity have zero Lebesgue measure. In Figure  1(C) and (D) the initial conditions with some biological meaning are in the shaded areas.
The proof of Theorem 1 modulo a conjecture is given in section 2. In that section the conjecture is well stated. That conjecture is supported by numerical computations.
The Selkov model of glycolysis is given by the differential system with a and b positive parameters. The parameter b is named phosphofructokinase and the parameter a is called hexokinase which is the activant from all the glycolytic cycle. For a detailed derivation of system (3) see [8]. Our main aim is to study the global phase portraits on the Poincaré disc for the system (3) in function of its parameters. This is the content of the second main result in the paper.
Theorem 2. The Selkov system (3), after a rescaling of the variables, can be written as with a > 0. The global phase portraits of this system for a ∈ R is topologically equivalent to the one of For system (2) the bifurcation values of the parameter a are a = −1, a = a 1 and a = 0. We note that the Selkov system (4) only reflects biological Figure 2. Phase portraits of the Selkov model. Only the phase portrait G, which corresponds to a > 0 for system (4), has biological meaning. As before the shaded areas correspond to the initial conditions of the orbits having a finite final evolution.
meaning when a > 0, then the unique finite equilibrium point is a global attractor.
The proof of Theorem 2 modulo a conjecture is given in section 3, where the conjecture is well stated. Again that conjecture is supported by numerical computations.
We must note that both conjectures are on the limit cycles of the Higgins-Selkov and Selkov models. These conjectures state that in both systems the number of limit cycles is at most one, and its existence depends on the values of the parameters. This limit cycle appears from a Hopf bifurcation (this is proved) and disappears in a graphic having a part at infinity (there is only numerical evidence of this last fact). As usual to control the existence or non-existence of limit cycles is one of the main problems of the qualitative theory of the differential equations in dimension two.

Proof of Theorem 1
Since the parameters k i for i = 1, 2, 3 are positive, we can do the rescaling of the variables and system (1) becomes system (2) where a = k 3 2 /(k 2 0 k 1 ) > 0, the prime denotes derivative with respect to the new time T and we have written x and y instead of X and Y . For completeness we shall study the phase portraits of system (2) for all a ∈ R, but for their applications to biology only a > 0 needs to be considered.

2.1.
Finite singular points and the Hopf bifurcation. System (2) has one finite singular point, namely the point (1, 1). See subsection 3.5 of the appendix for the definitions of the different kind of singular points, and Chapter 5 of [2], for the notions of weak focus and Lyapunov constant. Computing the eigenvalues of the Jacobian matrix at this singular point we get that they are (a − 1 ± √ 1 − 6a + a 2 )/2. The fixed point (1, 1) is It has a Hopf-bifurcation at a = 1 where the stability of the focus (1, 1) changes, for more details on Hopf-bifurcations see [3]. For a > 1 sufficiently DYNAMICS OF THE HIGGINS-SELKOV AND SELKOV SYSTEMS 5 small a stable limit cycle comes from the Hopf bifurcation (see Figure 1(D)) that ends in a graphic with some part at infinity for some value of a = a * ∈ (1.23, 1.24) (see Figure 1(E)), this value has been obtained by numerical integration of the differential system (2).

2.2.
Infinite singular points. Now we shall study the infinite singularities of system (2). In the local chart U 1 system (2) becomes Consider the case a = 0. The infinite singular points in the local chart U 1 are q 0 = (0, 0) and q 1 = (−a, 0). Computing the Jacobian matrix at the singular point q 0 we get that it is identically zero. Using the blow up technique (see subsection 3.5 of the appendix) we find that it is formed by two hyperbolic sectors and two parabolic ones, the straight line of infinity separates the hyperbolic sectors and the parabolic ones as it is described in the figures (C), (D), (E) and (F) of Figure 1. On the other hand, the eigenvalue of the Jacobian matrix at the singular point q 1 is a 2 with multiplicity two, so it is an unstable hyperbolic node.
If a = 0 then the straight lines y =constant are invariant and the curve 1 − xy 2 = 0 if filled of singular points, consequently the phase portrait of system (2) is given in Figure 1(B).
In the local chart U 2 system (2) becomes The origin of the local chart U 2 is an infinite singular point which is semihyperbolic (see subsection 3.5 of the appendix). When a = 0 we get that the origin of U 2 is a stable node if a < 0, and a topological saddle if a > 0, see Figure 1.

2.3.
On the periodic orbits. It is well known that if in the region bounded by a periodic orbit there are finitely many singular points, then the sum of the topological indices of these singular points is equal to one, see for instance Proposition 6.26 of [2]. Therefore, since if a < 0 the unique finite singular point is a saddle (topological index −1 see for instance Proposition 6.32 of [2]), then the differential system (2) has no periodic orbits if a < 0. Proof. Note that x ′ |x=0 = 1 > 0 and y ′ |y=0 = 0 and if there exists a periodic solution, since it must surround the singular point (1, 1), we get that if the limit cycle exists must be contained in the first quadrant {x > 0, y > 0}. The divergence of system (2) is We claim that the curve f (x, y) ∩ {x > 0, y > 0} is transversal for the flow of system (2) if a > 3 and transversal except on the point (2/3, 1) (in which is tangent) for a = 3. This claim together with the fact that the periodic orbit must surround the singular point (1, 1) will imply that if the periodic orbit exists must be contained in a simple connected region where the divergence is positive (or positive except in the point (2/3, 1) where it is zero) and this is not possible by the Bendixson criterium (see for instance Theorem 7.10 of [2]). Now we prove the claim. Solving f (x, y) = 0 we get Take the curve and consider the function Note that p a (y) is a fourth degree polynomial that we write as a 4 + a 3 y + a 2 y 2 + a 1 y 3 + a 0 y 4 with a 0 = −3, a 1 = a 2 = 0, a 3 = 4a and a 4 = −a 2 . Computing − 4a 0 a 3 2 a 2 3 + 18a 2 a 4 a 3 1 a 3 + 144a 0 a 2 a 2 4 a 2 1 − 80a 0 a 2 2 a 4 a 1 a 3 + 18a 0 a 2 a 3 3 a 1 − 4a 3 2 a 4 a 2 1 − 4a 3 1 a 3 3 + 16a 0 a 4 2 a 4 − 128a 2 0 a 2 2 a 2 4 + 144a 2 0 a 2 a 4 a 2 3 = 6912(a − 3)a 4 (a + 3), we get that D 4 > 0 for a > 3 and D 4 = 0 for a = 3. It follows from [7] that in the first case (i.e when a > 3), the polynomial p a (y) has four complex roots and so T a (x, y) > 0 on {x > 0, y > 0}; and for a = 3 it has two complex solutions and the double real solution y = 1, and so T 3 (x, y) ≥ 0, and only is zero at the point (2/3, 1) of (x, y) ∈ {x > 0, y > 0}. This proves the claim and concludes the proof of the lemma.
In short we have proved the following result. We have numerical evidence of the following conjecture: Conjecture 1. The differential system (2) has a unique periodic solution when a ∈ (1, a * ) where a * ∈ (1.23, 1.24). Moreover this periodic solution is a stable limit cycle which borns in a Hopf bifurcation at the singular point (1, 1) when a = 1 and ends in a graphic with a part at infinity when a = a * .

2.4.
The phase portraits in the Poincaré disc. Gluing all the finite and infinite information of system (2) together with the conjecture and the fact that the straight line y = 0 is invariant, we get the six phase portraits of Figure 1 described in the statement of Theorem 1. In particular, note that the bifurcation values of a are 0, 1 and a * , and that the biological meaning of the Higgins-Selkov model only takes place for the values of the parameter a ∈ (0, a * ], otherwise the orbits of system (2) that do not go to infinity have zero Lebesgue measure. This concludes the proof of Theorem 1.

Proof of Theorem 2
First we compute the normal form given in Theorem 2. Since b = 0 doing the rescaling of coordinates system (3) becomes the differential system (4) with a/b 2 replaced by a, writing (x, y) instead of (X, Y ), and where the prime denotes derivative with respect to the new time T . Now we shall study the dynamics of the differential system (4).

3.1.
Finite singular points and the Hopf bifurcation. System (4) has a unique finite singular point (1, 1/(1 + a)) when a = −1. Computing the eigenvalues of the Jacobian matrix at this point we get that it is • a saddle if a ∈ (−∞, −1), • an unstable hyperbolic node if a ∈ (−1, a 2 ] with • an unstable hyperbolic focus if a ∈ (a 2 , 0), • an unstable weak focus if a = 0 because the first Lyapunov constant is 1/8, • a stable hyperbolic focus if a ∈ (0, a 3 ) with
From the above it follows that system (4) has a Hopf bifurcation at a = 0 where the focus changes its kind of stability (see Figure 2(F)), from it bifurcates a stable limit cycle for a < 0 sufficiently small, see for more details on this Hopf bifurcation [3]. Doing numerical computations this stable limit cycle ends in a graphic with a part at infinity for the value a = a 1 ∈ (−0.036, −0.037), see Figure 2(E).

3.2.
Infinite singular points. In the local chart U 1 system (4) becomes The only infinite singular points in the local chart U 1 are (0, 0) and (−1, 0). Computing the eigenvalues of the Jacobian matrix at the origin we get that this singularity is semi-hyperbolic. Using Theorem 2.19 in [2] we conclude that it is a saddle. On the other hand, the eigenvalues of the Jacobian matrix at the point (−1, 0) are 1, 1 and so (−1, 0) is an unstable hyperbolic node.
In the local chart U 2 system (4) becomes The origin of the local chart U 2 is a linearly zero singular point. We apply the blow-up techniques to study it. We do a horizontal blow up. We consider the new variables (u, w) where w = v/u. In these new variables (6) can be writtenu = u 2 (1 + u + aw 2 + (a − 1)uw 2 − u 2 w 3 ), w = −uw(1 + aw 2 − uw 2 ). We eliminate the common factor u by making a rescaling of time and we get System (7) when a ≥ 0 has a unique singular point on u = 0 which is (u, w) = (0, 0). Computing the eigenvalues of the Jacobian matrix of system (7) at the origin we get that they are 1 and −1. Hence, the origin is a saddle. Doing now the blowing down, and passing from (u, w) to (u, v) we get that the origin of U 2 is the union of two hyperbolic sectors (see Figure 2(G)).
On the other hand if a ≤ 0 system (7) has three singular points on u = 0, which are (u, w) = (0, 0) and (u, w) = (0, ±1/ √ −a). Computing the eigenvalues of the Jacobian matrix at the point (u, w) = (0, 0) we get that they are 1, −1 and so it is a saddle. On the other hand computing the eigenvalues of the Jacobian matrix at the singular points (u, w) = (0, ±1/ √ −a) we get that they are 2 and 0. Hence these two singularities are semi-hyperbolic. Using Theorem 2.19 of [2] we get that are saddles if a < −1, Doing now the blowing down, passing from (u, w) to (u, v) we get that the origin of U 2 is formed by • two elliptic and two parabolic sectors (topological index 2) if a < −1, see Figure 2(A); • one elliptic, one hyperbolic and one parabolic sectors (topological index 1) if a = −1, see Figure 2(B); • two hyperbolic and one parabolic sectors (topological index 0) if a ∈ (−1, 0], see Figure 2(C)-(F).

3.3.
On the periodic orbits. Since if a < −1 the unique finite singular point is a saddle (topological index −1), then the differential system (4) has no periodic orbits if a < −1.
Since the the divergence of system (4) is negative in the strip R for a > a, and for a = a is negative except at the point (1/a, 2/a 2 ), by the Bendixson criterion, no periodic solutions can exist in R. This concludes the proof of the lemma.
We have numerical evidence of the following conjecture: Conjecture 2. The differential system (4) has a unique periodic solution when a ∈ (a 1 , 0) where a 1 ∈ (−0.036, −0.037). Moreover this periodic solution is a stable limit cycle which borns in a graphic with a part at infinity when a = a 1 and ends in a Hopf bifurcation at the singular point (1, 1/(1+a)) when a = 0.
3.4. The phase portraits in the Poincaré disc. Taking into account all the information on the finite and infinite singular points, the Hopf bifurcation, and the Conjecture 2 we conclude that the global phase portraits for the differential system (4) are the seven phase portraits of Figure 2 described in the statement of Theorem 2. In particular, note that the bifurcation values of a are −1, a 1 and 0, and that the biological meaning of the Selkov model is given when a > 0 by the phase portrait of Figure 2(G). Hence, when a > 0 the unique finite singular point is a global attractor. This concludes the proof of Theorem 2.

Appendix: Basic results
In this appendix we summarize some basic notations and results which are necessary for stating and proving the results presented in this paper.
3.5. Singular points of differential systems in plane. Letẋ = P (x, y), y = Q(x, y) be a differential system in the plane R 2 . An equilibrium point or a singular point of this differential system is a point (x 0 , y 0 ) ∈ R 2 such that P (x 0 , y 0 ) = Q(x 0 , y 0 ) = 0.
A singular point (x 0 , y 0 ) is hyperbolic if the eigenvalues of the Jacobian matrix have nonzero real part.
A singular point is semi-hyperbolic if one of the eigenvalues of the matrix (8) is zero.
A singular point is nilpotent if the eigenvalues of the matrix (8) are both zero but the matrix is not identically zero.
Finally, a singular point is linearly zero if the matrix (8) is identically zero.
The local phase portraits of the hyperbolic, semi-hyperbolic and nilpotent singular points are well-known see for instance the Theorems 2.15, 2.19 and 3.5 of [2]. While the study of the local phase portraits of the linearly zero singular points is done using special change of variables called blow-ups, see for instance [1].
3.6. Poincaré compactification. Let P (x 1 , x 2 ) and Q(x 1 , x 2 ) be real polynomials in the variables x 1 and x 2 , and let X = (P, Q) be a polynomial vector field of degree d in R n , i.e. d the maximum of the degrees of the polynomials P and Q. For all the details on the Poincaré compactification see Chapter 5 of [2].
The Poincaré sphere is S 2 = {y = (y 1 , y 2 , y 3 ) ∈ R 3 : 3 i=1 y 2 i = 1}. We denote by T y S 2 the tangent space at the point y of S 2 . We identify the tangent space T (0,0,1) S 2 with the plane R 2 where it is defined the polynomial vector field X. Consider the central projection f : T (0,0,1) S 2 −→ S 2 defined as follows: if q ∈ T (0,0,1) S 2 then f (q) is formed by the two intersection points of the straight line which connects the point q with the origin of coordinates with the sphere S 2 . Under the central projection f the infinity of the plane R 2 ≡ T (0,0,1) S 2 goes to the equator S 1 = {y ∈ S 2 : y 3 = 0} of S 2 .
Through the central projection f we obtain two copies Df • X of the polynomial vector field X on the sphere S 2 , one in the southern hemisphere and the other in the northern one. Let X ′ be the vector field on the sphere S 2 minus its equator S 1 formed by these two copies of X. The vector field X ′ can be extended from S 2 \S 1 to an analytic vector field p(X) on S 2 taking p(X) = y d+1 3 X ′ . Consider the projection π : R 3 −→ R 2 defined by π(y 1 , y 2 , y 3 ) = (y 1 , y 2 ). The projection π sends the closed northern hemisphere of S 2 into the Poincaré disc D. Thus the interior of D is diffeomorphic to R 2 and its boundary S 1 corresponds to the infinity of R 2 , and π(p(X)) is the extension of the polynomial vector field X from R 2 to the Poincaré disc D, called the Poincaré compactification of the polynomial vector field X.
In order to have the explicit expression of the Poincaré compactification p(X) for i = 1, 2, 3 we take local charts (U i , F i ) and (V i , G i ) on the sphere S 2 , defined by U i = {y ∈ S 2 : y i > 0}, V i = {y ∈ S 2 : y i < 0}, F i : U i → R 2 and G i : with 1 ≤ j 1 < j 2 ≤ 3 and j k = i for k = 1, 2. The expression of p(X) in the local chart (U 1 , F 1 ) is z d 2 (−z 1 P + Q, −z 2 P ), where P = P (1/z 2 , z 1 /z 2 ) and Q = Q(1/z 2 , z 1 /z 2 ), after a convenient rescaling in the independent variable.
Similarly the expression of p(X) in the local chart (U 2 , F 2 ) is z d 2 (−z 1 Q + P, −z 2 Q), where P = P (z 1 /z 2 , 1/z 2 ) and Q = Q(z 1 /z 2 , 1/z 2 ). The expression of p(X) in the chart (V i , G i ) is the same than in the chart (U i , F i ) multiplied by (−1) d for i = 1, 2.
An equilibrium point of π(p(X)) or X is called finite (respectively infinite) if it belongs to the interior (respectively to the boundary) of the Poincaré disc.
3.7. Topological equivalent polynomial vector fields. Two polynomial vector fields X and Y on R 2 are topologically equivalent if there exists a homeomorphism on the Poincaré disc D preserving the infinity S 1 and carrying orbits of the vector field π(p(X)) into orbits of the vector field π(p(Y )), either preserving or reversing the sense of all the orbits.
The separatrices of the Poincaré compactification π(p(X)) are the equilibrium points, the limit cycles, the orbits of the boundary of the hyperbolic sectors at a finite or infinite equilibrium point, and the orbits contained at the infinity S 1 . The set formed by all separatrices of π(p(X)) is denoted by Σ X and it is closed (see Neumann [5]).
An open connected component of D \ Σ X is a canonical region of π(p(X)) or X. The separatrix configuration of π(p(X)) or X is formed by the union of one orbit chosen from each canonical region with Σ X , and it is denoted by Σ ′ X . We say that two separatrix configurations Σ ′ X and Σ ′ Y are topologically equivalent if there exists a homeomorphism in D preserving the infinity S 1 carrying trajectories of Σ ′ X into trajectories of Σ ′ Y , either preserving or reversing the sense of all trajectories.
The topologically equivalence between two Poincaré compactified vector fields has been characterized by Markus [4], Neumann [5] and Peixoto [6] who proved that two Poincaré compactifications π(p(X)) and π(p(Y )) having finitely many separatrices are topologically equivalent if and only if their correspoinding separatrix configurations Σ ′ X and Σ ′ Y are topologically equivalent.