In distributive phosphorylation catalytic constants enable non-trivial dynamics

Ordered distributive double phosphorylation is a recurrent motif in intracellular signaling and control. It is either sequential (where the site phosphorylated last is dephosphorylated first) or cyclic (where the site phosphorylated first is dephosphorylated first). Sequential distributive double phosphorylation has been extensively studied and an inequality involving only the catalytic constants of kinase and phosphatase is known to be sufficient for multistationarity. As multistationarity is necessary for bistability it has been argued that these constants enable bistability. Here we show for cyclic distributive double phosphorylation that if its catalytic constants satisfy an analogous inequality, then Hopf bifurcations and hence sustained oscillations can occur. Hence we argue that in distributive double phosphorylation (sequential or distributive) the catalytic constants enable non-trivial dynamics. In fact, if the rate constant values in a network of cyclic distributive double phosphorylation satisfy this inequality, then a network of sequential distributive double phosphorylation with the same rate constant values will show multistationarity—albeit for different values of the total concentrations. For cyclic distributive double phosphorylation we further describe a procedure to generate rate constant values where Hopf bifurcations and hence sustained oscillations can occur. This may, for example, allow for an efficient sampling of oscillatory regions in parameter space. Our analysis is greatly simplified by the fact that it is possible to reduce the network of cyclic distributive double phosphorylation to what we call a network with a single extreme ray. We summarize key properties of these networks. Supplementary Information The online version contains supplementary material available at 10.1007/s00285-024-02114-8.


Introduction
Phosphorylation is a process where proteins are altered by adding and removing phosphate groups at designated binding sites.It is a recurrent motif in many large reaction networks involved in intracellular signaling and control [25].Often spatial effects are neglected and the time dynamics of the participating chemical species concentrations is described by ordinary differential equations.There exists a plethora of small ODE models that include phosphorylation at one or two binding sites together with the interaction of various regulating chemical species, see, for example, [21].These models exhibit a wide range of dynamical properties ranging from multistationarity and bistability to sustained oscillations [21].
Phosphorylation and dephosphorylation are catalyzed by two enzymes, a kinase and a phosphatase.As described in [4,23], this process can either be processive or distributive: if it is processive, then all available binding sites are phosphorylated or de-phosphorylated upon binding of protein and kinase/phosphatase.If the process is distributive, then at most one binding site is modified upon each binding of protein and kinase or phosphatase.Furthermore, multisite phosphorylation and dephosphorylation can occur at a random sequence of binding sites or at an ordered sequence of binding sites.An ordered mechanism is either sequential or cyclic [4,23].In a sequential mechanism the last site to be phosphorylated is dephosphorylated first, while in a cyclic mechanism the first site to be phosphorylated is also dephosphorylated first -as depicted in the reaction schemes of Fig. 1.In panel (a) the subscript ij denotes the state of the phosphorylation sites: 0 unphosphorylated, 1 phosphorylated (e.g. S 10 denotes those molecules of S, where the first site is phosphorylated and the second site is unphosphorylated).In panel (b) the subscript i denotes the number of attached phosphate groups (e.g. S 0 denotes unphosphorylated protein).
Networks of sequential phosphorylation have been studied extensively and it has been shown that already networks without any form of regulation can exhibit non-trivial dynamics.It is, for example, known that models of sequential and processive phosphorylation and dephosphorylation have a unique, globally attracting steady state [12].For sequential distributive phosphorylation and dephosphorylation multistationarity and bistability have been established [8,15,16].And for a mixed mechanism of sequential distributive phosphorylation and processive dephosphorylation at two binding sites Hopf bifurcations and sustained oscillations have been described in [24] and [10].For more information about the dynamics of multisite phosphorylation systems see the review [13] and the references therein.
Networks of cyclic phosphorylation have been studied in [26].The authors study in particular the following mass action network (1) derived from the reaction scheme of Fig. 1a (where the notation is the same as in Fig. 1a and KS 00 , KS 10 , F S 11 and F S 01 denote enzyme-substrate complexes): In [26] it is shown that for every admissible positive value of the total concentrations and all positive values of the reaction rate constants there exists a unique positive steady state.
The authors furthermore provide parameter values where a Hopf bifurcation occurs (with the total amount of kinase as a bifurcation parameter) and sustained oscillations emerge [26,Fig. 4].The authors however neither describe how the respective parameter values were found nor do they explain how parameter values leading to oscillations can be found.
Here we study the same network (1) and provide an answer to the latter question.In particular, we show that if the rate constants κ 3 , κ 6 , κ 9 and κ 12 satisfy the inequality then Hopf-bifurcations and sustained oscillations can occur.In enzyme catalysis these constants are known as catalytic constants (cf. the discussion in Section 4 or [9]).
To be more precise, we analyze the following irreversible subnetwork of (1) and show that if the values of the catalytic constants κ 3 , κ 6 , κ 9 and κ 12 satisfy (2), then there exist values for κ 1 , κ 4 , κ 7 , κ 10 and the total concentrations such that the steady state Jacobian of network (3) has a complex-conjugate pair of eigenvalues on the imaginary axis.Furthermore, if this eigenvalue pair crosses the imaginary axis as one of the parameters is varied, then a Hopf bifurcation occurs at this particular steady state (Theorem 3.4).We describe a derivative condition for this crossing and a procedure to find such parameter values.Based on a result by [2], we then argue that if a supercritical Hopf bifurcation occurs in network (3) and a stable limit cycle emerges, then the full network (1) will have a stable limit cycle as well (for appropriately chosen values of κ 2 , κ 5 , κ 8 and κ 11 ).This is Theorem 3.9.Finally, we compare our results for cyclic distributive double phosphorylation with those obtained in [9] for sequential distributive double phosphorylation.In [9] an inequality that is sufficient for multistationarity has been described.Remarkably this inequality involves the catalytic constants in the same way as our inequality (2).As multistationarity is necessary for bistability it is argued in [9] that these rate constants enable bistability in sequential and distributive double phosphorylation.Hence we conclude that in distributive phosphorylation (sequential or cyclic) these constants enable non-trivial dynamics.
The paper is organized as follows: to arrive at our results the ODEs derived from network (1) and (3) are analyzed.For this purpose we recall in Section 2 some well known facts about ODEs defined by reaction networks and introduce a special class of reaction networks that we call networks with a single extreme ray.In this section we also derive conditions for x 8 x 9 x 10 K F S 00 S 10 S 01 S 11 KS 00 KS 10 F S 01 F S 11 Table 1: Species and concentration variables for network (1) and (3). a simple Hopf bifurcation in such networks.In Section 3 we first analyze network (3) and verify these bifurcation conditions in Theorem 3.4.Then we turn to the full network (1) and Theorem 3.9.We also present the procedure to determine rate constant and total concentration values.In Section 4 we discuss inequality (2) in the light of the the results presented in [9].Appendix A and B contain some of the longer proofs of the results presented in Section 2. Appendix C, D and E contain information to reproduce the numerical results displayed in the figures throughout the paper.

Biochemical reaction networks with mass action kinetics
To establish our results we exploit the special structure of the Jacobian of a certain class of reaction networks that we call networks with a single extreme ray.We introduce this class here in full generality.To this end we first consider a general reaction network with n species and r reactions in Section 2.1 and recall the structure of the ODEs defined by such a general network.In Section 2.2 we discuss steady states and formally define networks with a single extreme ray.In Section 2.3 we present a formula for the Jacobian of a general reaction network at steady state.In Section 2.4 we discuss the steady state Jacobian of networks with a single extreme ray.Finally, in Section 2.5 we present conditions for simple Hopf bifurcations in networks of this kind.

Reaction networks with n species and r reactions
We briefly introduce the relevant notation, for a more detailed discussion we refer to the large body of literature on mass action networks, for example, [7] or [11].
To every chemical species we associate a variable x i denoting its concentration.For network (1) and (3) we use the association as given in Table 1.
Consider network (1), nodes like S 00 + K are called complexes.To every complex we associate a vector y ∈ R n representing the stoichiometry of the associated chemical species.The complex S 00 + K consists of one unit of S 00 and one unit of K, hence, in the ordering of Table 1 the vector y T = (1, 0, 1, 0, 0, 0, 0, 0, 0, 0) represents its stoichiometry.In a similar way one arrives at the ten complex vectors for networks (1) and (3) given in Table 3 of Appendix C.
To every reaction r (l) we associate the difference of the complexes at the tip and the tail of the reaction arrow.For example, to the reaction S 00 + K → S 00 K we associate the vector r (l) = y (2) − y (1) = (−1, −1, 0, 1, 0, 0, 0, 0, 0, 0) T (using the labeling of complexes given in Table 3 of Appendix C).All reaction vectors r (i) are collected as columns of the stoichiometric matrix S ∈ R n×r : S = r (1) . . .r (r) . (4) To every reaction r (l) we further associate a reaction rate function v l (k, x) describing the 'speed' of the reaction.We consider only mass action kinetics, hence the reaction rate function of reaction r (l) : y (i) → y (j) is given by the monomial function where k l is a parameter called the rate constant and x y is the customary shorthand notation for the product n l=1 x y l l of two n-vectors x and y.We collect the vectors at the tail of every reaction as columns of a matrix Y and note that this matrix may contain several copies of the same vector: Y = y (1) . . .y (r)  ( We collect all rate constants in a vector k T = (k 1 , . . ., k r ) and the monomials x y (i) in a vector φ(x) = (x y (1) , . . ., x y (r) ) T , where the vectors y (i) reference the columns of the matrix Y defined in (5).The reaction rate function is then defined using k and φ(x) as After an ordering of species and reactions is fixed, every reaction network with mass action kinetics defines a matrix S and a reaction rate function v(k, x) in a unique way.These objects in turn define the following system of ODEs: where ẋ denotes the vector of derivatives with respect to time.Example are abundant in the literature, see, for example, [7, Section 2] or [11].
Often the matrix S does not have full row rank.In this case let W be a matrix whose columns span ker(S T ), that is a full rank matrix W with W T S ≡ 0. Then W T ẋ ≡ 0 and for every solution x(t) with initial value x(0) = x 0 one obtains Hence, if rank(S) = s < n, then one obtains n − s conservation relations

Steady states
We are interested in points (k,x) that are solutions of If (k,x) is a solution of (10), then x is a steady state of (8) for the rate constants k.We proceed as in [6] and express the reaction rates v(k, x) at a steady state as a nonnegative combination of the extreme vectors of the pointed polyhedral cone ker(S) ∩ R r ≥0 .This idea goes back to Clarke and coworkers (cf.for example, [5,13]) and it is as follows: (k,x) satisfy (10), if and only if the corresponding v(k, x) is such that v(k, x) ∈ ker(S) ∩ R r ≥0 (cf.eg.[6]).
Convex polyhedral cones have a finite number of extreme vectors (up to a scalar positive multiplication [22]).Therefore, any element v of such a cone can be represented as a nonnegative linear combination of its extreme vectors {E 1 , . . ., where E is the matrix with columns E 1 , . . ., E l and λ T = (λ 1 , . . ., λ l ).
Remark 2.1 (The relative interior of ker(S) ∩ R r ≥0 ).(A) As explained in, for example, [27], a system (8) has positive solutions (k,x), if and only if the matrix E does not contain a zero row.Hence we will only consider networks where the matrix E are of this kind.
(B) We are only interested in positive values of k and x.Thus we are only interested in v(k, x) that are strictly positive and hence belong to the relative interior of the cone ker(S) ∩ R r ≥0 .Consequently we are only interested in those λ ∈ R l ≥0 that yield a positive Eλ.In the remainder of the paper we therefore restrict λ ∈ R l ≥0 to the set: This set has been introduced in [7], cf.[7,Remark 4].
(C) A more detailed discussion of the relation between positive solutions of ( 10), the cone ker(S) ∩ R r ≥0 and the generator matrix E can be found in [27], cf. in particular [27, Proposition 6].
In the analysis of network (1) we will later analyze subnetworks that fit the following definition: Definition 1 (Networks with a single extreme ray).Consider a reaction network with stoichiometric matrix S. If the cone ker(S) ∩ R r ≥0 is generated by a single positive vector then we call it a network with a single extreme ray.
Consequently, if (k,x) satisfy the steady state equation (10), then v(k, x) ∈ ker(S) ∩ R r ≥0 and there exists a vector λ ∈ Λ (E), such that (11) is satisfied.Thus one may parameterize all reaction rates at steady states via the equation Likewise, given some λ ∈ Λ (E) and a positive x one obtains a positive vector k such that (k,x) satisfy the steady state equation ( 10) by the following formula: where We observe that (i) the formula ( 15) is well defined as by assumption x > 0, (ii) that the vector k from ( 14) is positive since all entries of Eλ are positive (as, by assumption, the matrix E does not have any zero rows, cf.Remark 2.1) and (iii) that the formula ( 14) is obtained by solving the k-linear equation ( 13) for k.

The Jacobian at steady state
The equations ( 11) & ( 13) introduce a parametrization of the reaction rates v(k, x) at steady states in terms of the convex parameters h and λ.This parameterization can be used to parameterize the Jacobian at steady state: as explained in detail in [6] or [5], if (k,x) satisfy the steady state equation ( 10) and hence v(k, x) is such that (13) holds, then the Jacobian evaluated at that (k, x) is given by the following formula (cf.[6,Proposition 2]): where Y denotes the matrix introduced in (5).

The Jacobian of networks with a single extreme vector
For any reaction network where the cone ker(S)∩R r ≥0 is spanned by a single positive vector E the parameter λ is a positive scalar and formula ( 16) is equivalent to We use J λ (h) to denote the Jacobian in this special case.
Here we comment on the characteristic polynomial of general matrices J λ (h), that is on the polynomial det(µI − J λ (h)) of an n × n Jacobian of the form (17).We will assume that rank(J λ ) = s ≤ n and we will use the symbol a i (h, λ) to denote the coefficients of its characteristic polynomial.We will also discuss the case λ = 1, that is, the characteristic polynomial det(µI − J 1 (h)) of J 1 (h) with coefficients b i (h).For this purpose, the characteristic polynomial of the Jacobian J λ (h) is denoted by while the characteristic polynomial of J 1 (h) is denoted by Concerning these polynomials we have the following corollary of Lemma A.4 in Appendix A where we show the relationship between the coefficients a i (h, λ) and b i (h).
Corollary 2.2.Suppose the matrix E consists of a single positive column vector.Let J λ (h) be as in (17) with rank(J λ (h)) = rank(J 1 (h)) = s < n and let a i (λ, h), b i (h) be the coefficients of the characteristic polynomials in (18) and in (19).Then the coefficients a i (λ, h) and b i (h) satisfy: Moreover, the polynomial det(µI − J λ (h)) is given by the following formula: We observe the following consequences of Corollary 2.2: (II) In our setting λ > 0, hence the sign of the real part of the eigenvalues Re (III) In particular, the matrix J λ (h) has a purely imaginary pair of eigenvalues ±iλω(h) if and only if J 1 (h) has a purely imaginary pair ±iω(h).

Simple Hopf bifurcations in networks with a single extreme ray
We recall the definition of a simple Hopf bifurcation for a parameter dependent system of ODEs of the form ẋ = g p (x), where x ∈ R s , and g p (x) varies smoothly in p and x.Let x * ∈ R s be a steady state of the ODE system for some fixed value p 0 , that is, g p 0 (x * ) = 0. Furthermore, we assume that we have a smooth curve of steady states around p 0 : That is, g p (x(p)) = 0 for all p close enough to p 0 .Further let J(x(p), p) be the Jacobian of g p (x) evaluated at x(p).If, as p varies, a complex-conjugate pair of eigenvalues of J(x(p), p) crosses the imaginary axis, then there exists a Hopf bifurcation at (x(p 0 , p 0 )).A simple Hopf bifurcation occurs at (x(p 0 ), p 0 ), if no other eigenvalue crosses the imaginary axis at the same value p 0 .In this case, a limit cycle arises.If the Hopf bifurcation is supercritical, then stable periodic solutions are generated for nearby parameter values [19].
Similar to [10] and [6] we want to build on results described in [19] and [28] to establish Hopf bifurcations.As in previous work (cf.e.g.[6,19,28]), we will use a criterion based on the following Hurwitz determinants: Definition 2. The i-th Hurwitz matrix of a univariate polynomial p(z) = a 0 z s + a 1 z s−1 + • • • + a s is the following i × i matrix: , in which the (k, l)-th entry is a 2k−l as long as 0 ≤ 2k − l ≤ s, and 0 otherwise.The determinants det(H i ) are called Hurwitz determinants.
To every square matrix one can associate Hurwitz matrices via its characteristic polynomial in an analogous way.In the following we consider two families of Hurwitz matrices constructed according to Definition 2: matrices H l (λ, h) obtained from coefficients a i (λ, h) of ( 18) and matrices G l (h) of coefficients b i (h) of (19).Our analysis is greatly simplified by the relationship between the two families established in Proposition 2.4 below (a proof is given in Appendix B): Proposition 2.4.The Hurwitz determinants for the characteristic polynomials ( 18) and ( 19) satisfy the following equation The following proposition is a specialization of [6, Proposition 1] and [28,Theorem 2] to a network where ker(S)∩R r ≥0 is generated by a single positive vector.It implies that for such networks to detect a simple Hopf bifurcation one only needs to study the polynomial (19) and its Hurwitz matrices G i , i = 1, 2, . . ., s. Proposition 2.5.Let ẋ = Sv(k, x) be an ODE system model for a reaction network with mass action kinetics where rank(S) = s.Suppose that the matrix E in (11) consists of a single positive vector and let J λ (h) be the corresponding Jacobian in convex parameters.Further let the characteristic polynomials of J λ (h) and J 1 (h) be as in (18) and (19), respectively.If there exists a fixed value then (a) J 1 (h * ) has a single pair of purely imaginary eigenvalues, (b) J λ (h * ) has a single pair of purely imaginary eigenvalues for all λ > 0, (c) for the dynamical system ẋ = Sv(k, x) there exists a simple Hopf bifurcation at h = h * for all λ > 0 if there exists some l ∈ {1, 2, . . ., n} such that Proof.(c) The fact that J λ (h * l ) for all λ > 0 has a single pair of purely imaginary eigenvalues has been established in (b).Proposition 2.4 implies the following relationship between the two derivatives of det (H s−1 (h, λ)) and det (G s−1 (h, λ)) with respect to the h l : Thus, if (23) holds for some l ∈ {1, . . ., n}, then, for the same l, one has Therefore, by (a), (b) and by [28, Theorem 2&Remark 2], (c) follows as well.
Remark 2.6.The parameter λ is not suited as a bifurcation parameter in (23) as by Propo- 3 explaining that the existence of a pair of purely imaginary eigenvalues is independent of λ).

Analysis of networks of cyclic distributive double phosphorylation
To derive a dynamical model of network (1) we assign the variables given in Table 1 to the chemical species and obtain the following ODEs: Using the same variables x 1 , . . ., x 10 as in Table 1, we obtain the following ODEs: Both networks have the same set of three conservation relations, one for the total amount of kinase (c 1 ), phosphatase (c 2 ), and substrate (c 3 ), respectively: x 3 + x 4 + x 5 + x 6 + x 7 + x 8 + x 9 + x 10 = c 3

Steady states of network (3)
For the network (3) the matrix E consists of a single vector of all 1: hence network (3) is a network with a single extreme ray and we can apply the results of Section 2.4 and 2.5.Given E from ( 27) condition (13) consequently becomes (cf.eq. ( 49) of Appendix C for v(k, x)): These equations can be solved for x (in terms of k and λ) to obtain the following steady state parameterization: where x 1 and x 2 are arbitrary positive numbers.Similarly, solving for k one obtains where h i = 1 x i , i = 1, . . ., 10 (this is (14) for network (3) with v(k, x) given in Appendix C).

Simple Hopf bifurcations for network (3)
The Jacobian J λ (h) of network (3) computed via ( 16) is For J λ (h) given in (30) one has rank(J λ (h)) = 7. Hence its characteristic polynomial is of the following form (cf. supplementary file Cyc dd 2 coeffs charpoly.nband Corollary 2.2): where the coefficients b 1 (h), . . ., b 7 (h) are the coefficients of the characteristic polynomial of the matrix J 1 (h).With the next corollary we adapt Proposition 2.5 to the Jacobian of network (3).This allows us to work with the simpler coefficients b i (h).Then, for all λ > 0, the dynamical system (25a) -(25j) has a simple Hopf bifurcation at h = h * , if additionally ∃l ∈ {1, . . ., 10} such that The coefficients b 0 (h), . . ., b 5 (h) and b 7 (h), as well as the Hurwitz determinants det(G 2 (h)), . . ., det(G 5 (h)) contain only monomials with positive sign (cf.supplementary file Cyc dd 2 coeffs charpoly This establishes the following Lemma and Corollary: Lemma 3.2.Consider the coefficients of the characteristic polynomial (31) of J λ (h) given in (30) for λ = 1 and its Hurwitz determinants.For all h > 0 the following holds: The Hurwitz determinant det(G 6 (h)) contains monomials of both signs (cf.supplementary file cyc dd 2 detH6.txt), hence it can potentially be zero.To establish this, we consider det(G 6 (h)) as a polynomial in h 1 , h 2 , h 3 and h 6 only and study its Newton polytope.Using polymake [1,14] we compute the following hyperplane representation of this Newton polytope: We study the coefficients of det(G 6 ) as a polynomial in h 1 , h 2 , h 3 and h 6 and find that some of these contain the factor h 10 h 7 − h 8 h 9 .
In fact, visual inspection of all coefficients shows that all monomials with exponent vectors contained in the hyperplane have such a coefficient that factors h 10 h 7 − h 8 h 9 .Hence we want to make those monomials dominant.To achieve this we use the following transformation (based on the normal vector of the hyperplane): The result is the following degree 18 polynomial in t, with coefficients that are polynomials in h 4 , h 5 and h 7 , . . ., h 10 : As the constant coefficient is a sum of positive monomials one has D 6 (0) > 0. Thus, if then there exists a t 1 > 0 such that D 6 (t 1 ) = 0 and D 6 (t) < 0 for t > t 1 by the Intermediate Value Theorem.These observations are the basis for the following result: The existence of t 1 with D 6 (t) < 0, t > t 1 and D 6 (t 1 ) = 0 has been established above.By Lemma 3.2 b 1 , . . .b 5 and b 7 as well as det(G 2 ), . . ., det(G 5 ) are sums of positive monomials and thus, in particular, positive if evaluated at h * (t 1 ).Thus b 7 (t 1 ) and D i (t 1 ), i = 1, 2, . . ., 5 are positive.Theorem 3.4.Consider the dynamical system (25a) -(25j) arising from network (3) with Jacobian J λ (h) as in (30).Let h(t) be as in ( 35) and λ > 0. Fix the remaining h i = h * i such that the inequality (36) is satisfied.Then for t = t 1 computed as in Lemma 3.3, the dynamical system (25a) -(25j) undergoes a simple Hopf bifurcation for all λ > 0 if | h=h(t 1 ) = 0 (by the chain rule) and hence condition (33) is satisfied.Thus, by Corollary 3.1, the dynamical system (1) undergoes a Hopf bifurcation at t = t 1 in this case.Theorem 3.4 establishes the existence of simple Hopf bifurcations for the the dynamical system (25a) -(25j) derived from network (3).In the following remark we observe that inequality (2) and inequality (36) are equivalent: Remark 3.5.First recall that the h i have to satisfy (29) (as (κ, 1 h ) is a solution to the steady state equation ( 10)).That is, h 7 , . . ., h 10 can be represented in terms of κ 3 , κ 6 , κ 9 and κ 12 (and λ) as Using this one obtains for the left hand side of inequality (36): As λ 2 > 0 the inequality (36) is equivalent Finally, we make several remarks regarding the stability of the positive steady state x * (t) = 1 h(t) of the system (25a) -(25j) depending on t > 0.
Remark 3.6.(a) For sufficiently small t > 0, the positive steady state x * (t) is asymptotically stable.(This is true regardless of inequality (36)).
(b) Suppose that the inequality (36) is satisfied.Then for sufficiently large t > 0, x * (t) is unstable.

A procedure to locate simple Hopf bifurcations in network (3)
Following the proof of Lemma 3.3 we proceed as follows to find points that satisfy (32): Step 1: Choose positive values for h 7 , h 8 , h 9 , h 10 such that (36) is satisfied (e.g.h 7 = h 8 = h 10 = 1 and h 9 = 2).
Step 3: In p(t) set For the values chosen so far we obtain 497664 + 10946304t + 103721056t 2 + 579850652t Step 4: Approximate the positive real root(s) of p(t) = 0.For the values chosen so far we obtain t * ≈ 14.874.
We use Matcont to verify a Hopf bifurcation, cf.Fig. 2a.
Step 6: Choose some t > t * , compute vectors h and x = 1 h and use eq.( 14) and ( 26) to obtain rate constants and total concentrations.We have chosen t = 15 and hence obtain h T = (15, 15, 15, 1, 1, 15, 1, 1, 2, 1) and and the rate constants and total concentrations given in Table 2.   33) is satisfied as well.In practice we suggest the following approach: first determine a candidate point x * via Steps 1 to 5, second use ( 29) and ( 26) to determine the corresponding rate constants and total concentrations and third verify the existence of a simple Hopf bifurcation by using a numerical continuation software like MATCONT [18] to vary a rate constant or total concentration at this point.2 with λ = 1.Initial value x(0) as in (37) -apart from x 3 (0) and x 6 (0): to obtain an initial value near the steady state x given in (37) we choose x 3 (0) = 1.1 • 1 15 and x 6 (0) = 0.9 • 1 15 ).
Remark 3.8.In the course of Steps 1 -6 above all rate constants are fixed: (i) From ( 14) one obtains for network (3) the relation (29) between the κ i and the h i .
(iii) Choosing numerical values for h 4 and h 5 and assigning h 1 = h 2 = h 3 = h 6 = t is equivalent to choosing κ 1 , κ 4 , κ 7 and κ 10 (again up to the factor λ).In this case κ 1 and κ 7 are proportional to t 2 and κ 4 and κ 10 to t -as h 4 and h 5 are fixed to numerical values.

Lifting to the full network (1)
In [2] network modifications are described that preserve the existence of a stable positive limit cycle.This is the basis for the following result: Theorem 3.9.Consider networks (1) and ( 3).Assume the rate constant values of network (3) are such that the system (25a) -(25j) admits a stable limit cycle.Choose these values for the rate constants of network (1).For values of κ 2 , κ 5 , κ 8 , κ 11 small enough, there exists a stable limit cycle in the system (24a) -(24j) close to the limit cycle of the system (25a) -(25j).
Proof.In the language of [2], if, a network that admits a stable positive limit cycle (for some values of the rate constants and initial conditions) is modified by adding reactions that are in the span of the stoichiometric matrix, then the new network also admits a stable positive limit cycle [2, Theorem 1] (if the rate constants of the new reactions are chosen appropriately).
As the reaction vectors of reversible reactions are in the span of the stoichiometric matrix, the existence of a limit cycle in network (3) implies the existence of a limit cycle in the full network (1) for appropriately chosen rate constants of the backward reactions.
To illustrate Theorem 3.9, we use the ODEs (24a) -(24j) and the values of Table 2 (on page 16).Fig. 3 demonstrates the existence of a limit cycle in the full system (24a) -(24j)(for k b sufficiently small).
Remark 3.10 (Locating the limit cycle in Fig. 3a).For simplicity we choose κ 2 = κ 5 = κ 8 = κ 11 = k b .To obtain Fig. 3 the initial value given in the table of Fig. 3c was used.This point is 'close' to the limit cycle of the ODEs defined by network (3) (i.e. for k b = 0).It was obtained by solving the ODEs (24a) -(24j) using Matlab's ode15s with initial value given in Table 2 (on page 16) for a 'long' time (i.e. until T = 5000) and λ = 1.The point in the table of Fig. 3c corresponds to the last point of that first simulation.(c) Initial value for simulations displayed in Fig. 3.

Discussion
In this section we discuss inequality (2) in the light of results on multistationarity for a network of sequential distributive double phosphorylation described in [9].In Section 4.1 we introduce the corresponding reaction network and compare it to network (1).In Section 4.2 we briefly summarize the results presented in Section 3 and in Section 4.3 the multistationarity results of [9].We close by arguing our conclusion that in distributive double phosphorylation the catalytic constants enable non-trivial dynamics in Section 4.4.

Cyclic versus sequential distributive double phosphorylation
Sequential and distributive double phosphorylation can be described by the following mass action network (cf.e.g.[16] or [9]): Network ( 38) is structurally similar to network (1): both networks contain 12 reactions and the only difference is that network (1) contains two species of mono-phosphorylated protein (S 10 and S 01 ), while network(38) contains only one (S 1 ).Hence network (38) contains nine species, while network network (1) contains ten.
In particular, networks (1) and (38) contain the same four phosphorylation events: (i) the conversion of unphosphorylated protein to mono-phosphorylated protein catalyzed be the kinase K, (ii) the conversion of mono-phosphorylated protein to double-phosphorylated protein catalyzed by the same kinase K, (iii) the conversion of double-phosphorylated protein to mono-phosphorylated protein catalyzed by the phosphatase F and (iv) the conversion of mono-phosphorylated protein to unphosphorylated protein catalyzed by the same phosphatase F .As described in [9], in enzyme kinetics research it is customary to characterize such phosphorylation events by three constants, the Michaelis constant (K m ), the catalytic constant (k c ) and the equilibrium constant k eq of the respective enzyme substrate pair (see, for example, [3] for details on enzyme kinetics).
Of particular interest in the context of the present publication are the k c -values as these correspond to the rate constants involved in inequality (2): κ 3 is the k c -value of the kinase K with unphosphorylated substrate (S 00 or S 0 ), κ 6 of K with mono-phosphorylated substrate (S 10 or S 1 ), κ 9 of F with double-phosphorylated substrate (S 11 or S 2 ) and κ 12 of F with mono-phosphorylated substrate (S 10 or S 1 ).

Cyclic and distributive: emergence of oscillations
By Theorem 3.4, if these catalytic constants satisfy inequality (2), then there exists positive steady states of network (3) such that the Jacobian has a complex-conjugate pair of eigenvalues on the imaginary axis.This is necessary for a simple Hopf-bifurcation.If there is a supercritical simple Hopf bifurcation and a stable limit cycle emerges, then by Theorem 3.9 there is a stable limit cycle in network (1).Hence we say that for cyclic and distributive double phosphorylation the catalytic constants enable the emergence of oscillations.

Sequential and distributive: emergence of bistability
In [9] we have shown that the inequality (2) is sufficient for multistationarity in network (38).To be more precise, by [9,Theorem 5.1], if the catalytic constants satisfy inequality (2), then there exists values of the total concentrations of kinase, phosphatase and protein such that network (38) has three positive steady states -no matter what values the other rate constants take.As multistationarity is necessary for bistability, we say in [9], that the catalytic constants enable the emergence of bistability in sequential and distributive double phosphorylation.

Catalytic constants and non-trivial dynamics
In the previous subsections we have described how the catalytic constants of cyclic distributive double phosphorylation enable the emergence of oscillations, and how the catalytic constants of sequential distributive double phosphorylation enable the emergence of bistability.Hence we conclude that in distributive double phosphorylation the catalytic constants enable non-trivial dynamics.
As a consequence, if the rate constant are chosen according to the procedure of Section 3.3 and Theorem 3.9 and network (1) admits a stable limit cycle for these rate constants, then network (38) taken with the same rate constant values will show multistationarity -for some, usually different, value of the total concentrations.That is, if the catalytic constants satisfy (2), then a cyclic mechanism can show sustained oscillations, while a sequential mechanism equipped with the same rate constant values can show bistability.
As an example the procedure described in Section 3.3 together with Theorem 3.9 have been used to obtain the following rate constant values: These values satisfy inequality (2).Using these values in the ODEs derived from network (1), one detects simple Hopf bifurcations and oscillations as depicted in Fig. 4a where c 1 denotes total amount of kinase K, c 2 of phosphatase F and c 3 of substrate S.And for Fig. 4c following the results of [9] yields (using the same notation for the total concentrations)  where we have used that the sum of the differences of a cycle sum up to zero by Lemma B.2.Let the sum of the indices of each cycle c 1 , , . . .c f of a permutation σ be denoted by δ 1 , . . .,δ f , correspondingly.We use the fact that the cycles {c 1 , . . ., c f } in a permutation σ are pairwise disjoint.Thus it follows that By Lemma B.1, (47) and (48) we have Thus, det H l (λ, h) = λ l(l+1)/2 det G l (h) for l = 1, 2, . . ., s.

D Initial data for Fig. 4a & 4b
To generate Fig. 4a & 4b the ODEs derived from the full (reversible) reaction network (1) are used (i.e. the ODEs (24a) -(24j)).In both Fig. 4a and 4b we use the point displayed in Table 4 as initial value.
x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 x 10 1 7 1 7 1 7 1 1 , h 4 = h 5 = 1 and t = 1).We use Matlab's ode15s to solve the initial value problem defined by the above ODEs with initial value given in Table 4 and the backward constants κ 2 = κ 5 = κ 8 = κ 11 = 1 10 (and the remaining constants according to our procedure for the irreversible network).These backward constants are 'small enough' in the sense of the results of [2], that is, that the reversible network (1) has a limit cycle close to the limit cycle of the irreversible network.
However, on the one hand the steady state of Table 4 (of the irreversible network) is 'far enough' from a steady state of the reversible network in the following sense: the solution of the reversible system with initial value given in Table 4 approaches a limit cycle of the reversible system (if approximated with ode15s).If the point given in Table 4 was 'too close' to a steady state of the reversible system the solution with it as initial value would approach that steady state (if approximated with ode15s).
On the other hand, Fig. 4a is generated with Matcont using the point given in Table 4 as an initial guess of a steady state of the above ODEs (of the reversible network).And this point is close enough to a steady state of the reversible network for Matcont to converge to a steady state.

E Initial data for Fig. 4c
To obtain Fig. 4c we follow [9], where a system of nine ODEs is derived from network (38).As in this reference, we use the variables given in Table 5 to denote the species concentrations.As this network does not distinguish S 10 and S 01 , there is only one mono-phosphorylated from of the substrate (S 1 ) and S 2 is used to denote the double-phosphorylated substrate.Consequently, this network contains only nine species.It has, however, 12 reactions and the labeling is consistent with network (1).
Using the approach described in [9] we obtain the steady state depicted in Table 6.We use this point as a starting point for the numerical continuation in Matcont.Table 6: Starting point for the numerical continuation in Fig. 4c (a steady state of network (38) for the rate constant values of eq.(39) generated as described in [9]).

Figure 1 :
Figure 1: Reaction schemes describing distributive double phosphorylation of a protein S at two binding sites by a kinase K and a phosphatase F .Panel (a) cyclic double phosphorylation, panel (b) distributive.In panel (a) the subscript ij denotes the state of the phosphorylation sites: 0 unphosphorylated, 1 phosphorylated (e.g. S 10 denotes those molecules of S, where the first site is phosphorylated and the second site is unphosphorylated).In panel (b) the subscript i denotes the number of attached phosphate groups (e.g. S 0 denotes unphosphorylated protein).

Proof. If t = t 1
is as in Lemma 3.3, the conditions (32) of Corollary 3.1 are satisfied.If d det(G 6 ) dt | t=t 1 = 0, then at least one of the derivatives ∂ det(G 6 )∂h i

Figure 3 :
Figure 3: Simulation of network (1) for κ 2 = κ 5 = κ 8 = κ 11 = k b and different values k b (ODEs have been solved with ode15s (Mathworks) for x(0) as in the table of panel (c) and κ i , c i as in Table 2 (on page 16 with λ = 1).Panel (a): k b = 0 corresponds to the ODEs derived from network (3), k b = 0.05 to the ODES (24a) -(24j) for k b = 0.05.The oscillations indicate for k b = 0.05 a stable limit cycle close to the stable limit cycle for k b = 0. Panel (b): the stable limit cycle does not seem to exist for larger values of k b .
and Fig. 4b.And using these values in the ODEs derived from network (38), one obtains multistationarity as depicted in Fig. 4c.To create these figures, the same parameter values have been used in both ODE systems, albeit for different values of the total concentrations.For Fig. 4a and Fig. 4b the procedure of Section 3.3 yields Cyclic: numerical continuation of steady states.Steady state value of S 11 (x 6 ) vs. total concentration of Kinase K (c 1 ); Hopf points (H).Cyclic: numerical solutions of ODEs defined by network (1); plotting S 00 (x 3 ) and S 11 (x 6 ) vs. time t shows sustained oscillations.Sequential: numerical continuation of steady states; S 2 (x 6 ) vs. c 1 ; LP limit point, multistationarity for c 1 between LPs.

Table 2 :
(13) constants and total concentrations obtained by solving(13)for k.Remark 3.7.As a consequence of Lemma 3.3, any point x * obtained via Steps 1 to 5 is a candidate Hopf point.It is guaranteed to satisfy (32), but a simple Hopf bifurcation only occurs if the condition (