Unbounded solutions of models for glycolysis
aa r X i v : . [ q - b i o . M N ] M a r Unbounded solutions of models for glycolysis
Pia Brechmann and Alan D. RendallInstitut f¨ur MathematikJohannes Gutenberg-Universit¨atStaudingerweg 9D-55099 MainzGermany
Abstract
The Selkov oscillator, a simple description of glycolysis, is a system oftwo ordinary differential equations with mass action kinetics. In previouswork the authors established several properties of the solutions of thissystem. In the present paper we extend this to prove that this system hassolutions which diverge to infinity in an oscillatory manner at late times.This system was originally derived from another system with Michaelis-Menten kinetics. It is shown that the Michaelis-Menten system, like thatwith mass action, has solutions which diverge to infinity in a monotonemanner. It is also shown to admit subcritical Hopf bifurcations and thusunstable periodic solutions. We discuss to what extent the unboundedsolutions cast doubt on the biological relevance of the Selkov oscillatorand compare it with other models in the literature.
When trying to understand a biological system with the help of mathematicalmodelling it often happens that there are several different models for the samebiological situation in the literature. In view of this it is important to havecriteria for deciding between models. One strategy for identifying criteria ofthis type is to look at relatively simple examples in great detail. In order to dothis effectively it is necessary to have a sufficiently comprehensive understandingof the properties of solutions of the models being studied. In this paper, withthis strategy in mind, we look in detail at the dynamical properties of certainmodels for glycolysis.Glycolysis is part of the process by which living organisms extract energyfrom sugar [1]. A suitable model system for studying this phenomenon experi-mentally is yeast extract or a suspension of yeast cells. The first indication thatthis system might have interesting dynamical properties was given by dampedoscillations reported in [7]. Later it was discovered that a constant continuoussupply of sugar can lead to sustained oscillations (cf. [2]). Looking for the1ource of these oscillations revealed that they are produced by a small reactionnetwork describing the action of the enzyme phosphofructokinase. A mathe-matical model for this network was set up and studied by Higgins [11]. It wasfound by Selkov that this model was not adequate for describing the oscillationsand he introduced a modified one [19]. The starting point for the model ofSelkov is a reaction network with five chemical species. Assuming mass actionkinetics leads to a system of five ordinary differential equations. Using quasi-steady state assumptions this can be reduced to a system of two equations withnonlinearities of Michaelis-Menten type. For brevity we call it ’the Michaelis-Menten system’ in what follows. Setting one of the coefficients in this systemto zero leads to a further simplification, giving a system of two equations withmass 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 ofsolutions of the three systems just described. A number of properties of solutionsof the basic Selkov system were already established in [19] but for many yearsno further rigorous results on this subject were obtained. Important progresswas made in a paper of d’Onofrio [5] and a number of additional properties ofthe solutions were established in a recent paper of the authors [3]. In particularit was proved that for any values of the parameters there exist unboundedsolutions of this system which are eventually monotone in the sense that fora solution of this type both concentrations are monotone after a certain time.In [19] it is claimed that this system admits solutions which oscillate with anamplitude which grows without limit at late times. In what follows solutions ofthis type are referred to as ’solutions with unbounded oscillations’. The paper[19] provides no justification for the claim other than a mention of numericalsimulations, about which no details are given. Up to now there was no proofof the truth or falsity of this claim of [19]. One of the main results of thepresent paper is a proof of the existence of solutions of the basic Selkov systemwith unbounded oscillations. Our discovery of this proof was stimulated by thepaper [16], which belongs to the domain of theoretical chemistry. It deals witha system which turns out to be identical to the basic Selkov system when aparameter γ in the latter system takes the value two.In [16] a claim of the existence of solutions with unbounded oscillations is alsomade. It is supported by an intricate heuristic argument using matched asymp-totic expansions. It is not at all clear how this argument could be translated intoa rigorous one but it provided us with some ideas which, when combined withthe results of [3], do give a proof of the existence of solutions with unboundedoscillations. When written in dimensionless form the system contains one pa-rameter α . As claimed in [16], solutions with unbounded oscillations occur forprecisely one value α of α . When α is slightly less than α there exists a stableperiodic solution. As α approaches α from below the amplitude of the periodicsolution tends to infinity. One important element of this proof is to study thelimit of the system for α → ∞ after a suitable rescaling. The existence of α is then proved by a shooting argument. A monotonicity property, which wasapparently not previously known, is used to obtain the uniqueness of α .The presence of unbounded solutions, whether monotone or oscillatory, might2e seen as a feature which is unrealistic from the point of view of the biologicalapplications. The monotone unbounded solutions of the basic Selkov system arenot mentioned at all in [19]. That system is the limit of the Michaelis-Mentensystem when a parameter ν tends to zero. It is stated in [19] that solutionswith unbounded oscillations do not exist for ν >
0. On the other hand simu-lations reported in [12] suggest that the amplitude of periodic solutions of theMichaelis-Menten system diverges rapidly to infinity when a parameter is variedin a finite range. This indicates that, in contrast to the claim of Selkov, the ex-istence of unbounded oscillations is a phenomenon which may persist for ν > ν →
0. Theissue of the existence of solutions with unbounded oscillations in the case of theMichaelis-Menten system is not resolved in what follows but some partial resultsare obtained. In particular it is shown that for the Michaelis-Menten systemwith arbitrary parameters there are unbounded solutions which are eventuallymonotone and whose leading order asymptotics are identical to those found inthe basic Selkov system. It is also shown that for certain combinations of theparameters ( α, ν ) with ν >
In [19] a simple reaction network describing glycolysis is introduced. Assumingmass action kinetics for this network leads to a system of five ordinary differential3quations, system (4) of [19]. In a slightly modified notation this system is ds dt = v − k s x + k − x , (1) ds dt = k x − γk s γ e + γk − x − v s , (2) dx dt = − k s x + ( k − + k ) x + k s γ e − k − x , (3) dx dt = k s x − ( k − + k ) x , (4) dedt = − k s γ e + k − x . (5)In fact a factor γ was omitted in two places in [19] and this error has beencorrected here. All the parameters are positive and it is assumed that γ > e = e + x + x is a conserved quantity (total amount of enzyme) and this can be used toeliminate e from the first four evolution equations and discard the evolutionequation for e . This reduces the system to four equations.Dimensionless variables can be introduced by defining σ = k s k − + k , σ = (cid:18) k k − (cid:19) γ s , u = x e , u = x e , θ = e k k k − + k t. (6)This leads to the system dσ dθ = ν − k + k − k u σ + k − k u , (7) dσ dθ = η (cid:18) u − γ k − k σ γ (1 − u − u ) + γ k − k u − χσ (cid:19) , (8) ǫ du dθ = u − σ u + k − k + k − ( σ γ (1 − u − u ) − u ) , (9) ǫ du dθ = σ u − u (10)where ǫ = e k k ( k + k − ) , ν = v k e , η = k + k − k (cid:18) k k − (cid:19) γ , χ = v k e (cid:18) k − k (cid:19) γ . (11)Formally setting ǫ = 0 in the equations (9) and (10) gives u = σ u and u = σ γ σ γ + σ σ γ and substituting these relations into the evolution equationsfor σ and σ gives dσ dθ = ν − (cid:18) σ σ γ σ γ + σ σ γ (cid:19) , (12) dσ dθ = η (cid:18) σ σ γ σ γ + σ σ γ − χσ (cid:19) . (13)4s has been discussed in [3] geometric singular perturbation theory (GSPT) canbe used to show that solutions of (7)-(10) converge to solutions of (12)-(13) inthe limit ǫ → x = ν γ − χ γ σ , y = χν σ , α = ηχ γ +1 ν γ , β = ν γ − χ γ , τ = (cid:18) νχ (cid:19) γ θ. (14)Expressing the equations (12) and (13) in terms of these gives dxdτ = 1 − xy γ νy γ ( β + x ) , (15) dydτ = α (cid:20) xy γ νy γ ( β + x ) − y (cid:21) . (16)This system has a regular limit when ν tends to zero with α and β fixed. In thelimit we get the basic Selkov system, system (II) of [19], which is dxdτ = 1 − xy γ , (17) dydτ = αy ( xy γ − − . (18)It is the system of central interest in [19] and the dynamical properties of itssolutions are studied in detail in [3]. Of course (17)-(18) can be thought of asthe special case of (15)-(16) where ν = 0. The following proposition collects some of the properties of solutions of the basicSelkov system established in [3].
Proposition 1
The basic Selkov system (17)-(18) has the following properties.1. For each value of the parameter α the unique positive steady state hascoordinates (1 , α ∈ (cid:16) , γ − (cid:17) the positive steady state is asymptotically stable andthere exist no periodic solutions.3. For α = α = γ − a generic supercritical Hopf bifurcation occurs.4. For each value of the parameter α there exist positive numbers x and y such that if a solution satisfies x ( t ) ≥ x and y ( t ) ≤ y at some time t it satisfies˙ x ( t ) > y ( t ) < t →∞ x ( t ) = ∞ and lim t →∞ y ( t ) = 0.Using the standard theory of Hopf bifurcations it follows from statement 3.of the proposition that for any α slightly greater than α there exists a stableperiodic solution. A key question left open in [3] is that of what happens to theperiodic solution when α gets large. This question is answered in this section. Theorem 1
There exists a number α > α such that the basic Selkov system(17)-(18) has the following properties. 5. For α = α 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 α < α < α there exists a unique periodic solution and it is asymptoti-cally stable.3. For α > α each solution other than the steady state is unbounded and hasthe properties described in statement 4. of Proposition 1.4. As α tends to α from below the diameter of the image of the periodicsolution tends to infinity.We adopt some of the notation of [3]. There the Poincar´e compactificationof the basic Selkov system is computed and one of the resulting points at infinityis blown up. After this has been done there are four steady states at infinitycalled P , P , P and P . Their positions can be seen in Fig. 1 of [3]. Eachof the points P and P has a one-dimensional centre manifold with the flowon the centre manifolds being away from P and towards P . As a startingpoint for the proof of Theorem 1 we establish some further properties of thecentre manifolds of the points P and P , both of which are unique. Let L bethe segment of the line y = 1 where 0 < x ≤
1. We use the notation for thecomponents U i of the complement of the nullclines which can be seen in Fig. 2of [3]. Lemma 1
In the basic Selkov system the centre manifolds of P and P bothcontain a point of L in their closures. Proof
For a point on the centre manifold of P sufficiently near to P we have˙ x >
0. Hence the manifold initially lies in the region U . As long as x < U and both coordinates of a solution on the centre manifoldare monotone functions of time. Hence a solution on the centre manifold of P either reaches a point of L with x < t → ∞ . Similarly a solution on the centre manifold of P , when followed backwards in time, either reaches a point of L with x < t → −∞ . (cid:4) For a given value of α let ξ ( α ) be the x -coordinate of the point where thecentre manifold of P meets L if such a point exists and otherwise let ξ ( α ) = 1.Define ξ ( α ) similarly in terms of the centre manifold of P . Note that eachcentre manifold depends smoothly on the parameter α , in the sense that we canchoose 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 seenby considering the suspended system obtained by adjoining the equation ˙ α = 0to the basic Selkov system and noting that it has a two-dimensional centremanifold at the points corresponding to P and P . This manifold is foliated bycurves of constant α which are centre manifolds for the original system. Theirsmooth dependence on α follows from the smoothness of the two-dimensionalcentre manifold. Lemma 2
The function ξ − ξ describing the separation of the points wherethe centre manifolds of P and P reach y = 1 is continuous. Proof
Consider a value of α c for which ξ ( α c ) <
1. The centre manifold forthat value crosses L transversely and so, by the implicit function theorem, ξ is a smooth function of α close to α c . This also shows that the set of values6f α for which ξ ( α ) < α ∗ of α for which ξ ( α ∗ ) = 1 and a sequence α n satisfying lim n →∞ α n = α ∗ . It will be shownthat lim n →∞ ξ ( α n ) = 1. Together with the information already obtained thisimplies that ξ is continuous everywhere. The desired statement will be provedby contradiction. If ξ ( α n ) did not converge to one then by passing to a sub-sequence we could assume that lim n →∞ ξ ( α n ) = ξ s <
1. Consider now thesequence of solutions of the basic Selkov system with x n (0) = ξ ( α n ), y n (0) = 1and α = α n and the solution with x s (0) = ξ s , y s (0) = 1 and α = α ∗ . Weare 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 thecentre manifold of P for α = α ∗ . If ( x s , y s ) lies to the left of the centre manifoldthen 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 thesolutions ( x n , y n ) would do the same, a contradiction. The conclusion is thatthe solution ( x s , y s ) lies on the centre manifold and hence ξ ( α ∗ ) <
1, in con-tradiction to the definition of α ∗ . It has thus been proved that ξ is continuous.A similar argument shows that ξ is continuous. Hence ξ − ξ is continuous. (cid:4) Lemma 3
The function ξ − ξ describing the separation of the points wherethe centre manifolds of P and P reach y = 1 is positive for 0 < α ≤ α andnegative for α sufficiently large. There exists an α with ξ ( α ) = ξ ( α ). Proof
Suppose that for a given value of α we have ( ξ − ξ )( α ) ≤
0. The regionof the Poincar´e compactification bounded by the parts of the centre manifoldsof P and P ending on L and the part of L between them and above the centremanifold of P is invariant under evolution backwards in time. Consider thesolution obtained by backward time evolution of a point in this region otherthan the steady state. By Poincar´e-Bendixson theory its α -limit set must bea steady state or a periodic solution. If α ≤ α this leads to a contradiction,because in that case no periodic solutions exist and the positive steady state is asink. Thus we can conclude that the function ξ − ξ is positive for 0 < α ≤ α .Next we investigate the behaviour of solutions for α large. The followingcalculations were inspired by a transformation introduced in [16] in the case γ = 2. It is given by µ = α − γ , ˜ x = α γ − γ x , ˜ y = α − γ y and ˜ τ = ατ . Theequations become d ˜ xd ˜ τ = µ − ˜ x ˜ y γ , (19) d ˜ yd ˜ τ = ˜ x ˜ y γ − ˜ y (20)and we are interested in the limit µ →
0. A Poincar´e compactification of thissystem was carried out in [16]. After a suitable rescaling this leads to a systemin the standard form of a fast-slow system in GSPT. (For background on GSPTwe refer to [14].) Unfortunately in this system the important propery of normalhyperbolicity breaks down at the point corresponding to P . It turns out thatthis problem can be got around by using the transformations introduced in [3]7o treat the behaviour of solutions for x large. These can be summed up bydefining ¯ y = x − γ y γ and ¯ z = x − γ y − γ − γ and choosing a time coordinate s satisfying dsdτ = γ xy γ − . This transforms the basic Selkov system into system(12)-(13) of [3]. Now introduce ǫ = α − and ¯ w = α (¯ z − s by a prime, we get the system¯ y ′ = − γ ¯ y ¯ w − ¯ yǫ − [(1 + ǫ ¯ w ) γ − − γǫ ¯ w ] + ¯ y γ +1 − ¯ y γ (1 + ǫ ¯ w ) γ +1 , (21) ǫ ¯ w ′ = γ ( γ −
1) ¯ w − ( γ − y γ (1 + ǫ ¯ w )+( γ − ǫ − [(1 + ǫ ¯ w ) γ +1 − − ǫ ( γ + 1) ¯ w ] − ¯ y γ − (1 + ǫ ¯ w ) γ +2 + γ ¯ y γ (1 + ǫ ¯ w ) . (22)Note that, due to cancellations in the expressions in square brackets this systemis regular at ǫ = 0 and in fact the apparently singular term even vanishes as ǫ →
0. The critical manifold has the equation γ ( γ −
1) ¯ w = ¯ y γ − − ¯ y γ = ¯ y γ − (1 − ¯ y ).The derivative of the right hand side of the equation (22) with respect to ¯ w ,evaluated at ǫ = 0, is γ ( γ − d ¯ yds = − γγ − y γ (1 − ¯ y ) . (23)On the critical manifold there are two steady states, a source and a sink. Theyare connected by a heteroclinic orbit. For ǫ small and positive the criticalmanifold perturbs to a one-dimensional invariant manifold. All steady stateswhich exist must lie on that manifold. The steady state at ¯ y = 1 is hyperbolicand so perturbs to a hyperbolic source. The steady state at ¯ y = 0 continues toexist and there are no others. It follows that there is also a connection betweenthe positive steady state of the Selkov system and the point P on the boundaryfor α sufficiently large. In other words, when α is sufficiently large the centremanifold of P converges to the positive steady state in the past. This meansthat ξ ( α ) = 1. On the other hand, since the positive steady state is a source inthis case the centre manifold of P cannot converge to the positive steady state.We conclude that ξ ( α ) < ξ − ξ is negative. By the intermediatevalue theorem there exists some α with ξ ( α ) = ξ ( α ). Note that in the endthe equations (19)-(20) were not needed in the proof but we judged it useful toinclude them so as to give an indication of how the argument was found. (cid:4) It turns out that the value of α for which the centre manifolds of P and P meet is unique. This follows from a monotonicity property of the dependenceof the centre manifolds on α . Lemma 4
The function ξ − ξ describing the separation of the points wherethe centre manifolds of P and P reach y = 1 is strictly decreasing and has aunique zero. Proof
For this proof it is convenient to think of y as a function of x for a givensolution. Suppose that a solution y ( x ) for a parameter α crosses a solution ˆ y ( x )for a parameter ˆ α < α . Then the (negative) slope of ˆ y is smaller in magnitude8han that of y . Thus if ˆ y is larger than y for some x it must remain so forall larger x . Similarly, if ˆ y is smaller than y for some x it must remain so forall smaller x . The leading order approximation to the centre manifold of P isgiven by ¯ z = 1 + ν ¯ y γ − + . . . where ν = αγ ( γ − . This translates (in terms ofvariables used in [3]) to Z = Y γ − γ + ν Y γ − γ . . . and x = y − γ − γν + . . . .Putting these things together shows that when α is reduced the intersectionof the centre manifold of P with the line y = 1 moves to the left. To obtaininformation about the position of the centre manifold of P in its dependenceon α it is necessary to determine one more order in the expansion of the centremanifold than was done in [3]. The result is X = Z γ +1 − γαZ γ +1 + . . . . Inthe original variables this gives x = y − γ − γαy − γ + . . . . 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 centremanifold of P with the line y = 1 moves to the right. This implies that thefunction ξ − ξ is strictly decreasing and cannot have more than one zero. (cid:4) Proof of Theorem 1
By Lemma 3 and Lemma 4 there exists a unique α > α for which the centre manifolds of P and P coincide. With this informationthe first statement of Theorem 1 follows immediately from the first statementof Theorem 3 of [3]. For α < α < α the positive steady state is unstable andthere is no heteroclinic cycle at infinity. It follows from the Poincar´e-Bendixsontheorem that the ω -limit set of a solution which starts near the steady statebut 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 3 of [3].Thus the second statement of Theorem 1 holds. If α > α then there is againno heteroclinic cycle at infinity. The α -limit set of the solution on the centremanifold of P must then, by the Poincar´e-Bendixson theorem, be either aperiodic solution or the positive steady state. Moreover, if a periodic solutionexists then only the first possibility can occur. Since, however, it follows from[3] that any periodic solution which exists is stable the first possibility is ruledout. There can be no periodic solution and the third case of Theorem 3 of [3]must be realised. This completes the proof of the third statement. Finally, thefourth statement will be proved by contradiction. Let β i be a sequence tendingto α from below. For a given i the system with parameter β i has a uniqueperiodic 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 convergentsubsequence. 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 throughthe point (1 , z ∗ ), which is a periodic solution of the system with parameter value α . This contradicts the fact that there are no such solutions. (cid:4) In the system (15)-(16) the x -axis is an invariant manifold of the flow and thevector field is directed toward positive values of x on the y -axis. For each fixed9hoice of the parameters with ν < (cid:16) βν − ν , (cid:17) . For ν ≥ J = − (1 − ν ) βν − γ (cid:16) − ν βν (cid:17) α (1 − ν ) βν α (cid:16) γ (cid:16) − ν βν (cid:17) − (cid:17) . (24)The determinant of J is α (1 − ν ) βν which is always positive. Thus the stability ofthe steady state is determined by the trace of J , which is α (cid:20) γ (cid:18) − ν βν (cid:19) − (cid:21) − (1 − ν ) βν . (25)If γ ≤ βν − ν then the trace of J is negative for all values of α and the steadystate is always stable. If γ > βν − ν define α = (1 − ν ) γ (1 − ν ) − (1+ βν ) . Then for α < α the trace of J is negative and the steady state is asymptotically stable while for α > α the trace of J is positive and the steady state is a source. For α = α there is a pair of imaginary eigenvalues. If we consider the real part of theeigenvalues as a function of α then it passes through zero when α = α and itsderivative with respect to α at that point is non-zero. Thus a Hopf bifurcationoccurs.In the limiting case ν = 0 it was shown in [3] that the Hopf bifurcation issupercritical so that there exists a stable periodic solution for any α slightlygreater than α . The computation of the Lyapunov number required to obtainthis conclusion becomes considerably more complicated for ν >
0. Rather thantrying to do this in general we will confine ourselves to obtaining some infor-mation for restricted sets of parameters. The Lyapunov number of the Hopfbifurcation is a function of the parameters α , β , γ and ν and we are interestedin its sign. A general formula for this quantity is given in Section 4.4 of [18]. Itis of the form − π b ∆ / 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 inthat case the Hopf bifurcation is supercritical. For ν small and positive f is stillnegative and the bifurcation supercritical. It will now be proved that there alsoexist parameters for which f is positive, so that there exists a subcritical Hopfbifurcation. In that case there exist unstable periodic solutions for α slightlyless than α . Note for comparison that it was shown in [3] that for ν = 0 un-stable periodic solutions never exist. It suffices to treat the case β = 0 sincean example with β small and positive follows by continuity. Since we are onlylooking for some example we can also restrict to the case γ = 2.10ith a suitable normalization the function f is of the following form. α (1 − ν ) ( − a + 2 αa a )+2(1 − ν )( α a − αa ( a + a ))+ α (1 − ν ) ( a a − αa ) + 2 α (1 − ν ) ( α a − a a )+4(1 − ν )( − a + α a a ) + 4(2 αa − α a a ) (26)+(2 α (1 − ν ) + 2(1 − ν ) )( − α a a + a a ) + (1 − ν ) [2 α − (1 − ν )] × [3( − α (1 − ν ) a + 2 a ) + 2(1 − ν )( − a + αa ) + α ((1 − ν ) a − a )]Here the notation a ij is taken from Perko [18]. In order that there exist abifurcation a restriction on ν must be satisfied and in the case γ = 2 it isgiven by ν < . Consider now the limit ν → . Since α = (1 − ν ) − ν at thebifurcation point it tends to infinity in this limit. The highest power of α in theabove expression is α and two terms containing α cancel. Substituting in theexpression for the bifurcation point gives a function depending on ν alone andwe want to examine its behaviour near ν = . To do this it suffices to retainonly those terms in the above expression which contain a power of α which isat least two. It is also the case that the expressions for a and a contain afactor of 1 − ν . Thus to order (1 − ν ) − we get the expression14 α [ − αa ) a − αa ) + 8 a a + 3 a − a ] + . . . (27)The expression in square brackets tends to a positive value as ν → . Thus theleading term in the expression for the Lyapunov number is positive for ν closeto its limiting value. This proves the desired statement. In [3] it was investigated using the Poincar´e compactification in which wayssolutions of (17)-(18) can tend to infinity for large times. Here we want to carryout corresponding calculations for (15)-(16). A useful preliminary step is tointroduce a new time coordinate T satisfying dτdT = 1 + νy γ ( β + x ). Then weget the system dxdT = 1 − xy γ + νy γ ( β + x ) , (28) dydT = α [ xy γ − y − νy γ +1 ( β + x )] . (29)This makes the right hand side into a polynomial while leaving the phase portraitunchanged.The phase portrait is more complicated than that in the case of mass actionkinetics. A schematic picture of it is given in Fig. 1 and its properties aresummarized in the following lemma which is the analogue of Lemma 2 in [3].11 γ γ γ P P P P P P P P Figure 1: Poincar´e compactification.
Lemma 5
Suppose that ν <
1. There is a smooth mapping of the closure of thepositive quadrant into itself mapping the axes into themselves with the followingproperties. The restriction of φ to the open quadrant is a diffeomorphism ontoits image. This image is a region whose closure is a compact set boundedby intervals [0 , x ] and [0 , y ] on the x - and y -axes and four smooth curves γ i , ≤ i ≤
4. The curve γ joins the point P = (0 , y ) with a point P in thepositive quadrant. γ joins the point P with the point P . For 3 ≤ i ≤ γ i joins the point P i − with the point P i and P = ( x , φ in such a way that P and P i , ≤ i ≤
4, are steady statesand the γ i and the image of the x -axis under φ are invariant manifolds. Thereare further steady states P and P i +1 , ≤ i ≤
3, on the boundary belonging tothe interior of γ and γ i +1 , ≤ i ≤
3, respectively.To analyse the case where x becomes large (Case 1 in the terminology of [3])introduce the variables Y = yx , Z = x . Define a new time variable t satisfying dtdT = Z − γ − . The result of the transformation is dYdt = αY γ Z + Y γ +1 Z − αY Z γ +1 − Y Z γ +2 − νY γ +1 ( α + Z )(1 + βZ ) , (30) dZdt = Y γ Z − Z γ +3 − νY γ Z (1 + βZ ) . (31)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 doa quasihomogeneous directional blow-up. An appropriate transformation canbe obtained by using a Newton polygon as in [4]. The coefficients are α = γ and β = γ −
1. (These are the same values as occurred in the blow-up of thecorresponding point for the model (17)-(18).) Thus we use variables ¯ y and ¯ z satisfying ( Y, Z ) = (¯ y γ , ¯ y γ − ¯ z ). In addition we introduce a new time coordinate12 satisfying dsdt = γ − ¯ y γ − . The system becomes d ¯ yds = α ¯ y ¯ z + ¯ y γ +1 ¯ z − α ¯ y ¯ z γ +1 − ¯ y γ ¯ z γ +2 − ν ( α + ¯ y γ − ¯ z )¯ y (1 + β ¯ y γ − ¯ z ) , (32) d ¯ zds = − α ( γ − z + ¯ y γ ¯ z + α ( γ − z γ +2 − ¯ y γ − ¯ z γ +3 + ν [( γ − α − ¯ y γ − ¯ z ]¯ y ¯ z (1 + β ¯ y γ − ¯ z ) . (33)Both axes are invariant under the flow. There is a steady state at the originand one at the point (0 , P . The linearization at theorigin is identically zero.Next the centre manifold of P will be studied. Introducing w = ¯ z − w = ρ ¯ y with ρ = − αν α for γ = 2 and ρ = − νγ for γ >
2. Consider now the case γ = 2, where¯ y ′ = α ¯ y ¯ z − α ¯ y ¯ z γ +1 − αν ¯ y − ¯ y γ + O (¯ y γ +1 ) . (34)Substituting the asymptotic expansion for ¯ z in terms of ¯ y which holds on thecentre manifold into this relation gives¯ y ′ = αρ ¯ y − αρ ¯ y − αν ¯ y − ¯ y + O (¯ y ) = − y + O (¯ y ) . (35) Lemma 6
For γ > y ′ = − γγ − ¯ y γ + o (¯ y γ ) holds on the centremanifold of P . Proof
We use the relation¯ z ′ = − ¯ y γ − + ( γ − − α ¯ z + α ¯ z γ +2 + αν ¯ y ¯ z ] + O (¯ y γ ) . (36)Substituting this into the evolution equation for ¯ y gives¯ y ′ = − γγ − y γ − γ − y ¯ z − ¯ z ′ + O (¯ y γ +1 ) . (37)With this it is possible to adapt the argument used to analyse the centre man-ifold of P in [3] to get the desired conclusion as follows. Suppose that weknow that ¯ y ′ = O (¯ y k ) for some k with 2 ≤ k ≤ γ −
1. Then it follows that¯ z ′ = O (¯ y k +1 ). Hence ¯ y ′ = O (¯ y k +1 ). After finitely many steps we get the secondconclusion of the lemma. (cid:4) We see that the flow on the centre manifold is towards P and since the non-zero eigenvalue of the linearization at that point is positive P is a topologicalsaddle. In fact the flow on the boundary is everywhere away from P . We nextblow up the origin in the coordinates (¯ y, ¯ z ). This time the procedure described in[4] leads to the choice of coefficients α = β = 1. Blow-ups in the two coordinatedirections are required. The only terms in the resulting equations which willbe written explicitly are those which have a direct influence on the analysis13hich follows. In the first transformed system, with ¯ y = ˜ y and ¯ z = ˜ y ˜ z , theequations are ˜ y ′ = ˜ y [ − α ( ν − ˜ z )˜ y + · · · ] , (38)˜ z ′ = ˜ y [ γα ( ν − ˜ z )˜ z + · · · ] . (39)A change of time coordinate eliminates the common factor ˜ y . On the boundarythere is a steady state at the point (0 , ν ), which corresponds to P . The origin ofthis coordinate system corresponds to P . The terms which have been retainedsuffice to determine the steady states on the boundary and the linearization ofthe system at those points. The point P also appears in the second transformedsystem but since it can be analysed in the chart corresponding to the firsttransformed system the second transformed system, with ¯ y = ˜ y ˜ z and ¯ z = ˜ z ,is only needed to analyse the steady state P at the origin of that coordinatesystem. For this purpose the only terms which need to be retained are˜ y ′ = ˜ z [ γα ˜ y + · · · ] , (40)˜ z ′ = ˜ z [ − α ( γ − z + · · · ] . (41)The common factor ˜ z can be eliminated by a change of time coordinate. Theorigin of this coordinate system corresponds to P . 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 will be studied in the case γ = 2. We do notexpect that the case γ > γ = 2 has been worked out. Thecentre subspace is parallel to the ˜ y -axis. We have ˜ z ′ = O (˜ y ) on the centremanifold and this implies that if ˜ z = ν + w then2 ανw = [ ν (1 − ν ) + 2 αν + 2 αβν ]˜ y + . . . (42)It follows that provided ν < P is awayfrom P . 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 quasihomogeneousdirectional blow-up of the origin in the ( Y, Z ) coordinate system. In this case(
Y, Z ) = (¯ y ¯ z γ , ¯ z γ − ). The time coordinate is transformed using the relation dsdt = γ − ¯ z γ − . The resulting system is¯ y ′ = ( γ − α ¯ y γ − α ¯ y − αν ¯ y γ +1 ¯ z (1 + β ¯ z γ − )+¯ y γ +1 ¯ z γ − ¯ y ¯ z γ − − ν ¯ y γ +1 ¯ z γ (1 + β ¯ z γ − )] − γ [¯ y γ +1 ¯ z γ − ¯ y ¯ z γ − − ν ¯ y γ +1 ¯ z γ (1 + β ¯ z γ − )] , (43)¯ z ′ = ¯ y γ ¯ z γ +1 − ¯ z γ − ν ¯ y γ ¯ z γ +1 (1 + β ¯ z γ − ) . (44)There is a steady state at the point (1 ,
0) but since it is just another representa-tion of P it does not need to be analysed further. The origin of this coordinatesystem corresponds to P . The ¯ z -axis is a centre manifold for P and the flowthere is towards P . Since the non-zero eigenvalue of the linearization at P isnegative it can be concluded that P is a sink.14aving completed the analysis of the case where x gets large we now turn tothe case where where y gets large (Case 2 in the terminology of [3]), with newvariables X = xy and Z = y . The result is dXdT = 1 Z γ +1 [ Z γ +2 − XZ + νZ ( X + βZ ) − αX Z + αXZ γ +1 + ανX ( X + βZ )] , (45) dZdT = 1 Z γ +1 [ − αXZ + αZ γ +2 + ανZ ( X + βZ )] . (46)The common factor Z γ +1 can be removed by a suitable change of time coordinatesatisfying dtdT = Z − γ − . The linearization of the resulting system about theorigin is identically zero so that it is again necessary to do a blow-up. In thiscase a calculation using a Newton polygon gives the exponents α = 1 and β = 1.The transformation in the X direction uses the relation ( X, Z ) = (¯ x , ¯ x ¯ z ). Theresulting system is d ¯ x dt = ¯ x [¯ x γ +11 ¯ z γ +21 − ¯ x ¯ z + ν ¯ x ¯ z (1 + β ¯ z ) − α ¯ x ¯ z + α ¯ x γ +11 ¯ z γ +11 + αν ¯ x (1 + β ¯ z )] , (47) d ¯ z dt = ¯ x [ − ¯ x γ ¯ z γ +31 + ¯ z − ν ¯ z (1 + β ¯ z )] . (48)The origin of this coordinate system corresponds to P . By a change of timecoordinate we can remove the factor ¯ x . The linearization of the system whichresults at the origin has one positive eigenvalue and the ¯ z -axis is invariant anddefines a centre manifold at that point. It can be concluded that P is a source.If ν < (cid:16) , − νβν (cid:17) which corresponds to thepoint P . That point is a hyperbolic saddle whose stable manifold is the ¯ z -axis.The transformation in the Z direction uses the relation ( X, Z ) = (¯ x ¯ z , ¯ z ).The resulting system is d ¯ x dt = ¯ z [¯ z γ − ¯ x + ν ( β + ¯ x )] , (49) d ¯ z dt = − αXZ + αZ γ +2 + ανZ ( X + βZ )= ¯ z [ − α ¯ x ¯ z + α ¯ z γ +12 + αν ¯ z ( β + ¯ x )] . (50)The origin of this coordinate system is P . By a change of time coordinate wecan remove the factor ¯ z . In the system which results there is inflow on the ¯ z -axis while the ¯ x -axis is invariant and corresponds to the ¯ z -axis in the previoussystem. Note that the point P is not a steady state.The facts which have now been collected imply strong restrictions on thepossible ω -limit sets of solutions. The only points of the boundary which theycan contain are those on the part connecting P and P . Poincar´e-Bendixsontheory implies that the ω -limit set of a positive solution must be either a point15which can only be the positive steady state, P or P ), a periodic solutionor a heteroclinic cycle joining P and P . The last of these can only occur ifthe centre manifolds of P and P coincide. Note that any periodic solution orheteroclinic cycle must contain the positive steady state in its interior.We have the following analogue of Theorem 1 of [3]. Theorem 2
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)) , (51) y ( τ ) = y e − ατ (1 + o (1)) . (52)There exists a solution, unique up to time translation, which has the asymptoticbehaviour x ( τ ) = τ (1 + o (1)) , (53) y ( τ ) = τ − γ − (1 + o (1)) . (54) Proof
This theorem can be proved in the same way as Theorem 1 of [3]. Theonly extra element is that it is necessary to use the fact that for this type ofsolution the time coordinates τ and T are asymptotically equal. The parameter ν does not contribute to the leading order asymptotics. (cid:4) To understand the global phase portrait it is helpful to understand the geometryof the nullclines N and N given by ˙ x = 0 and ˙ y = 0, respectively. We restrictconsideration to the case ν < N and N intersect in a single point.The equation for N can be expressed in the equivalent forms y = (cid:20) − βν + (1 − ν ) x (cid:21) γ , x = 11 − ν ( y − γ + βν ) . (55)Thus on N the coordinate y can be expressed as a smooth function of x witha smooth inverse. Note, however, that while the second function is defined forall positive y the first is only defined for x > βν − ν . The equation for N can beexpressed in the form x = y − γ + βνy − νy . (56)Thus on N the coordinate x can be expressed as a (locally defined) smoothfunction of y . Since x can be written as a function of y in both cases and thereis only one point of intersection it is clear that the complement of N ∪ N isa union of the four connected components defined by the signs of ˙ x and ˙ y . Aschematic picture of the null clines is given in Fig. 2. As in [3] we denote theregions with the sign combinations (+ , − ), (+ , +), ( − , +) and ( − , − ) by U , U ,16 ν βν − ν N N U U U U Figure 2: Nullclines. U and U , respectively. Where ˙ y = 0 and ˙ x = 0 we can use the fact that N is a graph over the y -axis to conclude that a solution can only pass from U to U and U to U and not the other way round. Similarly the fact that N is agraph over the y -axis implies that a solution can only pass from U to U and U to U and not the other way round. Thus the possible passages between theregions U i are just as in the case with mass action kinetics. Let L be the part ofthe horizontal line segment joining the positive steady state to the y -axis withthe endpoint on the axis excluded. The part of L excluding the steady state iscontained in U . Lemma 7
In the Michaelis-Menten system each of the centre manifolds of P and P contains a point of L in its closure. Proof
A point on the centre manifold of P which is sufficiently close to P liesin the region U . If we follow a solution which lies on this manifold forwards intime then it must either tend to the positive steady state as t → ∞ or it mustenter U after a finite time and in the latter case it must enter U . Once it hasdone so it must either tend to the positive steady state as t → ∞ or it mustmeet L after a finite time. Similarly a solution on the centre manifold of P which starts close to P must, when followed backwards in time, either convergeto the positive steady state as t → −∞ or meet L after a finite time. (cid:4) Denote the x -coordinates of the points of L in the closure of the centremanifolds of P and P for given value of α and ν by ξ ( α, ν ) and ξ ( α, ν ). Lemma 8
The function ξ − ξ describing the separation of the points wherethe centre manifolds of P and P reach y = 1 is continuous. Proof
The proof is similar to that of Lemma 2. The essential facts which mustbe used are that any solution which approaches P close enough in the past timedirection and does not lie on the centre manifold of P cannot remain close to P forever and the corresponding statement with ’past’ replaced by ’future’ and P by P . (cid:4) Lemma 9
There are pairs of parameters ( α, ν ) for which the function ξ − ξ describing the separation of the points where the centre manifolds of P and P y = 1 is negative. For ν > < α ≤ α either there exists anunstable periodic solution or ξ − ξ is positive. If there is some α ≤ α for whichno periodic solutions exist then there exists an α with ξ ( α , ν ) = ξ ( α , ν ). Proof
For α sufficiently large the positive steady state is a source and thus ξ ( α, ν ) <
1. Thus in order to prove the first part of the lemma it sufficesto show that for some ( α, ν ) we have ξ ( α, ν ) = 1. To prove this we proceedas in the case of mass action kinetics. First the system is transformed to thecoordinates (¯ y, ¯ z ) and then the quantities ǫ and ¯ w are introduced. The righthand side of each equation in the Michaelis-Menten case is the sum of theright hand side of the corresponding equation in the mass action case and anexpression which can be written as ǫ − ν times a function which is regular inthe limit ǫ →
0. Fixing ν and letting ǫ tend to zero would cause this term toexplode. Instead we let ǫ and ν tend to zero in such a way that ν = ǫ . Thenthe second summand behaves in a smooth manner as ǫ → (cid:4) Note that for the choice of parameters in the first part of Lemma 9 there is aheteroclinic orbit joining the positive steady state to the point P . It follows thatfor these values of the parameters no periodic solutions exist. This is because aperiodic solution would have to contain the positive steady state in its interiorand therefore would have to cross the heteroclinic orbit.There is no straightforward generalization of the monotonicity result ofLemma 4 to the Michaelis-Menten case. The proof of monotonicity fails forthe centre manifold of P since it may pass through the region U . For thisreason even in a case where the existence of a zero of ξ − ξ can be proved wedo not get its uniqueness, Moreover, we do not get the analogue of the stabilitystatement in the mass action case. It is possible to do a calculation analogousto that done to determine the stability of the heteroclinic cycle in [3]. Unfortu-nately in the estimate for the return map the power γ is replaced by the powerone and this gives no information about stability. The following theorem sumsup the results obtained. Theorem 3
The Michaelis-Menten system (15)-(16) has the following proper-ties.1. For each choice of the parameters ( α, ν ) with ν < (cid:16) βν − ν , (cid:17) .2. If γ ≤ βν − ν the steady state is stable. Otherwise if α < α = (1 − ν ) γ (1 − ν ) − (1+ βν ) it is stable and for α > α unstable.3. For α = α a generic Hopf bifurcation occurs. Parameters can be chosen soas to make it supercritical or subcritical.4. For given ( α, ν ) there exist positive numbers x and y such that if a solutionsatisfies x ( t ) ≥ x and y ( t ) ≤ y at some time t then it has the late timeasymptotics described in Theorem 2. 18. For γ = 2 there exists a choice of positive parameters α and ν for which allsolutions other than the steady state have the late time asymptotics describedin Theorem 2.6. If for γ = 2 and given ν there exist no periodic solutions for α sufficientlysmall then there exits a heteroclinic cycle passing through the steady states P and P in that order.It should be noted that it has not been shown here whether the case describedin point 6. ever occurs. It has been shown that the basic Selkov system admits solutions with unboundedoscillations and that the diameter of the image of a periodic solution can tendto infinity as α approaches a finite limit, thus completing the results of [3] onthat system and rigorously confirming a claim made in [19]. Note that somestatements related to this issue have been made in [8] but that that referencedoes not contain rigorous proofs of those statements. One remaining questionis that of the rate with which the diameter of the image of the periodic solutiontends to infinity as the critical parameter value α is approached. A suggestionfor this has been made in [16] for the case γ = 2 but there is neither a rigorousproof that this suggestion is correct nor a generalization of the statement tohigher values of γ .It was also investigated which properties of the basic Selkov system persistin 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 allparameter values. The question of whether the five-dimensional system fromwhich the Michaelis-Menten system itself was derived has unbounded solutionsremains open. It was shown that for suitable parameter values unstable periodicsolutions of the Michaelis-Menten system exist. It was left open whether thereexist unbounded oscillatory solutions or periodic solutions whose images havearbitrarily large diameter for bounded ranges of the parameters.The unbounded solutions cast doubt on the suitability of the Selkov modelfor describing glycolytic oscillations. An alternative model often preferred tothe Selkov model is that of Goldbeter and Lefever [9]. There the amplitudeof the periodic solutions created in a Hopf bifurcation increases to a maximumbefore decreasing again to zero at a point where the periodic solutions vanishagain in a second Hopf bifurcation. It has been proved by d’Onofrio [6], on thebasis of an analysis of a more general class of systems in [17], that all solutionsof the Goldbeter-Lefever model are bounded. Further aspects of the dynamicsof solutions of that model have been studied in [6] and a sophisticated analysisof some of its properties has been carried out in [13].The questions of the origin of the unbounded solutions and how they couldbe eliminated by modifying the system have been discussed in [15]. The origin ofthe unbounded growth can be seen in the constant source term in the equation19or x . This corresponds to an unlimited supply of the substrate. In [15] this iscalled the pooled chemical approximation. If this is replaced by a mechanismwhere the substrate is formed from a precursor which itself is limited in quantitythen the oscillations only grow within a finite time period before decaying again.An alternative modification is to introduce an additional uncatalysed conversionof the substrate into the product. The resulting system is called the (cubic)autocatalator. According to the analysis of [15] this leads to a situation similarto that described above for the Goldbeter-Lefever model and the unboundedoscillations are absent. Some aspects of this type of model have been analysedrigorously in [10].The Selkov system and other related ones are model cases for understandingoscillations in biological and chemical systems. The study of equations of thistype raises a number of issues. In what ways can heuristic and numerical resultsbe made into rigorous theorems? How can we understand the relations betweenthe choices made in modelling and the relevance of the resulting models for theapplications? The present paper is intended as a contribution to the clarificationof these issues. Acknowledgements
One of the authors (ADR) is grateful to James Sneyd andJ´anos T´oth for helpful discussions.