Dynamics of the Selkov oscillator
aa r X i v : . [ q - b i o . S C ] M a r Dynamics of the Selkov oscillator
Pia Brechmann and Alan D. RendallInstitut f¨ur MathematikJohannes Gutenberg-Universit¨atStaudingerweg 9D-55099 MainzGermany
Abstract
A classical example of a mathematical model for oscillations in a bi-ological system is the Selkov oscillator, which is a simple description ofglycolysis. It is a system of two ordinary differential equations which,when expressed in dimensionless variables, depends on two parameters.Surprisingly it appears that no complete rigorous analysis of the dynam-ics of this model has ever been given. In this paper several properties ofthe dynamics of solutions of the model are established. With a view tostudying unbounded solutions a thorough analysis of the Poincar´e com-pactification of the system is given. It is proved that for any values ofthe parameters there are solutions which tend to infinity at late timesand are eventually monotone. It is shown that when the unique steadystate is stable any bounded solution converges to the steady state at latetimes. When the steady state is unstable it is shown that for given valuesof the parameters either there is a unique periodic solution to which allbounded solutions other than the steady state converge at late times orthere is no periodic solution and all solutions other than the steady stateare unbounded. In the latter case each unbounded solution which tends toinfinity is eventually monotone and each unbounded solution which doesnot tend to infinity has the property that each variable takes on arbitrarilylarge and small values at arbitrarily late times.
Glycolysis is a part of the process by which living organisms extract energyfrom sugar. It has been observed that under suitable circumstances the rate atwhich products of glycolysis accumulate shows oscillations in time although theinput rate of sugar to the system is constant. Soon after these experimentalobservations had been made Higgins [8] used mathematical modelling to obtaina deeper understanding of the process. Later it was observed by Selkov [16]that the model of Higgins did not show sustained oscillations for biologicallycorrect parameter values and he introduced an alternative mathematical model,1hich we refer to as the Selkov model. This system of two ordinary differentialequations is the subject of what follows. (We comment further on the work ofHiggins in the last section of the paper.) It was shown in [16] that the Selkovmodel exhibits a Hopf bifurcation and thus there exist parameter values forwhich it has a periodic solution. Results on the uniqueness and stability ofperiodic solutions of the Selkov model were obtained by d’Onofrio [3].It turns out that solutions of the Selkov model can be unbounded and thus itbecomes relevant to investigate the behaviour of the system near infinity. Thiscan be done by means of a Poincar´e compactification. This leads to two newsteady states at infinity which are not hyperbolic. It is relatively easy to deter-mine the behaviour of solutions close to one of these steady states, which has aone-dimensional centre manifold, and this was done in [16]. The other steadystate is more complicated since the linearization about that point vanishes iden-tically and it was not treated fully in [16]. To do so it is necessary to definesome blow-ups of the singularity and this is done in what follows. It is sug-gested in [16], on the basis of numerical calculations, that for certain parametervalues there are unbounded solutions which approach infinity in an oscillatoryfashion. Analytically these might be interpreted as solutions which approach aheteroclinic cycle in (a suitable blow-up of) the Poincar´e compactification.Section 2 presents the basic facts about the model. In particular it is shownthat the model has a unique steady state for fixed values of the parameters andthat there is a Hopf bifurcation. The first Lyapunov coefficient is calculatedand shown to be negative. Thus the bifurcation is generic and supercriticaland there exist stable periodic solutions for parameter values close to thoseat the bifurcation. The subject of Section 3 is the Poincar´e compactificationof the system. In order to obtain a compactification where the dimension ofthe centre manifold of each steady state is at most one suitable blow-ups arecarried out. This allows the qualitative nature of the dynamics near the steadystates on the boundary of the compactification to be determined. In additionto this local information it is shown that if there is a heteroclinic cycle in thecompactification it is asymptotically stable.Section 4 goes on to study the global properties of the phase portrait. It isshown that for any values of the parameters there are solutions which are un-bounded at late times and eventually monotone. The leading order asymptoticsof these solutions is determined. The results of d’Onofrio on the uniquenessand stability of periodic solutions are reviewed and completed. It is proved thatwhen the steady state is stable there exist no periodic solutions and all boundedsolutions converge to the steady state at late times. It is also proved that whenthe steady state is unstable one of two mutually exclusive possibilities occurs.Either exactly one periodic solution exists and all bounded solutions convergeto it at late times or no periodic solution exists and all solutions ( x ( t ) , y ( t ))other than the steady state are unbounded. In the latter case for a given un-bounded solution either lim t →∞ x ( t ) = ∞ and lim t →∞ y ( t ) = 0 and x ( t ) and y ( t ) are eventually monotone or lim sup t →∞ x ( t ) = lim sup t →∞ y ( t ) = ∞ andlim inf t →∞ x ( t ) = lim inf t →∞ y ( t ) = 0. The final section discusses what remainsto be discovered about the Selkov model and how this model relates to other2odels of glycolysis in the literature, in particular to a higher-dimensional modelconsidered in [16].This paper is based in part on the MSc thesis of the first author [1]. The system considered in what follows is the system (II) of [16]. It is dxdτ = 1 − xy γ , (1) dydτ = αy ( xy γ − − . (2)It is written in dimensionless variables. The quantities x and y represent dimen-sionless concentrations of ATP and ADP, respectively, while τ is a dimensionlesstime variable. The aim is to study the future evolution of solutions of theseequations for which x and y are positive. The parameter α is a positive realnumber. The parameter γ is a number greater than one and to avoid technicalcomplications with differentiability properties it will be assumed here to be aninteger.The x -axis is an invariant manifold of the flow of the system (1)-(2) and thevector field is directed towards positive values of x on the y -axis. For each fixed α there is a unique positive steady state and it satisfies x = y = 1. Linearizingthe system about this point leads to the Jacobi matrix J = (cid:20) − − γα α ( γ − (cid:21) . (3)The determinant of J is α which is always positive. Thus the stability of thesteady state is determined by the trace of J , which is α ( γ − −
1. Let α = γ − .Then for α < α the trace of J is negative and the steady state is asymptoticallystable while for α > α the trace of J is positive and the steady state is a source.The steady state is hyperbolic for all α = α . For α = α there is a pair ofnon-zero imaginary eigenvalues. If we consider the real part of the eigenvaluesas a function of α then it passes through zero when α = α and its derivativewith respect to α at that point is non-zero. Thus a Hopf bifurcation occurs.More information can be obtained by computing the first Lyapunov number σ of the bifurcation and it turns out that it is feasible to do this by hand using theformula given in [15]. The result is σ = − π ( γ − ( γ ( γ −
1) + 1) <
0. Thus theHopf bifurcation is non-degenerate and supercritical. This implies that for α slightly larger than α there exists a unique periodic solution close to the steadystate and that this periodic solution is stable and hyperbolic. For α = α thesteady state is not hyperbolic but it is topologically equivalent to a hyperbolicsink, since this is a property satisfied by all generic Hopf bifurcations.3 The Poincar´e compactification
The aim of this section is to investigate the ways in which the solutions of theSelkov system can tend to infinity at late times. The first step in doing this isto introduce two coordinate transformations, which we refer to as Case 1 andCase 2, respectively. In Case 1 we define Y = yx and Z = x . These variablesare appropriate for investigating the case where x becomes large. The originalvariables can be recovered using the relations x = Z and y = YZ . The result ofthe transformation is the system dYdτ = 1 Z γ ( αY γ + Y γ +1 − αY Z γ − Y Z γ +1 ) , (4) dZdτ = 1 Z γ − ( Y γ − Z γ +1 ) . (5)We introduce a new time variable t by dtdτ = Z − γ and obtain the system dYdt = αY γ + Y γ +1 − αY Z γ − Y Z γ +1 , (6) dZdt = Y γ Z − Z γ +2 . (7)In Case 2 we define X = xy and Z = y . These variables are appropriatefor investigating the case where y becomes large. The original variables can berecovered using the relations x = XZ and y = Z . The result of the transformationis the system dXdτ = 1 Z γ ( Z γ +1 − X − αX + αXZ γ ) , (8) dZdτ = 1 Z γ − ( − αX + αZ γ ) . (9)Introducing a new time variable as in Case 1 gives dXdt = Z γ +1 − X − αX + αXZ γ , (10) dZdt = − αXZ + αZ γ +1 . (11)Case 2 is easier to analyse than Case 1 and so we begin with that. The systemextends smoothly to the boundary where Z = 0. This boundary is an invariantmanifold for the extended flow and there X is decreasing except when X = 0.The origin in the new coordinates is a steady state, which we denote in whatfollows by P . The linearization of the right hand side of the equations at P hasrank one and so there is a one-dimensional centre manifold there. The centresubspace is given by X = 0 and the centre manifold is of the form X = h ( Z )where h vanishes faster than linearly at the origin. The invariance of the centremanifold under the flow implies that ˙ X = h ′ ( Z ) ˙ Z . Putting this informationinto the evolution equations gives Z γ +1 − h ( Z ) − α ( h ( Z )) + α ( h ( Z )) Z γ = h ′ ( Z )[ − α ( h ( Z )) Z + αZ γ +1 ] .
4f we expand this about Z = 0 we see that each of the terms other than Z γ +1 and h ( Z ) is negligible compared to one of these. It can be concluded that h ( Z ) = Z γ +1 + o ( Z γ +1 ). Substituting this relation into the evolution equationfor Z shows that on the centre manifold dZdt = αZ γ +1 + o ( Z γ +1 ). Hence the flowon the centre manifold is away from the steady state. For X = 0 the flow is intothe positive region while between the line Z = 0 and the centre manifold theflow is topologically equivalent to part of a hyperbolic saddle, as follows fromthe reduction theorem [14]. In particular, no solution starting in the positiveregion can tend to the point P at late times.We now turn to Case 1. The system extends smoothly to the boundarieswhere Y = 0 and Z = 0 and these are invariant. The origin is a steady state.Otherwise the variable Z is decreasing when Y = 0 and the variable Y isincreasing when Z = 0. The linearization at the steady state is identically zeroand gives no information on the dynamics in a neighbourhood of that point.To go further a blow-up procedure is applied to the steady state. One way ofdoing this is to pass to polar coordinates. It turns out that a new steady state isproduced where the linearization is zero. Introducing polar coordinates a secondtime then leads to a successful analysis of the dynamics in the case γ = 2, as wasshown in [1]. It seems, however, that for γ > γ in a unified way. It isnecessary to do one such blow-up for each of the coordinates. They depend ontwo integers α and β which are determined with the help of a Newton diagram[5]. In the present case α = γ and β = γ − Y, Z ) (¯ y, ¯ z ) with ( Y, Z ) =(¯ y γ , ¯ y γ − ¯ z ). We introduce a new time coordinate by dsdt = γ ¯ y γ − γ . Then thesystem becomes d ¯ yds = α ¯ y + ¯ y γ +1 − α ¯ y ¯ z γ − ¯ y γ ¯ z γ +1 , (12) d ¯ zds = − α ( γ − z − ( γ − y γ ¯ z + α ( γ − z γ +1 − ¯ y γ − ¯ z γ +2 + γ ¯ y γ ¯ z. (13)Evidently the origin of coordinates is a hyperbolic saddle which we denote by P . There is one other steady state P on the boundary where (¯ y, ¯ z ) = (0 , P there is a one-dimensional centre manifold, which will now be investigated.Introduce a new variable w by the relation ¯ z = 1 + w . We get the system¯ y ′ = − αγ ¯ yw − ¯ y γ (1 + w ) γ +1 − α ¯ y [(1 + w ) γ − − γw ] + ¯ y γ +1 ,w ′ = αγ ( γ − w − ¯ y γ − − ( γ − y γ (1 + w ) (14)+ α ( γ − w ) γ +1 − − ( γ + 1) w ] − ¯ y γ − [(1 + w ) γ +2 − γ ¯ y γ (1 + w ) . Important properties of the centre manifold are collected in the following lemma.
Lemma 1
On the centre manifold of P the relations w = ν ¯ y γ − + o (¯ y γ − )and ¯ y ′ = − ν ¯ y γ + o (¯ y γ ) hold with ν = αγ ( γ − and ν = γγ − .5 roof Consider first the case γ ≥
3. To start with it will be proved that w = O (¯ y γ − ) , w ′ = O (¯ y γ ) . (15)For γ ≥ w = 0. Thus on the centre manifold w = h (¯ y ) = O (¯ y ) with h ′ (¯ y ) = o (1) and hence ¯ y ′ = O (¯ y ). Differentiatingthe equation of the centre manifold gives w ′ = h ′ (¯ y ) y ′ and it follows that w ′ = O (¯ y ). This completes the proof of (15) if γ = 3. If γ ≥ w and obtain thestatement that w = O (¯ y ). This implies that ¯ y ′ = O (¯ y ) and that w ′ = O (¯ y ).Thus the proof of (15) is complete if γ = 4. For a general γ we can repeatthe argument just given until (15) is obtained. Substituting the informationobtained so far into the evolution equation for w shows that w ′ = ( γ − (cid:26) αγw − γ − y γ − + O (¯ y γ ) (cid:27) . (16)Combining this with (15) gives the first assertion of Lemma 1. Substituting theinformation already obtained back into the evolution equation for ¯ y gives¯ y ′ = − αγ ¯ yw − ¯ y γ + O (¯ y γ +1 ) (17)which implies the second statement.It remains to treat the case γ = 2. In that case the centre subspace is ofthe form ¯ y = µw with µ = 2 α . It is convenient to introduce a new coordinate v = ¯ y − µw , so that the centre subspace becomes v = 0. The centre manifoldis of the form v = h (¯ y ) = O (¯ y ). Thus w = µ ¯ y + O (¯ y ). This already givesthe first assertion of Lemma 1 in the case γ = 2 since in that case ν = µ .Substituting the equation of the centre manifold into the evolution equation for¯ y gives ¯ y ′ = − α ¯ yw − ¯ y + O (¯ y ) (18)which implies the second assertion in the case γ = 2. (cid:4) The second blow-up uses the transformation (
Y, Z ) (¯ y, ¯ z ) with ( Y, Z ) =(¯ y ¯ z γ , ¯ z γ − ). We introduce a new time coordinate by dsdt = γ − ¯ z γ − γ . Then thesystem becomes d ¯ yds = − ¯ y γ +1 ¯ z γ + ¯ y ¯ z γ − + α ( γ − y γ − α ( γ − y, (19) d ¯ zds = ¯ y γ ¯ z γ +1 − ¯ z γ . (20)On the boundary ¯ z = 0 there is a steady state with coordinates (1 , P and so it does notneed to be analysed further. We denote the steady state at the origin of thesecoordinates by P . At P the ¯ z -axis is a centre manifold. Any other centremanifold at that point is of the form ¯ y = h (¯ z ) where the function h vanishes toall orders for ¯ z →
0. 6 P P P γ γ Figure 1: Poincar´e compactification.This analysis can be summarized as follows (cf. Fig. 1).
Lemma 2
There is a smooth mapping φ of the closed positive quadrant intoitself mapping the axes into themselves with the following properties. Therestriction of φ to the open quadrant is a diffeomorphism onto its image. Thisimage is a region whose closure is a compact set bounded by intervals [0 , x ]and [0 , y ] on the x - and y -axes and two smooth curves γ and γ . The curve γ joins the point P = (0 , y ) with a point P in the positive quadrant. The curve γ joins the point P with the point P = ( x , φ in such a way that P , P and P are steady states and γ , γ and the imageof the x -axis are invariant manifolds. There is precisely one further steady state P on the boundary of the image of φ and it belongs to the interior of γ .The notation in Lemma 2 has been chosen to fit with that in the precedinganalysis. Using that analysis we obtain the following picture of the flow of thecompactified system. Within the closure of the image of φ the points P , P and P are saddles while P is a sink. The flow on γ is towards P . The flowon γ is away from P , i.e. towards P and P . The flow on the horizontal axisis towards P . The centre manifold of P enters the image of φ and the flowon it is away from P . The flow on the centre manifold of P is towards P .It is conceivable that the centre manifold of P might approach P , in whichcase it would coincide with the centre manifold of P . In that case we say thatthere is an unbounded heteroclinic cycle. Note that a solution remains boundedtowards the future on any finite time interval. For it is obvious that x remainsbounded on such an interval and that y remains bounded follows by consideringthe Poincar´e compactification.It will now be shown that if a heteroclinic cycle at infinity exists it is asymp-totically stable. This is based on a study of the passage of a solution closeto the steady states belonging to the cycle. In fact it will be shown that thiscan be reduced to the study of the corresponding passages for some simplifieddynamical systems, which will be examined first. The next lemma concerns the7ase of a hyperbolic steady state. Lemma 3
Consider the dynamical system ˙ x = − ax , ˙ y = by for positive con-stants a and b . The solution which starts at the point (1 , y ) reaches the line y = 1 at the point ( x ,
1) with x = ( y ) ab . Proof
The dynamical system can be solved explicitly. Substituting the initialand final conditions into the solution gives an expression for the time T requiredfor the solution to go from the first to the second point. Substituting theexpression for T back into the solution gives the final result. (cid:4) Treating a non-hyperbolic steady state requires the following more compli-cated statement.
Lemma 4
Consider the dynamical system˙ x = − ax (1 + ǫr ( y, ǫ )) , (21)˙ y = by k (1 + ǫr ( y, ǫ )) , (22)where a > b > k ≥ r , r are bounded. Suppose that the solutionwhich starts at the point (1 , y ) reaches the line y = 1 at the point ( x , η > ǫ and y sufficiently small the inequalities x ≤ exp( − [(1 − η ) a/ ( b ( k − y − k ) and y ≤ C ( − log x ) − k hold for a constant C . Proof
Suppose that the solution crosses the lines x = 1 and y = 1 for t = 0and t = T , respectively. Given δ > ǫ > ǫ ≤ ǫ implies | ǫr | ≤ δ and | ǫr | ≤ δ . It follows that x ≤ e − a (1 − δ ) T and y − k ≤ − b (1 − k )(1 + δ ) T . Hence T ≥ b ( k − δ ) ( y − k −
1) and if y is smallenough x ≤ exp {− a (1 − δ ) b ( k − δ ) y − k } . For δ small enough − δ (1+ δ ) ≥ − η andwe get the desired inequality for x . On the other hand x ≥ e − a (1+ δ ) T and y − k ≥ − b (1 − k )(1 − δ ) T . Hence T ≤ b ( k − − δ ) ( y − k −
1) and x ≥ exp {− a (1+ δ ) b ( k − − δ ) y − k } . Hence y ≤ C ( − log x ) − k .The reduction of the case of more general dynamical systems to these explicitcases is a consequence of the following Lemma. Lemma 5 (i) Consider a two-dimensional dynamical system which has a steadystate at a point ( x ∗ , y ∗ ) which is a hyperbolic saddle, the eigenvalues of whoselinearization at that point are − a and b , where a and b are positive. Let ( x ′ , y ′ )and ( x ′ , y ′ ) be points on the stable and unstable manifolds of ( x ∗ , y ∗ ) whichare sufficiently close to ( x ∗ , y ∗ ). Fix curves β and β through ( x ′ , y ′ ) and( x ′ , y ′ ) which are transverse to the stable and unstable manifolds respectively.Let ( x , y ) be a point on β sufficiently close to, but not equal to, ( x ′ , y ′ ). Bythe Grobman-Hartman theorem this solution intersects β and let ( x , y ) bethe first such point it reaches. Let s and s be any smooth coordinates on β and β which vanish on the corresponding invariant manifolds and increase inthe direction of the points ( x , y ) and ( x , y ). Then for s sufficiently smallthe inequality s ≤ C ( s ) ab holds for a positive constant C .(ii) Consider a two-dimensional dynamical system which has a steady state at apoint ( x ∗ , y ∗ ) whose centre and stable manifolds have dimension one. Supposethat the non-zero eigenvalue is equal to − a , the flow on the centre manifoldis away from ( x ∗ , y ∗ ) and the leading term in the evolution equation along the8entre manifold is bs k for a coordinate s on that manifold. Define points, curvesand coordinates on the curves as in part (i) except that the unstable manifoldis replaced by the centre manifold and the Grobman-Hartman theorem by thereduction theorem. Then for any η > s ≤ C exp( − [(1 − η ) a/b ( k − s − k ) and s ≤ C ( − log s ) − k hold for s sufficiently small and apositive constant C . Proof (i) It follows from Sternberg’s theorem [17] that there is a smooth dif-feomorphism which transforms the situation described to the model situationin Lemma 3 and this implies the desired statement. Note in particular that asuitable scaling of the coordinates ensures that after the transformations thedistance between the steady state and the points where the transverse sectionscross the axis is one. Sternberg’s theorem requires the hypothesis that there areno resonances. In the given case a resonance would mean an equation of theform b = − na or a = − nb for a positive integer n , which is clearly impossible.(ii) In this case we must replace Sternberg’s theorem by Takens’ theorem [18],allowing us to reduce the original system to the model system by a diffeomor-phism of arbitrarily high finite differentiability. When there is only one non-zeroeigenvalue the condition that there are no resonances is trivially fulfilled. In thiscase scaling the coordinate y allows ǫ to be made arbitrarily small. Lemma 6
When for a given value of the parameter α in the Selkov system aheteroclinic cycle at infinity exists this cycle is stable. Proof
Consider curves β i with the following properties. Let β be a curvetransverse to the unstable manifold of P and sufficiently close to P . Similarlylet β be transverse to the stable manifold of P and close to P , β transverse tothe centre manifold of P and close to P , β transverse to the centre manifoldof P and close to P , β transverse to the unstable manifold of P and close to P and β transverse to the stable manifold of P and close to P . For each i let s i be a coordinate on β i which is zero on the relevant invariant manifold andincreases towards the interior of the image of φ . Each solution which crosses β for a sufficiently small positive value of s crosses each β i at some parametervalue s i = f i − ( s i − ). It then crosses β again at some value s = f ( s ). Theaim is to obtain an estimate for the composition f of the f i which shows thatthe iterates of f tend to zero. The simplest mappings to estimate are f , f and f . These maps are smooth at zero and so by Taylor’s theorem there is aconstant C so that f ( s ) ≤ Cs , f ( s ) ≤ Cs and f ( s ) ≤ Cs . It is alsorelatively easy to estimate f since P is hyperbolic. Using part (i) of Lemma 5we get the estimate s ≤ Cs γ − . The mapping f can be estimated using part(ii) of Lemma 5 with the result that s ≤ Ce − cs − γ for a constant c . A similarestimate holds for s . To estimate f we use the second conclusion of Lemma 4.This gives s ≤ C ( − log s ) − γ − . Composing f ◦ f ◦ f and simplifying leadsto an estimate of the form s ≤ Cs γγ − . Putting all these estimates togethergives an inequality of the form f ( s ) ≤ Cs γ . Hence if s is small the distance tothe cycle is reduced by a fixed factor with each return and the cycle is stable.9 U U U N , − N , − N , + N , + Figure 2: Nullclines.
To understand the global phase portrait of the Selkov system it is helpful toconsider the nullclines N and N which are given by the equations xy γ = 1and xy γ − = 1, respectively. The relevant geometry is illustrated in Figure 2.The complement U of N ∩ N has four connected components which can bedistinguished by the signs of the right hand sides of the evolution equations. Wedenote them by U , U , U and U for the combinations of signs (+ , − ), (+ , +),( − , +) and ( − , − ) respectively. The set N \ N has two connected components N , + and N , − which are distinguished by the sign of x −
1. Similarly N \ N two connected components N , + and N , − . A solution which starts at a pointof U for t = t must either stay in U for t ≥ t or it must reach a point of N , + after a finite time, after which it immediately enters U . If it stays in U then either it is bounded and then it converges to (1 ,
1) for t → ∞ or it isunbounded and in the latter case lim t →∞ x ( t ) = ∞ and lim t →∞ y ( t ) = 0. Asolution which starts in U for t = t must reach a point of N , + after a finitetime, after which it immediately enters U . If a solution which starts in U for t = t were unbounded while remaining in U then it would approach P inthe Poincar´e compactification for t → ∞ and this has already been ruled outin the last section. Thus it must either converge to (1 ,
1) as t → ∞ or reach N , − after a finite time, after which it immediately enters U . A solution whichstarts in U for t = t must reach a point of N , − after a finite time, afterwhich it immediately enters U . Putting all this information together we getthe following result Lemma 7
Each solution of the Selkov system has one of the following behavioursin the future.(i) it converges to (1 ,
1) for t → ∞ while staying in U or U (ii) lim t →∞ x ( t ) = ∞ and lim t →∞ y ( t ) = 0(iii) it cycles indefinitely between the regions U i .A similar analysis can be done in the past time direction. Note that it is notpossible that a solution is always in U before a certain time because no solution10an approach the curve γ in the Poincar´e compactification towards the past.Thus we have the following result Lemma 8
Each solution of the Selkov system has one of the following behavioursin the past.(i) it converges to (1 ,
1) for t → −∞ while staying in U or U (ii) it starts on the y -axis(iii) lim t →−∞ x ( t ) = 0 and lim t →−∞ y ( t ) = ∞ (iv) it cycles indefinitely between the regions U i .If a solution cycles indefinitely in the future then it passes through pointsof the form (1 , y ) with y < t i . Let y i = y ( t i ). It follows from Poincar´e-Bendixson theory that the sequence y i ismonotone. If the sequence is constant then the solution is periodic. Otherwisethe sequence is strictly monotone. If it is monotone decreasing then it must tendto some y − , ∗ >
0. That this number is strictly positive follows from the factthat in the Poincar´e compactification all solutions which start with x = 1 and y sufficiently small converge to P for t → ∞ . Poincar´e-Bendixson theory furthershows that either (1 , y − , ∗ ) lies on a periodic solution or on a heteroclinic cycleat infinity. If the sequence y n is monotone increasing then y − , ∗ lies on a periodicsolution. This discussion tells us that to classify the asymptotic behaviour ofall solutions of the Selkov system for given values of γ and α it suffices to knowhow many periodic solutions there are, whether there is a heteroclinic cycleat infinity and what the stability properties of the periodic solutions and theheteroclinic cycle are. Note that every periodic solution must contain a pointof the form (1 , y ) with y < x gets large and y gets small at late times and which havean asymptotic behaviour which can be determined. This is the content of thefollowing theorem. These are exactly the solutions which belong to Case (ii) inLemma 7. Theorem 1
There exists a positive number ǫ > x (0) = x and y (0) = y which satisfies x > ǫ − and x y γ < ǫ has the late-time asymptotics x ( τ ) = τ (1 + o (1)) , (23) y ( τ ) = y e − ατ (1 + o (1)) . (24)for a constant y . There exists a solution, unique up to time translation, whichhas the asymptotic behaviour x ( τ ) = τ (1 + o (1)) , (25) y ( τ ) = τ − γ − (1 + o (1)) . (26) Proof
Any solution which starts close to the point P converges to that pointas t → ∞ . Using this information in the evolution equations for ¯ y and ¯ z allowsthem to be integrated in leading order. First we have d ¯ zds = − ¯ z γ (1 + o (1)) (27)11nd hence ¯ z ( s ) = [( γ − s ] − γ − (1 + o (1)) . (28)Putting this information back into the evolution equation allows the error term o (1) to be improved to O ( s − γ − ). The relation d ¯ yds = ¯ y ( − ( γ − α + o (1)) (29)shows that ¯ y decays exponentially as t → ∞ . Using this fact and substitutingthe improved version of (28) into the evolution equation for ¯ y gives d ¯ yds = ¯ y (cid:20) − ( γ − α + 1 γ − s − + q ( s ) (cid:21) (30)where the function q is integrable on intervals of the form [ s , ∞ ). Integratingthis last relation gives ¯ y ( s ) = As γ − e − α ( γ − s (1 + o (1)) (31)for a constant A . It follows from the estimates obtained up to now that theimage of each such solution is tangent to the ¯ z -axis at P and is therefore acentre manifold at P . It remains to transform these formulae from the variables(¯ y, ¯ z, s ) to the variables ( x, y, τ ). Note that dsdt = γ − Z γ and dτdt = Z γ . Thecoordinate s is only defined up to an additive constant and if this constant ischosen appropriately it follows that τ = ( γ − s . Using the expressions (27) and(29) and the definitions of the coordinate transformations gives the first partof the theorem. The solution mentioned in the second part of the theorem is asolution on the centre manifold of P . Integrating the equation for ¯ y in Lemma1 gives ¯ y ( s ) = ( ν ( γ − s ) − γ − (1 + o (1)) = ( γs ) − γ − (1 + o (1)). Putting thisinto the equation for w given there leads to w ( s ) = ν ( γs ) − (1 + o (1)). Now dsdt = γ − ( γs ) − γ (1 + o (1)). On the other hand dτdt = Z γ = ( γs ) − γ (1 + o (1)).Thus we can choose s so that τ = γs . Using the definitions of the coordinatetransformation gives the second part of the theorem. (cid:4) To obtain more information about periodic solutions it is helpful to do thefollowing changes of variables. Introduce a time variable by dtdτ = y γ and define v = dydt . Let f ( y ) = 1 − α ( γ − y − γ and g ( y ) = α ( y − y − γ . Denote by F theprimitive of f which vanishes for y = 1. Introduce variables by ˜ x = y − y = − v − F ( y ). Dropping the tildes gives the system˙ x = − y − F ( x ) , (32)˙ y = g ( x ) . (33)We are interested in solutions of this for which x > − y . It will be shown that under suitable circumstances thesystem (32)-(33) defined on a region of the form ( x , ∞ ) × R has at most oneperiodic solution and that if such a solution exists it is asymptotically stable.12his implies a corresponding result for the Selkov system. Next we discuss someauxiliary results used in the proof of these statements. Theorem 2
Suppose that the system (32)-(33) satisfies the following conditions:(i) xg ( x ) > x = 0 and g ′ (0) > F (0) = 0 and f (0) < f ( x ) = F ′ ( x )(iii) There exists a real number α such that the function f ( x ) = f ( x ) + αg ( x )has zeroes x < x > f ( x ) ≤ x < x < x and there exists x ∗ with f ( x ∗ ) > x ∗ ,
0) in the interior of L .(iv) All the limit cycles are contained in the interval x ≤ x ≤ x , where x < x < < x < x and the function f /g does not decrease for x ≤ x ≤ x and x ≤ x ≤ x .(v) all the limit cycles contain the interval x ≤ x ≤ x on the x -axis.Then the system has at most one limit cycle and if one exists it is stable.This theorem is closely related to a result stated in Kuang and Freedman [13],which is in turn related to a result of Cherkas and Zhilevich [2]. The differencesbetween Theorem 2 and Theorem 3.1 in [13], apart from the different notation,are as follows. For simplicity the function φ ( y ) in [13] has been set to y and theconstant β has been set to zero. The condition (4) of the theorem of [13] hasbeen replaced by the condition (4) ′ . The condition g ′ (0) > x and x are simple zeroes have been replaced by the assumption concerning the existenceof the point x ∗ . This is because we could not see how to prove that Theorem3.2 of [13] follows from Theorem 3.1, a claim made without proof in the paper. Proof of Theorem 2
We give only an outline of the proof. More details can befound in [13] and [1]. The unique steady state of the system is at the origin. Itfollows from the condition f (0) < L and L . It is necessarily the case that the steady state is in the interior of theJordan curves defined by the solutions and that one of them is in the interiorof the other. Suppose that L is in the interior of L . To examine the stabilityof the solutions the Poincar´e stability criterion [15] will be used. It involves thequantity obtained by integrating the divergence of the vector field defining thedynamical system over the solution. We call this the Poincar´e quantity. Let A and A be the points on L and L with x = x and y > D and D bethe points with with x = x and y <
0. Let E and E be the points on L withthe same y coordinates as A and D . Integrate the differential form ( f /g ) dy overthe boundary of the region R bounded by the line segments AE and DE andthe parts of L and L joining A with D and E with E . This results in the sumof two integrals along the curved parts of the boundary. On the other hand, byStokes’ theorem this integral is equal to the integral of a non-negative quantityover R . Thus we get an inequality relating the integrals of f along certainparts of L and L . We can do a similar construction starting from x = x with points B, B , C, C , K, K corresponding to A, A , D, D , E, E . Next wereplace ( f /g ) dy in this construction by f / ( y + F ( x )) dx and integrate overregions with corners C, D, C , D and A, B, A , B . Putting these computationstogether shows that the integral of f over L is equal to its integral over L L joining E with A , B with K , K with C and D with E . In fact these extra contributions are non-negative. Thisis because f is non-negative on the relevant intervals. Moreover it is positiveon the interval [ x ∗ , x ], due to the monotonicity of f /g . It follows that theintegral of f over L is strictly greater than its integral over L . It can beshown that the integrals of g over L and L are zero so that the integral of f over L is strictly greater than its integral over L . The divergence of thevector field defining the system is equal to − f and so the Poincar´e quantity for L is less than that for L . Since g ′ (0) > t → ∞ . As a consequence it is a limit cycle and we choose it as L . It followsthat L is stable from the inside and that the Poincar´e quantity is non-positiveand hence the integral of − f over L also non-positive. Hence the Poincar´equantity for L must be negative. If L were stable on the outside we couldchoose L to be the closest limit cycle outside it. L would have to be unstableon the inside, a contradiction. Thus L is semistable (stable on one side andunstable on the other). It can be shown by perturbing the system using themethod of rotated vector fields [6] that if this were the case there would be asystem where the above contradiction is obtained. It follows that no secondlimit cycle can exist and if one limit cycle does exist it is stable. (cid:4) This theorem can be used to prove another where the opaque conditioninvolving the function f is removed from the hypotheses. It is a modificationof a theorem of Zhang [19]. Theorem 3
Suppose that the system (32)-(33) satisfies the conditions (i) and(ii) of Theorem 2 and that all limit cycles are contained in the interval a < x < b where a < < b and f ( x ) /g ( x ) is non-decreasing when x increases in a < x < < x < b . Then the conclusion of Theorem 2 holds. Proof
The condition that f /g is non-decreasing is equivalent to the conditionthat f /g is non-decreasing. The requirement on the region where this functionis non-decreasing in the hypotheses of Theorem 3 implies that in the hypothesesof Theorem 2. We need to show that α can be chosen so that f has zeroes withthe properties required in condition (iii) of Theorem 2. Following [19] let ˜ x bethe smallest value of x attained at a point of L and define α = − f (˜ x ) g (˜ x ) . Then f (˜ x ) is zero by construction and ˜ x is a zero of f with ˜ x <
0. Let x be thelargest value of x for which f ( x ) = 0 for ˜ x ≤ x ≤ x . Since f (0) < x <
0. We will now show that f has a second zero x > x whichis in the interior of L . If this were not the case then L would lie completely inthe region where f is non-positive, as will now be shown. For f ( x ) g ( x ) = 0 andwe have assumed that ddx (cid:16) f ( x ) g ( x ) (cid:17) is non-negative on the interval [ x , f ( x ) g ( x ) ≥ f ( x ) ≤ f ( x ) g ( x ) < f ( x ) g ( x ) is negative on an interval beginning at x and extending beyond x = 0.This establishes the claim concerning L . It can be concluded that the integralof f over L is negative. For under the given assumptions f is non-positive14verywhere on L and negative at some point of L . The Poincar´e criterionthen implies that L is unstable in the interior. This contradicts the fact thatthere are solutions starting near the steady state which converge to L . Thusin reality f has a zero x > L . In fact f must becomepositive for some x ∗ > x ∗ ,
0) in the interior of L since otherwise thecontradiction would still occur. Choosing x to be the smallest x > f is negative on (0 , x ) ensures that the hypotheses of Theorem 3 are satisfied. (cid:4) Theorem 3 can be used to obtain a result about the Selkov system followingd’Onofrio.
Theorem 4 [3] If the steady state in the Selkov model is unstable for a givenvalue of α then there is at most one limit cycle and if one exists it is asymptot-ically stable.In his proof of this theorem d’Onofrio cites the Theorem of Kuang andFreedman. However the theorem he formulates (Proposition 1 of his paper)is not equivalent to the corresponding result stated by Kuang and Freedman(Theorem 3.2 of their paper). In [13] the assumptions include the conditionthat φ ( y ) + F ( x ) is defined for all x ∈ ( −∞ , ∞ ). In the case φ ( y ) = y of interesthere this means that F would have to be defined on the whole real line. Thiscondition is not assumed in d’Onofrio’s Proposition 1 and indeed it does nothold in the situation of the application to the Selkov model. It is this apparentgap which motivates our discussion of Theorem 2 and Theorem 3 above. Thatdiscussion shows that the extra assumption is not necessary. When this hasbeen clarified a calculation done by d’Onofrio showing the monotonicity of f /g in the case of the Selkov model suffices to show that Theorem 4 follows fromTheorem 3.Theorem 4 can be used to obtain information about the long-time behaviourof solutions in the case that the steady state is unstable. Using Poincar´e-Bendixson theory it follows that the ω -limit set in the Poincar´e compactificationof any solution which cycles indefinitely in the future (case (iii) of Lemma 7)is either the unique periodic solution or the heteroclinic cycle at infinity. Weknow as a result of Lemma 6 and Theorem 4 that both the periodic solutionand the heteroclinic cycle are stable whenever they exist. If both existed thenby continuity there would have to be a periodic solution between them and thiswould contradict Theorem 4. Thus we obtain the following result Theorem 5
In the case α > α exactly one of the following three situationsoccurs.(i) The centre manifolds of P and P coincide so that there is a heterocliniccycle at infinity. Any solution which starts below this centre manifold convergesto P as t → ∞ while any solution other than the steady state which startsabove this manifold converges to the heteroclinic cycle at infinity as t → ∞ .(ii) The centre manifolds of P and P do not coincide. There exists a uniqueperiodic solution. Any solution which starts below the centre manifold of P converges to P as t → ∞ while any solution other than the steady state whichstarts above this manifold converges to the periodic solution as t → ∞ .(iii) The centre manifolds of P and P do not coincide. Any solution other15han the steady state which does not lie on the centre manifold of P convergesto P as t → ∞ .We next turn to the case α ≤ α . For this case d’Onofrio proved a result(Proposition 5.3 of his paper) as an application of a Theorem of Hwang andTsai [9]. It was shown in [1] that one of the conditions in his result is satisfiedautomatically for α ≤ α , thus showing that the result can be generalized. Wedo so in Theorem 7 below. The following theorem is closely related to Theorem2.1 of [9]. Theorem 6
Suppose that in the system (32)-(33) defined on the region ( r , r ) × R the following conditions hold:(i) there exists λ ∈ ( r , r ) such that g ′ ( λ ) > x − λ ) g ( x ) > x ∈ ( r , r ) \ { λ } .(ii) there exist constants a, b ∈ R such that the function f ( x )+ ag ( x )+ bg ( x ) F ( x )is non-negative on ( r , r ) and not identically zero.Then the system has no periodic solutions.The relation of this statement to that of Theorem 2.1 of [9] is that y hasbeen replaced by − y , the function φ has been taken to be equal to one andthe function π ( y ) equal to y . Under these circumstances conditions (A1), (A2)and (A4) of the theorem of [9] are satisfied automatically and the remainingconditions reduce to (i) and (ii) above. In the case of the Selkov system theinterval can be chosen to be ( − , ∞ ) and condition (i) is satisfied by λ = 0.Thus it only remains to consider the condition (ii) and in fact it suffices tochoose b = 0.Substituting the values of f and g in the case of the Selkov system into thepositivity condition in (ii) gives( x + 1) γ + aα ( x + 1) + α (1 − γ − a ) ≥ . (34)Let β = α ( γ − β ≤
1. Choose a = − γ ( γ − x + 1) γ + aα ( x + 1) + α (1 − γ − a )= ( x + 1) γ − αγ ( γ − x + 1) + α ( γ − (35)= ( x + 1) γ − βγ ( x + 1) + β ( γ − . (36)The last function has a global minimum in (0 , ∞ ). This occurs at the point x + 1 = β γ − . Substituting this back into the function shows that the minimumis ( γ − β (1 − β γ − ) ≥
0. Thus this choice of a has the desired property andTheorem 6 can be applied to the Selkov system. Theorem 7
If the positive steady state of the Selkov sytem is stable then thereare no periodic orbits and every solution either tends to the positive steady stateor tends to P or P at late times. Proof
The first statement follows from Theorem 6. There cannot be a hetero-clinic cycle at infinity since the fact that both the cycle and the steady stateare stable would mean that there would have to be a periodic solution betweenthem. The only remaining possibilities for the ω -limit set of a solution are thoselisted in the theorem. 16 Further remarks
The main question which has been left open by the analysis in this paper iswhether there are parameter values for which there are solutions of the Selkovmodel which exhibit unbounded oscillations. This is the question, whether case(i) in Theorem 5 ever occurs. To prove this it would suffice to show that thereis a value of α for which cases (ii) and (iii) of Theorem 5 are ruled out. To ruleout case (ii) for a given value of α it would be enough to show that the systemadmits a Lyapunov function or a Dulac function.The model (1)-(2) which we have studied in this paper was originally derivedby Selkov from a five-dimensional model, the system (4) in [16], by a singularlimiting process. This limit will now be considered. The variables s , s , x , x and e in the larger model are the concentrations of the substrate ATP, theproduct ADP, the activated enzyme, the complex between the activated enzymeand the substrate and the inactive enzyme, respectively. With some assump-tions about the nature of these reactions and assuming mass action kinetics thefollowing system of equations is obtained. ds dt = v − k s x + k − x , (37) ds dt = k x − γk s γ e + γk − x − ks , (38) dx dt = − k s x + ( k − + k ) x + k s γ e − k − x , (39) dx dt = k s x − ( k − + k ) x , (40) dedt = − k s γ e + k − x . (41)Note that e = e + x + x is a conserved quantity (total amount of enzyme)and this can be used to eliminate e from the first four evolution equations anddiscard the evolution equation for e . This reduces the system to four equations.Next dimensionless variables are introduced by defining σ = k s k − + k s , σ = (cid:16) k k − (cid:17) γ s , u = x e , u = x e , θ = e k k k − + k t . This leads to the system dσ dθ = ν − k + k − k u σ + k − k u , (42) dσ dθ = η (cid:18) u − γ k − k σ γ (1 − u − u ) + γ k − k u − χσ (cid:19) , (43) ǫ du dθ = u − σ u + k − k + k − ( σ γ (1 − u − u ) − u ) , (44) ǫ du dθ = σ u − u . (45)Explicit expressions for the parameters ǫ , ν , η and χ can be found in [1] or[10]. Formally setting ǫ = 0 in the equations (44) and (45) gives u = σ u and17 = σ γ σ γ + σ σ γ and substituting these relations into the evolution equationsfor σ and σ gives dσ dθ = ν − (cid:18) σ σ γ σ γ + σ σ γ (cid:19) , (46) dσ dθ = η (cid:18) σ σ γ σ γ + σ σ γ − χσ (cid:19) . (47)It can be shown that solutions of this system of two equations can be approx-imated by solutions of the four-dimensional system (42)-(45) using geometricsingular perturbation theory [12]. To do this we need to examine the transverseeigenvalues. This means computing the matrix of partial derivatives of the righthand sides of the evolution equations for u and u with respect to the variables u and u and evaluating the result for ǫ = 0. A computation shows (for detailssee [1]) that this matrix has determinant k − k + k − ( σ σ γ + σ γ + 1) > − σ − k − k + k − ( σ γ + 1) − <
0. Hence both eigenvalues have negative real partsand the limit is well-behaved.Selkov claims that starting from the assumptions of Higgins and supposingon biological grounds that certain quantities are small leads to system (46)-(47)with γ = 1. This system has the unfortunate property that it does not admitperiodic solutions. This can be proved using the fact that σ γ + σ σ γ σ σ is a Dulacfunction.Consider now the additional rescaling x = ν γ − χ γ σ , y = χν σ and τ = (cid:16) νχ (cid:17) γ θ .In the equations for dxdτ and dydτ with parameters ν , χ and η we can replace η by α = ηχ γ +1 ν γ . Then the limit ν → χ is of orderone, ν is small and η is small. Geometric singular perturbation theory thenallows us to conclude that the system (42)-(45) possesses a periodic solutionwhich is hyperbolic and stable for ǫ small.Note that there is an alternative mathematical model of glycolysis due toGoldbeter and Lefever [7] which appears to explain some of the experimentalobservations better than that of Selkov. For a discussion of this issue we referto [10]. The model of [7] exhibits intricate dynamical behaviour which has beeninvestigated using advanced methods of geometric singular perturbation theoryin [11]. The Goldbeter-Lefever model was also studied mathematically in [4]. Acknowledgement
We thank Hussein Obeid for drawing our attention to thepossible relevance of quasihomogeneous directional blow-ups for this problem.18 eferences [1] Brechmann, P. 2017 Der Higgins-Selkov-Oszillator - ein mathematischesModell der Glykolyse. MSc thesis, Johannes Gutenberg-Universit¨at Mainz.[2] Cherkas, L. A. and Zhilevich, L. I. 1970 Some tests for the absence oruniqueness of limit cycles. Diff. Eq. 6, 891–897.[3] d’Onofrio, A. 2010 Uniqueness and global attractivity of glycolytic oscilla-tions suggested by Selkov’s model. J. Math. Chem. 48, 339-346.[4] d’Onofrio, A. 2011 Globally attractive oscillations in open monosubstrateallosteric enzyme reactions. J. Math. Chem. 49, 531-545.[5] Dumortier, F., Llibre, J. and Art´es, J. C. 2006 Qualitative theory of planardifferential systems. Springer, Berlin.[6] Duff, G. F. D. 1953 Limit-cycles and rotated vector fields. Ann. Math. 57,15–31.[7] Goldbeter, A. and Lefever, R. 1972 Dissipative structures for an allostericmodel. Biophys. J. 12, 1302-1315.[8] Higgins, J. 1964 A chemical mechanism for oscillation of glycolytic inter-mediates in yeast cells. Proc. Natl. Acad. Sci. (USA) 51, 989-994.[9] Hwang, T.-W. and Tsai, H.-J. 2005 Uniqueness of limit cycles in theoreticalmodels of certain oscillating chemical reactions. J. Phys. A 38, 8211-8223.[10] Keener, J. and Sneyd, J. 2009 Mathematical physiology. I: Cellular Physi-ology. Springer, Berlin.[11] Kosiuk, I. and Szmolyan, P. 2011 Scaling in singular perturbation problems:blowing up a relaxation oscillator. SIAM J. Appl. Dyn. Sys. 10, 1307–1343.[12] Kuehn, C. 2015 Multiple time scale dynamics. Springer, Berlin.[13] Kuang, Y. and Freedman, H. I. 1988 Uniqueness of limit cycles in Gause-type models of predator-prey systems. Math. Biosci. 88, 67-84.[14] Kuznetsov, J. A. 1995 Elements of Applied Bifurcation Theory. Springer,Berlin.[15] Perko, L. 2001 Differential Equations and Dynamical Systems. Springer,Berlin.[16] Selkov, E. E. 1968 Self-oscillations in glycolysis. I. A simple kinetic model.Eur. J. Biochem. 4, 79–86.[17] Sternberg, S. 1958 On the structure of local homeomorphisms of Euclidean nn