Unbounded solutions of models for glycolysis

The Selkov oscillator, a simple description of glycolysis, is a system of two ordinary differential equations with mass action kinetics. In previous work the authors established several properties of the solutions of this system. In the present paper we extend this to prove that this system has solutions which diverge to infinity in an oscillatory manner at late times. This is done with the help of a Poincaré compactification of the system and a shooting argument. This system was originally derived from another system with Michaelis–Menten kinetics. A Poincaré compactification of the latter system is carried out and this is used to show that the Michaelis–Menten system, like that with mass action, has solutions which diverge to infinity in a monotone manner. It is also shown to admit subcritical Hopf bifurcations and thus unstable periodic solutions. We discuss to what extent the unbounded solutions cast doubt on the biological relevance of the Selkov oscillator and compare it with other models for the same biological system in the literature.


Introduction
When trying to understand a biological system with the help of mathematical modelling it often happens that there are several different models for the same biological situation in the literature. In view of this it is important to have criteria for deciding between models. One strategy for identifying criteria of this type is to look at relatively simple examples in great detail. In order to do this effectively it is necessary to have a sufficiently comprehensive understanding of the properties of solutions of the B Alan D. Rendall rendall@uni-mainz.de 1 Institut für Mathematik Johannes Gutenberg-Universität, Staudingerweg 9, 55099 Mainz, Germany models being studied. In this paper, with this strategy in mind, we look in detail at the dynamical properties of certain models for glycolysis.
Glycolysis is part of the process by which living organisms extract energy from sugar (Alberts et al. 2008). A suitable model system for studying this phenomenon experimentally is yeast extract or a suspension of yeast cells. The first indication that this system might have interesting dynamical properties was given by damped oscillations reported in Duysens and Amesz (1957). Later it was discovered that a constant continuous supply of sugar can lead to sustained oscillations (cf. Boiteux et al. 1975). Looking for the source of these oscillations revealed that they are produced by a small reaction network describing the action of the enzyme phosphofructokinase. A mathematical model for this network was set up and studied by Higgins (1964). It was found by Selkov that this model was not adequate for describing the oscillations and he introduced a modified one (Selkov 1968). The starting point for the model of Selkov is a reaction network with five chemical species. Assuming mass action kinetics leads to a system of five ordinary differential equations. Using quasi-steady state assumptions this can be reduced to a system of two equations with nonlinearities of Michaelis-Menten type. For brevity we call it 'the Michaelis-Menten system' in what follows. Setting one of the coefficients in this system to zero leads to a further simplification, giving a system of two equations with mass action kinetics, which we call the 'basic Selkov system' in what follows.
The aim of this paper is to obtain a better understanding of the dynamics of solutions of the three systems just described. A number of properties of solutions of the basic Selkov system were already established in Selkov (1968) but for many years no further rigorous results on this subject were obtained. Important progress was made in a paper of d'Onofrio (2010) and a number of additional properties of the solutions were established in a recent paper of the authors (Brechmann and Rendall 2018). In particular it was proved that for any values of the parameters there exist unbounded solutions of this system which are eventually monotone in the sense that for a solution of this type both concentrations are monotone after a certain time. In Selkov (1968) it is claimed that this system admits solutions which oscillate with an amplitude which grows without limit at late times. In what follows solutions of this type are referred to as 'solutions with unbounded oscillations'. The paper Selkov (1968) provides no justification for the claim other than a mention of numerical simulations, about which no details are given. Up to now there was no proof of the truth or falsity of this claim of Selkov (1968). One of the main results of the present paper is a proof of the existence of solutions of the basic Selkov system with unbounded oscillations. Our discovery of this proof was stimulated by the paper (Merkin et al. 1987), which belongs to the domain of theoretical chemistry. It deals with a system which turns out to be identical to the basic Selkov system when a parameter γ in the latter system takes the value two.
In Merkin et al. (1987) a claim of the existence of solutions with unbounded oscillations is also made. It is supported by an intricate heuristic argument using matched asymptotic expansions. It is not at all clear how this argument could be translated into a rigorous one but it provided us with some ideas which, when combined with the results of Brechmann and Rendall (2018), do give a proof of the existence of solutions with unbounded oscillations. When written in dimensionless form the system contains one parameter α. As claimed in Merkin et al. (1987), solutions with unbounded oscillations occur for precisely one value α 1 of α. When α is slightly less than α 1 there exists a stable periodic solution. As α approaches α 1 from below the amplitude of the periodic solution tends to infinity. One important element of this proof is to study the limit of the system for α → ∞ after a suitable rescaling. The existence of α 1 is then proved by a shooting argument. A monotonicity property, which was apparently not previously known, is used to obtain the uniqueness of α 1 .
The presence of unbounded solutions, whether monotone or oscillatory, might be seen as a feature which is unrealistic from the point of view of the biological applications. The monotone unbounded solutions of the basic Selkov system are not mentioned at all in Selkov (1968). That system is the limit of the Michaelis-Menten system when a parameter ν tends to zero. It is stated in Selkov (1968) that solutions with unbounded oscillations do not exist for ν > 0. On the other hand simulations reported in Keener and Sneyd (2009) suggest that the amplitude of periodic solutions of the Michaelis-Menten system diverges rapidly to infinity when a parameter is varied in a finite range. This indicates that, in contrast to the claim of Selkov, the existence of unbounded oscillations is a phenomenon which may persist for ν > 0. If this is true then the presence of these biologically problematic solutions of the basic Selkov system is not just an artefact of taking the limit ν → 0. The issue of the existence of solutions with unbounded oscillations in the case of the Michaelis-Menten system is not resolved in what follows but some partial results are obtained. In particular it is shown that for the Michaelis-Menten system with arbitrary parameters there are unbounded solutions which are eventually monotone and whose leading order asymptotics are identical to those found in the basic Selkov system. It is also shown that for certain combinations of the parameters (α, ν) with ν > 0 all positive solutions except the steady state have these late-time asymptotics. It turns out that there are parameter values for which there exist unstable periodic solutions of the Michaelis-Menten system. This is in contrast to the basic Selkov system where it was proved in Brechmann and Rendall (2018) that all periodic solutions are asymptotically stable.
The structure of the paper is as follows. The various systems considered in the paper are defined in Sect. 2. In Sect. 3, after some necessary results on the basic Selkov system proved in Brechmann and Rendall (2018) have been recalled, the existence of solutions with unbounded oscillations is proved. Similarities and differences between the properties of solutions of the basic Selkov system and the Michaelis-Menten system are discussed in the next three sections. Section 4 discusses the Hopf bifurcation exhibited by the Michaelis-Menten system. Its Poincaré compactification is computed in Sect. 5. Global properties of the Michaelis-Menten system are discussed in Sect. 6. The paper ends with a conclusion and outlook.

Survey of the systems considered
In Selkov (1968) a simple reaction network describing glycolysis is introduced. Assuming mass action kinetics for this network leads to a system of five ordinary differential equations, system (4) of Selkov (1968). In a slightly modified notation this system is In fact a factor γ was omitted in two places in Selkov (1968) and this error has been corrected here. All the parameters are positive and it is assumed that γ > 1, which encodes the biological property of cooperativity. Note that e 0 = e + x 1 + x 2 is a conserved quantity (total amount of enzyme) and this can be used to eliminate e from the first four evolution equations and discard the evolution equation for e. This reduces the system to four equations. Dimensionless variables can be introduced by defining This leads to the system Formally setting = 0 in the Eqs.
(3) and (4) gives u 2 = σ 1 u 1 and u 1 = σ γ 2 1+σ γ 2 +σ 1 σ γ 2 and substituting these relations into the evolution equations for σ 1 and σ 2 gives As has been discussed in Brechmann and Rendall (2018) geometric singular perturbation theory (GSPT) can be used to show that solutions of (1)-(4) converge to solutions of (5) and (6) in the limit → 0. In Selkov (1968) a further simplification of this system is introduced. Consider the rescaled quantities Expressing the Eqs. (5) and (6) in terms of these gives This system has a regular limit when ν tends to zero with α and β fixed. In the limit we get the basic Selkov system, system (II) of Selkov (1968), which is It is the system of central interest in Selkov (1968) and the dynamical properties of its solutions are studied in detail in Brechmann and Rendall (2018). Of course (9) and (10) can be thought of as the special case of (7) and (8) where ν = 0. Note that it follows from (7) and (8) that d dτ (αx + y) = α(1 − y). Since any solution with positive initial data remains positive as long as it exists this relation shows that it remains bounded on any finite time interval. Hence, when maximally extended, it exists globally in the future.

The basic Selkov system
The following proposition collects some of the properties of solutions of the basic Selkov system established in Brechmann and Rendall (2018).
Proposition 1 The basic Selkov system (9) and (10) has the following properties.
1. For each value of the parameter α the unique positive steady state P 0 has coordinates (1, 1). 2. For each α ∈ 0, 1 γ −1 the positive steady state P 0 is asymptotically stable and there exist no periodic solutions.
Using the standard theory of Hopf bifurcations it follows from statement 3. of the proposition that for any α slightly greater than α 0 there exists a stable periodic solution.
For the convenience of the reader we recall one of the main results of Brechmann and Rendall (2018) (Theorem 3 of that paper). The points P i in the statement are those occurring in Fig. 1 and whose definition is explained below.
Theorem 1 In the case α > α 0 of the system (9) and (10) exactly one of the following three situations occurs 1. The centre manifolds of P 1 and P 3 coincide so that there is a heteroclinic cycle at infinity. Any solution which starts below this centre manifold converges to P 4 as t → ∞ while any solution other than the steady state which starts above this manifold converges to the heteroclinic cycle at infinity as t → ∞. 2. The centre manifolds of P 1 and P 3 do not coincide. There exists a unique periodic solution. Any solution which starts below the centre manifold of P 3 converges to P 4 as t → ∞ while any solution other than the steady state which starts above this manifold converges to the periodic solution as t → ∞. 3. The centre manifolds of P 1 and P 3 do not coincide. Any solution other than the steady state which does not lie on the centre manifold of P 3 converges to P 4 as t → ∞.
A key question left open in Brechmann and Rendall (2018) is that of what happens to the periodic solution when α gets large. This question is answered in this section.
Theorem 2 There exists a number α 1 > α 0 such that the basic Selkov system (9) and (10) has the following properties.
1. For α = α 1 there exist solutions with the properties that lim inf t→∞ x(t) = lim inf t→∞ y(t) = 0 and lim sup t→∞ x(t) = lim sup t→∞ y(t) = ∞. 2. For α 0 < α < α 1 there exists a unique periodic solution and it is asymptotically stable. 3. For α > α 1 each solution other than the steady state P 0 is unbounded and has the properties described in statement 4. of Proposition 1. 4. As α tends to α 1 from below the diameter of the image of the periodic solution tends to infinity.
We adopt some of the notation of Brechmann and Rendall (2018). There the Poincaré compactification of the basic Selkov system is computed and one of the resulting points at infinity is blown up. After this has been done there are four steady states at infinity called P 1 , P 2 , P 3 and P 4 . Their positions can be seen in Fig. 1, which is a modification of Fig. 1 of Brechmann and Rendall (2018). Each of the points P 1 and P 3 has a one-dimensional centre manifold with the flow on the centre manifolds being away from P 1 and towards P 3 . As a starting point for the proof of Theorem 2 we  (9) and (10) Nullclines of the basic Selkov system (9) and (10) establish some further properties of the centre manifolds of the points P 1 and P 3 , both of which are unique. Let L be the segment of the line y = 1 where 0 < x ≤ 1. We use the notation for the components U i of the complement of the nullclines N 1 = N + 1 ∪ N − 1 and N 2 = N + 2 ∪ N − 2 which can be seen in Fig. 2, a modification of Fig. 2 of Brechmann and Rendall (2018). Here N 1 and N 2 are the zero sets of dx dτ and dy dτ , respectively.

Lemma 1
In the basic Selkov system the centre manifolds of P 1 and P 3 both contain a point of L in their closures.
Proof For a point on the centre manifold of P 1 sufficiently near to P 1 we haveẋ > 0. Hence the manifold initially lies in the region U 1 . As long as x < 1 it must remain in U 1 and both coordinates of a solution on the centre manifold are monotone functions of time. Hence a solution on the centre manifold of P 1 either reaches a point of L with x < 1 after a finite time or it tends to the positive steady state as t → ∞. Similarly a solution on the centre manifold of P 3 , when followed backwards in time, either reaches a point of L with x < 1 after a finite time or it tends to the positive steady state as t → −∞.
For a given value of α let ξ 1 (α) be the x-coordinate of the point where the centre manifold of P 1 meets L if such a point exists and otherwise let ξ 1 (α) = 1. Define ξ 2 (α) similarly in terms of the centre manifold of P 3 . Note that each centre manifold depends smoothly on the parameter α, in the sense that we can choose initial data for solutions on the centre manifold for different values of α in such a way that the solutions depend smoothly on α. This can be seen by considering the suspended system obtained by adjoining the equationα = 0 to the basic Selkov system and noting that it has a two-dimensional centre manifold at the points corresponding to P 1 and P 3 . This manifold is foliated by curves of constant α which are centre manifolds for the original system. Their smooth dependence on α follows from the smoothness of the two-dimensional centre manifold.

Lemma 2
The function ξ 1 −ξ 2 describing the separation of the points where the centre manifolds of P 1 and P 3 reach y = 1 is continuous.
Proof Consider a value of α c for which ξ 1 (α c ) < 1. The centre manifold for that value crosses L transversely and so, by the implicit function theorem, ξ 1 is a smooth function of α close to α c . This also shows that the set of values of α for which ξ 1 (α) < 1 is open. Consider now a value α * of α for which ξ 1 (α * ) = 1 and a sequence α n satisfying lim n→∞ α n = α * . It will be shown that lim n→∞ ξ 1 (α n ) = 1. Together with the information already obtained this implies that ξ 1 is continuous everywhere. The desired statement will be proved by contradiction. If ξ 1 (α n ) did not converge to one then by passing to a subsequence we could assume that lim n→∞ ξ 1 (α n ) = ξ s < 1. Consider now the sequence of solutions of the basic Selkov system with x n (0) = ξ 1 (α n ), y n (0) = 1 and α = α n and the solution with x s (0) = ξ s , y s (0) = 1 and α = α * . We are interested in these solutions for t ≤ 0. The sequence (x n , y n ) converges to (x s , y s ) uniformly on compact time intervals. We claim that (x s , y s ) lies on the centre manifold of P 1 for α = α * . If (x s , y s ) lies to the left of the centre manifold then it reaches negative values of x for finite negative values of t. Then for n sufficiently large the solutions (x n , y n ) would do the same, a contradiction. If (x s , y s ) lies to the right of the centre manifold then it must reach values of x greater than ξ s for finite negative values of t. Then for n sufficiently large the solutions (x n , y n ) would do the same, a contradiction. The conclusion is that the solution (x s , y s ) lies on the centre manifold and hence ξ 1 (α * ) < 1, in contradiction to the definition of α * . It has thus been proved that ξ 1 is continuous. A similar argument shows that ξ 2 is continuous. Hence ξ 1 − ξ 2 is continuous.
Proof Suppose that for a given value of α we have (ξ 1 − ξ 2 )(α) ≤ 0. The region of the Poincaré compactification bounded by the parts of the centre manifolds of P 1 and P 3 ending on L and the part of L between them and above the centre manifold of P 3 is invariant under evolution backwards in time. Consider the solution obtained by backward time evolution of a point in this region other than the steady state. By Poincaré-Bendixson theory its α-limit set must be a steady state or a periodic solution. If α ≤ α 0 this leads to a contradiction, because in that case no periodic solutions exist and the positive steady state is a sink. Thus we can conclude that the function ξ 1 − ξ 2 is positive for 0 < α ≤ α 0 .
Next we investigate the behaviour of solutions for α large. The following calculations were inspired by a transformation introduced in Merkin et al. (1987) in the case and we are interested in the limit μ → 0. A Poincaré compactification of this system was carried out in Merkin et al. (1987). After a suitable rescaling this leads to a system in the standard form of a fast-slow system in GSPT. [For background on GSPT we refer to Kuehn (2015).] Unfortunately in this system the important property of normal hyperbolicity breaks down at the point corresponding to P 3 . It turns out that this problem can be got around by using the transformations introduced in Brechmann and Rendall (2018) to treat the behaviour of solutions for x large. These can be summed and choosing a time coordinate s satisfying ds dτ = 1 γ x y γ −1 . This transforms the basic Selkov system into system (12) and (13) of Brechmann and Rendall (2018). Now introduce = α −1 andw = α(z−1). Then, denoting the derivative with respect to s by a prime, we get the system Note that, due to cancellations in the expressions in square brackets, this system depends smoothly on at = 0 and in fact the apparently singular term even vanishes as → 0. The critical manifold has the equation . The derivative of the right hand side of the Eq. (14) with respect tow, evaluated at = 0, is γ (γ − 1). Thus the critical manifold is normally hyperbolic repelling.
[For the terminology see Kuehn (2015).] This implies that when the system (13) and (14) is restricted to the slow manifold its limit as → 0 becomes regular. The evolution equation on the critical manifold is On the critical manifold there are two steady states, a source and a sink. They are connected by a heteroclinic orbit. For small and positive the critical manifold perturbs to a one-dimensional invariant manifold, which is the restriction of the slow manifold to that value of . The points (1, 0) and (0, 0) are steady states of (13) and (14) for all values of . We claim that for small they lie on the slow manifold and hence on the invariant manifold mentioned above. If we consider the extended system obtained by adjoining the equation = 0 to the system (13) and (14) then the slow manifold is a centre manifold of each point of the critical manifold with respect to the extended system. Since for a steady state of a dynamical system any steady state sufficiently close to it must lie on any centre manifold of the original point the claim follows. For (1, 0, ) and (0, 0, ) are close to (1, 0, 0) and (0, 0, 0), respectively when is small. The steady state atȳ = 1 is hyperbolic and so perturbs to a hyperbolic source. The steady state atȳ = 0 continues to exist. There are no other steady states between these two because we know that there is only one positive steady state for > 0 due to the uniqueness of P 0 . It follows that there is also a connection between the positive steady state P 0 of the Selkov system and the point P 3 on the boundary for α sufficiently large.
(The direction of the flow on the connecting orbit is determined by the fact that P 0 is a hyperbolic source.) In other words, when α is sufficiently large the centre manifold of P 3 converges to the positive steady state in the past. This means that ξ 2 (α) = 1. On the other hand, since the positive steady state is a source in this case the centre manifold of P 1 cannot converge to the positive steady state. We conclude that ξ 1 (α) < 1 and that ξ 1 − ξ 2 is negative. By the intermediate value theorem there exists some α 1 with ξ 1 (α 1 ) = ξ 2 (α 1 ). Note that in the end the Eqs. (11)-(12) were not needed in the proof but we judged it useful to include them so as to give an indication of how the argument was found.
It turns out that the value of α for which the centre manifolds of P 1 and P 3 meet is unique. This follows from a monotonicity property of the dependence of the centre manifolds on α.

Lemma 4
The function ξ 1 −ξ 2 describing the separation of the points where the centre manifolds of P 1 and P 3 reach y = 1 is strictly decreasing and has a unique zero.
Proof For this proof it is convenient to think of y as a function of x for a given solution. Suppose that a solution y(x) for a parameter α and a solutionŷ(x) for a parameter α < α satisfy y(x 1 ) =ŷ(x 1 ) for some x 1 . Thenŷ (x 1 ) < 0 and |ŷ (x 1 )| < |y (x 1 )|.
The leading order approximation to the centre manifold of P 3 is given byz = 1 + ν 1ȳ γ −1 + . . . where ν 1 = 1 αγ (γ −1) . This translates [in terms of variables Y =ȳ γ , Z =ȳ γ −1z used in Brechmann and Rendall (2018) . . . and x = y 1−γ − γ ν 1 + . . .. Putting these things together shows that when α is reduced the intersection of the centre manifold of P 3 with the line y = 1 moves to the left. To obtain information about the position of the centre manifold of P 1 in its dependence on α it is necessary to determine one more order in the expansion of the centre manifold than was done in Brechmann and Rendall (2018). The result is X = Z γ +1 −γ αZ 2γ +1 +. . ..
In the original variables this gives x = y −γ − γ αy −2γ + . . .. When α is reduced x becomes larger for fixed y. This also means that y becomes larger for fixed x and this propagates to larger values of x. Thus the intersection of the centre manifold of P 1 with the line y = 1 moves to the right. This implies that the function ξ 1 − ξ 2 is strictly decreasing and cannot have more than one zero.
Proof of Theorem 2 By Lemmas 3 and 4 there exists a unique α 1 > α 0 for which the centre manifolds of P 1 and P 3 coincide. With this information the first statement of Theorem 2 follows immediately from the first statement of Theorem 1. For α 0 < α < α 1 the positive steady state is unstable and there is no heteroclinic cycle at infinity. It follows from the Poincaré-Bendixson theorem that the ω-limit set of a solution which starts near the steady state but is not the steady state itself must be a periodic solution. In particular, a periodic solution exists and we are in the second case of Theorem 1. Thus the second statement of Theorem 2 holds. If α > α 1 then there is again no heteroclinic cycle at infinity. The α-limit set of the solution on the centre manifold of P 3 must then, by the Poincaré-Bendixson theorem, be either a periodic solution or the positive steady state. Moreover, if a periodic solution exists then only the first possibility can occur. Since, however, it follows from Brechmann and Rendall (2018) that any periodic solution which exists is stable the first possibility is ruled out. There can be no periodic solution and the third case of Theorem 1 of must be realised. This completes the proof of the third statement. Finally, the fourth statement will be proved by contradiction. Let β i be a sequence tending to α 1 from below. For a given i the system with parameter β i has a unique periodic solution and there is a unique point in its image of the form (1, z i ) with z i > 1. If this sequence did not tend to infinity then it would have a convergent subsequence. Thus after passing to a subsequence z i tends to a finite limit z * . The periodic solutions through the points (1, z i ) converge to a solution through the point (1, z * ), which is a periodic solution of the system with parameter value α 1 . This contradicts the fact that there are no such solutions.

The Michaelis-Menten system
In the system (7) and (8) the x-axis is an invariant manifold of the flow and the vector field is directed towards positive values of x on the y-axis. For each fixed choice of the parameters with ν < 1 there is a unique positive steady state at 1+βν 1−ν , 1 . For ν ≥ 1 there is no positive steady state. Linearizing the system about the steady state leads to the Jacobian The determinant of J is α (1−ν) 2 1+βν which is always positive. Thus the stability of the steady state is determined by the trace of J , which is If γ ≤ 1+βν 1−ν then the trace of J is negative for all values of α and the steady state is always stable. 1+βν) . Then for α < α 0 the trace of J is negative and the steady state is asymptotically stable while for α > α 0 the trace of J is positive and the steady state is a source. For α = α 0 there is a pair of imaginary eigenvalues. If we consider the real part of the eigenvalues as a function of α then it passes through zero when α = α 0 and its derivative with respect to α at that point is non-zero. Thus a Hopf bifurcation occurs.
In the limiting case ν = 0 it was shown in Brechmann and Rendall (2018) that the Hopf bifurcation is supercritical so that there exists a stable periodic solution for any α slightly greater than α 0 . The computation of the Lyapunov number required to obtain this conclusion becomes considerably more complicated for ν > 0. Rather than trying to do this in general we will confine ourselves to obtaining some information for restricted sets of parameters. The Lyapunov number of the Hopf bifurcation is a function of the parameters α, β, γ and ν and we are interested in its sign. A general formula for this quantity is given in Sect. 4.4 of Perko (2001). It is of the form −3π 2b 3/2 f , where the first factor is positive in the present case and f is a function of (α, β, γ, ν) which is negative when ν = 0. This shows that in that case the Hopf bifurcation is supercritical. For ν small and positive f is still negative and the bifurcation supercritical. It will now be proved that there also exist parameters for which f is positive, so that there exists a subcritical Hopf bifurcation. In that case there exist unstable periodic solutions for α slightly less than α 0 . Note for comparison that it was shown in Brechmann and Rendall (2018) that for ν = 0 unstable periodic solutions never exist. It suffices to treat the case β = 0 since an example with β small and positive follows by continuity. Since we are only looking for some example we can also restrict to the case γ = 2.
With a suitable normalization the function f is of the following form.
The expression in square brackets tends to a positive value as ν → 1 2 . Thus the leading term in the expression for the Lyapunov number is positive for ν close to its limiting value. This proves that there are parameters for which the Hopf bifurcation is subcritical and thus that there exist unstable periodic solutions.

The Poincaré compactification
In Brechmann and Rendall (2018) it was investigated using the Poincaré compactification in which ways solutions of (9) and (10) can tend to infinity for large times. Here we want to carry out corresponding calculations for (7) and (8). A useful preliminary step is to introduce a new time coordinate T satisfying dτ dT = 1 + ν y γ (β + x). Then we get the system This makes the right hand side into a polynomial while leaving the phase portrait unchanged.
The phase portrait is more complicated than that in the case of mass action kinetics. A schematic picture of it is given in Fig. 3 and its properties are summarized in the following lemma which is the analogue of Lemma 2 in Brechmann and Rendall (2018).

Lemma 5 Suppose that ν < 1. There is a smooth mapping of the closure of the positive quadrant into itself mapping the axes into themselves with the following properties. The restriction of φ to the open quadrant is a diffeomorphism onto its image. This image is a region whose closure is a compact set bounded by intervals [0,
x 0 ] and [0, y 0 ] on the x-and y-axes and four smooth curves γ i , 1 ≤ i ≤ 4. The curve γ 1 joins the point P 1 = (0, y 0 ) with a point P 3 in the positive quadrant. γ 2 joins the point P 3 with the point P 4 . For 3 ≤ i ≤ 4 the curve γ i joins the point P 2i−2 with the point P 2i and P 8 = (x 0 , 0). The image of the dynamical system can be rescaled so as to extend smoothly to the closure of the image of φ in such a way that P 3 and P 2i , 2 ≤ i ≤ 4, are steady states and the γ i and the image of the x-axis under φ are invariant manifolds.
The statements of Lemma 5 are proved in the remainder of this section. To analyse the case where x becomes large [Case 1 in the terminology of Brechmann and Rendall (2018) Both axes are invariant under the flow and the flow there is towards the origin. The linearization of the system about the origin is identically zero. Thus we do a quasihomogeneous directional blow-up. An appropriate transformation can be obtained by using a Newton polygon as in Dumortier et al. (2006). The coefficients areα = γ and β = γ −1. (These are the same values as occurred in the blow-up of the corresponding point for the model (9) and (10). The notation has been changed compared to that of the sources quoted by the addition of a tilde to avoid confusion with other uses of the same letters elsewhere in the present paper.) Thus we use variablesȳ andz satisfying (Y , Z ) = (ȳ γ ,ȳ γ −1z ). In addition we introduce a new time coordinate s satisfying ds dt = γ −1ȳγ 2 −1 . The system becomes Both axes are invariant under the flow. There is a steady state at the origin and one at the point (0, 1), which corresponds to P 7 . The linearization at the origin is identically zero.
Next the centre manifold of P 7 will be studied. Introducing w =z − 1 moves the steady state to origin. The centre subspace is given w = ρȳ with ρ = 1−αν 2α for γ = 2 and ρ = − ν γ for γ > 2.
Proof In the case γ = 2 we havē Substituting the relationz = 1 + ρȳ + O(ȳ 2 ) which holds on the centre manifold into this relation gives and this completes the proof for γ = 2. For γ > 2 we use the relation
We see that the flow on the centre manifold is towards P 7 and since the non-zero eigenvalue of the linearization at that point is positive P 7 is a topological saddle. In fact the flow on the boundary is everywhere away from P 7 . We next blow up the origin in the coordinates (ȳ,z). This time the procedure described in Dumortier et al. (2006) leads to the choice of coefficientsα =β = 1. Blow-ups in the two coordinate directions are required. The only terms in the resulting equations which will be written explicitly are those which have a direct influence on the analysis which follows. In the first transformed system, withȳ =ỹ 1 andz =ỹ 1z1 , the equations arẽ A change of time coordinate eliminates the common factorỹ 1 . On the boundary there is a steady state at the point (0, ν), which corresponds to P 5 . The origin of this coordinate system corresponds to P 4 . The terms which have been retained suffice to determine the steady states on the boundary and the linearization of the system at those points. The point P 5 also appears in the second transformed system but since it can be analysed in the chart corresponding to the first transformed system the second transformed system, withȳ =ỹ 2z2 andz =z 2 , is only needed to analyse the steady state P 6 at the origin of that coordinate system. For this purpose the only terms which need to be retained areỹ The common factorz 2 can be eliminated by a change of time coordinate. The origin of this coordinate system corresponds to P 6 . We see that in both cases, after a suitable change of time coordinate, the origin is a hyperbolic saddle.
Next the centre manifold of P 5 will be studied in the case γ = 2. We do not expect that the case γ > 2 is essentially different but since the algebra becomes significantly more complicated only the case γ = 2 has been worked out. The centre subspace is parallel to theỹ 1 -axis. We havez 1 = O(ỹ 3 1 ) on the centre manifold and this implies that ifz 1 = ν + w then 2ανw = [ν 2 (1 − ν) + 2αν 4 + 2αβν 3 ]ỹ 2 1 + . . .
It follows that provided ν < 1 the flow on the centre manifold of P 5 is away from P 5 . For the rest of the discussion we return to the case of general γ .
In the case where x gets large it remains to do one further quasihomogeneous directional blow-up of the origin in the (Y , Z ) coordinate system. In this case (Y , Z ) = (ȳz γ ,z γ −1 ). The time coordinate is transformed using the relation ds dt = 1 γ −1z γ 2 −1 . The resulting system is There is a steady state at the point (1, 0) but since it is just another representation of P 7 it does not need to be analysed further. The origin of this coordinate system corresponds to P 8 . Thez-axis is a centre manifold for P 8 and the flow there is towards P 8 . Since the non-zero eigenvalue of the linearization at P 8 is negative it can be concluded that P 8 is a sink.
Having completed the analysis of the case where x gets large we now turn to the case where where y gets large [Case 2 in the terminology of Brechmann and Rendall (2018)], with new variables X = x y and Z = 1 y . The result is The common factor 1 Z γ +1 can be removed by a suitable change of time coordinate satisfying dt dT = Z −γ −1 . The linearization of the resulting system about the origin is identically zero so that it is again necessary to do a blow-up. In this case a calculation using a Newton polygon gives the exponentsα = 1 andβ = 1. The transformation in the X direction uses the relation (X , Z ) = (x 1 ,x 1z1 ). The resulting system is The origin of this coordinate system corresponds to P 3 . By a change of time coordinate we can remove the factorx 1 . The linearization of the system which results at the origin has one positive eigenvalue and thez 1 -axis is invariant and defines a centre manifold at that point. It can be concluded that P 3 is a source. If ν < 1 there is a steady state at the point 0, 1−ν βν which corresponds to the point P 2 . That point is a hyperbolic saddle whose stable manifold is thez 1 -axis.
The transformation in the Z direction uses the relation (X , Z ) = (x 2z2 ,z 2 ). The resulting system is The origin of this coordinate system is P 1 . By a change of time coordinate we can remove the factorz 2 . In the system which results there is inflow on thez 2 -axis while thex 2 -axis is invariant and corresponds to thez 1 -axis in the previous system. Note that the point P 1 is not a steady state. The facts which have now been collected imply strong restrictions on the possible ω-limit sets of solutions. The only points of the boundary which they can contain are those on the part connecting P 5 and P 7 . Poincaré-Bendixson theory implies that the ω-limit set of a positive solution must be either a point (which can only be the positive steady state, P 7 or P 8 ), a periodic solution or a heteroclinic cycle joining P 5 and P 7 .
The last of these can only occur if the centre manifolds of P 5 and P 7 coincide. Note that any periodic solution or heteroclinic cycle must contain the positive steady state in its interior.
We have the following analogue of Theorem 1 of Brechmann and Rendall (2018).
Theorem 3 There exists a positive number > 0 such that any solution of the Michaelis-Menten system (7) and (8) with initial data x(0) = x 0 and y(0) = y 0 which satisfies x 0 > −1 and x 0 y γ 0 < has the late-time asymptotics for a constant y 1 . There exists a solution, unique up to time translation, which has the asymptotic behaviour (1)).
Proof The proof of this theorem is very similar to that of Theorem 1 of Brechmann and Rendall (2018), whose basic structure we now recall. Any solution which starts close enough to the point P 8 converges to that point as t → ∞. Using this information in the system (16) and (17) allows these equations forȳ(s) andz(s) to be integrated to leading order in the limit s → ∞. The resulting asymptotic expressions can then be transformed back to the original variables x(τ ) and y(τ ). The only extra element is that, while in the original proof only one change of time coordinate was required, in the present case we must first transform from s to T and then from T to τ . Since for this type of solution the time coordinates τ and T are asymptotically equal this extra element does not change the final answer. In particular, the parameter ν does not contribute to the leading order asymptotics. This gives the first statement of the theorem. The solution mentioned in the second statement of the theorem is a solution on the centre manifold of P 7 . Integrate the equation forȳ in Lemma 6 in leading order in the limit s → ∞ and substitute the result into the equation forz. This provides asymptotic expressions forȳ(s) andz(s). Transforming these back to the variables x(τ ) and y(τ ) gives the second part of the theorem.

The global phase portrait
To understand the global phase portrait it is helpful to understand the geometry of the nullclines N 1 and N 2 given byẋ = 0 andẏ = 0, respectively. We restrict consideration to the case ν < 1 where N 1 and N 2 intersect in a single point. The equation for N 1 can be expressed in the equivalent forms

Fig. 4
Nullclines of the Michaelis-Menten system (7) and (8) Thus on N 1 the coordinate y can be expressed as a smooth function of x with a smooth inverse. Note, however, that while the second function is defined for all positive y the first is only defined for x > βν 1−ν . The equation for N 2 can be expressed in the form Thus on N 2 the coordinate x can be expressed as a (locally defined) smooth function of y. Since x can be written as a function of y in both cases and there is only one point of intersection it is clear that the complement of N 1 ∪ N 2 is a union of the four connected components defined by the signs ofẋ andẏ. A schematic picture of the null clines is given in Fig. 4. As in Fig. 2 we denote the regions with the sign combinations (+, −), (+, +), (−, +) and (−, −) by U 1 , U 2 , U 3 and U 4 , respectively. Whereẏ = 0 anḋ x = 0 we can use the fact that N 2 is a graph over the y-axis to conclude that a solution can only pass from U 3 to U 4 and U 1 to U 2 and not the other way round. Similarly the fact that N 1 is a graph over the y-axis implies that a solution can only pass from U 4 to U 1 and U 2 to U 3 and not the other way round. Thus the possible passages between the regions U i are just as in the case with mass action kinetics. Let L be the part of the horizontal line segment joining the positive steady state to the y-axis with the endpoint on the axis excluded. The part of L excluding the steady state is contained in U 1 .

Lemma 7
In the Michaelis-Menten system each of the centre manifolds of P 5 and P 7 contains a point of L in its closure.
Proof A point on the centre manifold of P 5 which is sufficiently close to P 5 lies in the region U 3 . If we follow a solution which lies on this manifold forwards in time then it must either tend to the positive steady state as t → ∞ or it must enter U 4 after a finite time and in the latter case it must enter U 1 . Once it has done so it must either tend to the positive steady state as t → ∞ or it must meet L after a finite time. Similarly a solution on the centre manifold of P 7 which starts close to P 7 must, when followed backwards in time, either converge to the positive steady state as t → −∞ or meet L after a finite time.
Denote the x-coordinates of the points of L in the closure of the centre manifolds of P 5 and P 7 for given value of α and ν by ξ 1 (α, ν) and ξ 2 (α, ν).

Lemma 8
The function ξ 1 −ξ 2 describing the separation of the points where the centre manifolds of P 5 and P 7 reach y = 1 is continuous.
Proof The proof is similar to that of Lemma 2. The essential facts which must be used are that any solution which approaches P 5 close enough in the past time direction and does not lie on the centre manifold of P 5 cannot remain close to P 5 forever and the corresponding statement with 'past' replaced by 'future' and P 5 by P 7 .
Proof For α sufficiently large the positive steady state is a source and thus ξ 1 (α, ν) < 1. Thus in order to prove the first part of the lemma it suffices to show that for some (α, ν) we have ξ 2 (α, ν) = 1. To prove this we proceed as in the case of mass action kinetics. First the system is transformed to the coordinates (ȳ,z) and then the quantities and w are introduced. The right hand side of each equation in the Michaelis-Menten case is the sum of the right hand side of the corresponding equation in the mass action case and an expression which can be written as −1 ν times a function which is regular in the limit → 0. Fixing ν and letting tend to zero would cause this term to explode. Instead we let and ν tend to zero in such a way that ν = 2 . Then the second summand behaves in a smooth manner as → 0 and in fact tends to zero. Thus proceeding in the same way as in the proof of Lemma 3 gives the first conclusion. For the second part we can again proceed as in the proof of Lemma 3. The difference is that while in the mass action case we knew that there was no unstable periodic solution in the Michaelis-Menten case we have to assume it.
Note that for the choice of parameters in the first part of Lemma 9 there is a heteroclinic orbit joining the positive steady state to the point P 7 . It follows that for these values of the parameters no periodic solutions exist. This is because a periodic solution would have to contain the positive steady state in its interior and therefore would have to cross the heteroclinic orbit.
There is no straightforward generalization of the monotonicity result of Lemma 4 to the Michaelis-Menten case. The proof of monotonicity fails for the centre manifold of P 5 since it may pass through the region U 4 . For this reason even in a case where the existence of a zero of ξ 1 −ξ 2 can be proved we do not get its uniqueness, Moreover, we do not get the analogue of the stability statement in the mass action case. It is possible to do a calculation analogous to that done to determine the stability of the heteroclinic cycle in Brechmann and Rendall (2018). Unfortunately in the estimate for the return map the power γ is replaced by the power one and this gives no information about stability. The following theorem sums up the results obtained. (7) and (8)  It should be noted that it has not been shown here whether the case described in point 6. ever occurs.

Conclusions and outlook
It has been shown that the basic Selkov system admits solutions with unbounded oscillations and that the diameter of the image of a periodic solution can tend to infinity as α approaches a finite limit, thus completing the results of Brechmann and Rendall (2018) on that system and rigorously confirming a claim made in Selkov (1968). Note that some statements related to this issue have been made in Erneux (2018) but that reference does not contain rigorous proofs of those statements. One remaining question is that of the rate with which the diameter of the image of the periodic solution tends to infinity as the critical parameter value α 1 is approached. A suggestion for this has been made in Merkin et al. (1987) for the case γ = 2 but there is neither a rigorous proof that this suggestion is correct nor a generalization of the statement to higher values of γ . It was also investigated which properties of the basic Selkov system persist in the Michaelis-Menten system from which Selkov derived his basic model. Partial results were obtained and it was shown in particular that the Michaelis-Menten system has unbounded solutions which are eventually monotone for all parameter values. The question of whether the five-dimensional system from which the Michaelis-Menten system itself was derived has unbounded solutions remains open. It was shown that for suitable parameter values unstable periodic solutions of the Michaelis-Menten system exist. It was left open whether there exist unbounded oscillatory solutions or periodic solutions whose images have arbitrarily large diameter for bounded ranges of the parameters.
The unbounded solutions cast doubt on the suitability of the Selkov model for describing glycolytic oscillations. An alternative model often preferred to the Selkov model is that of Goldbeter and Lefever (1972). There the amplitude of the periodic solutions created in a Hopf bifurcation increases to a maximum before decreasing again to zero at a point where the periodic solutions vanish again in a second Hopf bifurcation. It has been proved by d'Onofrio (2011), on the basis of an analysis of a more general class of systems in Othmer and Aldridge (1978), that all solutions of the Goldbeter-Lefever model are bounded. Further aspects of the dynamics of solutions of that model have been studied in d'Onofrio (2011) and a sophisticated analysis of some of its properties has been carried out in Kosiuk and Szmolyan (2011).
The questions of the origin of the unbounded solutions and how they could be eliminated by modifying the system have been discussed in Merkin et al. (1986). The origin of the unbounded growth can be seen in the constant source term in the equation for x. This corresponds to an unlimited supply of the substrate. In Merkin et al. (1986) this is called the pooled chemical approximation. If this is replaced by a mechanism where the substrate is formed from a precursor which itself is limited in quantity then the oscillations only grow within a finite time period before decaying again. An alternative modification is to introduce an additional uncatalysed conversion of the substrate into the product. The resulting system is called the (cubic) autocatalator. According to the analysis of Merkin et al. (1986) this leads to a situation similar to that described above for the Goldbeter-Lefever model and the unbounded oscillations are absent. Some aspects of this type of model have been analysed rigorously in Gucwa and Szmolyan (2009).
The Selkov system and other related ones are model cases for understanding oscillations in biological and chemical systems. The study of equations of this type raises a number of issues. In what ways can heuristic and numerical results be made into rigorous theorems? How can we understand the relations between the choices made in modelling and the relevance of the resulting models for the applications? The present paper is intended as a contribution to the clarification of these issues.
Acknowledgements One of the authors (ADR) is grateful to James Sneyd and János Tóth for helpful discussions.
Funding Open Access funding enabled and organized by Projekt DEAL.
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/.