On the existence of Hopf bifurcations in the sequential and distributive double phosphorylation cycle

: Protein phosphorylation cycles are important mechanisms of the post translational modiﬁcation of a protein and as such an integral part of intracellular signaling and control. We consider the sequential phosphorylation and dephosphorylation of a protein at two binding sites. While it is known that proteins where phosphorylation is processive and dephosphorylation is distributive admit oscillations (for some value of the rate constants and total concentrations) it is not known whether or not this is the case if both phosphorylation and dephosphorylation are distributive. We study simpliﬁed mass action models of sequential and distributive phosphorylation and show that for each of those there do not exist rate constants and total concentrations where a Hopf bifurcation occurs. To arrive at this result we use convex parameters to parametrize the steady state and Hurwitz matrices.


Introduction
Protein phosphorylation cycles consist of three proteins, a substrate S and two enzymes K and F. One enzyme, the kinase K, attaches phosphate groups to the substrate and hence phosphorylates the substrate while the other, the phosphatase F, removes phosphate groups and hence dephosphorylates the substrate. Protein phosphorylation cycles are important mechanisms of post translational modification of a protein and as such an integral part of intracellular signaling and control [1]. Often phosphorylation and dephosphorylation follow a sequential and distributive mechanism as depicted in Figure 1a: in each encounter of S and either K or F exactly one binding site is (de)phosphorylated. If either phosphorylation or dephosphorylation follows a processive mechanism, then at least two binding sites are (de)phosphorylated in each encounter of S and either K or F (cf. Figure 1b & 1c). Here we study the sequential and distributive phosphorylation of a protein S at two binding sites as depicted in Figure 1a. Such a study of the behavior of an important biochemical module is of particular interest in the light of studies elucidating the complex behavior of signaling pathways composed of such modules [2].
Mathematical models of both processive and distributive phosphorylation have been extensively studied and it is known that they admit complex dynamics (see e.g. [1,[3][4][5] and the many references therein). The mass action model of the sequential and distributive phosphorylation cycle depicted in Figure 1a is arguably one of the -mathematically -best studied and challenging systems of post translational modification: both multistationarity (the existence of at least two positive steady states) and bistability (the existence of two locally stable positive steady states) have been established (cf. for example [6,7] for multistationarity, [8] for bistability). In fact it has been shown that this mass action model admits at most three positive steady states [9,10].  In contrast purely processive phosphorylation cycles have a unique stable steady state [11], while the mixed cycles depicted in Figure 1b and 1c have a unique steady state that might not be stable, and admit oscillations [12,13]. Distributive steps therefore seem to be involved with the emergence of oscillations, in particular as in more involved combinations of distributive and processive steps, oscillations have been reported as well [1,14].
Interestingly, oscillations have not been reported for the cycle depicted in Figure 1a, despite considerable effort by different research groups. One avenue to establish oscillations is to determine values of rate constants and concentration variables where a supercritical Hopf bifurcation occurs, see e.g., [15]. In [16] Hopf bifurcation points have been located in the parameter space of a variety of models from the systems biology literature. In [16,Section 5.36] a simplification of the mass action model of Figure 1a is examined. The authors provide steady state concentration values and rate constants of a candidate Hopf bifurcation point. (The rate constants are provided indirectly since they can be computed from the 'convex parameters' the authors use, cf. Section 3.2). Applying Theorem 1, however, shows that at this point no Hopf bifurcation occurs: the candidate point satisfies the conditions on the Hurwitz determinants but it fails the condition on the constant coefficient of the characteristic polynomial (in our case the last nonzero coefficient of the characteristic polynomial). Thus it cannot be a point of Hopf bifurcation, for more details see the discussion in Section 6. In summary, to the best of our knowledge, for the mass action model of the phosphorylation cycle in Figure 1a nor simplifications of it, neither Hopf bifurcations nor oscillations have been reported to date.
In this paper we work towards solving the problem of existence or non-existence of Hopf bifurcations for the double phosphorylation network. We follow the strategy outlined in [12] and as in [16] we use convex parameters (see Section 3.2). As in [16] we consider simplifications of the mass action model derived from Figure 1a. Specifically we consider all four mass action networks derived from Figure 1a that contain only two (out of four possible) enzyme-substrate complexes. This is done for two reasons: first, the four networks are biochemically interesting and hence worth studying by themselves (cf. Remark 5). Second, while it is trivial to see that three of the models do not admit Hopf bifurcations, the fourth model displays a nice structure that can be explored, even when all parameters treated as unknown, by performing symbolic computations. Treatment of the full model is currently out of reach due to the computational complexity.
The paper is organized as follows: in Section 2 we introduce the notation. We further recall convex parameters and a criterion for purely imaginary roots in Section 3. In Section 4 we motivate the mass action models studied in this paper. In Section 5 we prove that for none of the four mass action models Hopf bifurcations are possible: for each mass action model we show that there do not exist parameter values such that a Hopf bifurcation occurs (this is Theorem 2). The proof of the theorem relies on large symbolic computations that we present in the Maple Supplementary File 'SupplMat1.mw' (see also 'SupplMat1.pdf' for a pdf version). In Section 6 we comment in more detail on the candidate point for Hopf bifurcation presented in [16].

Notation
We consider systems of n chemical species A 1 , . . . , A n and r chemical reactions of the form where the integer numbers α i j ≥ 0, β i j ≥ 0 are the stoichiometric coefficients and κ j > 0 the rate constants. We use x i to denote the concentration of species A i . Throughout this paper we will assume that all reactions are governed by mass action kinetics, that is, the reaction rate is proportional to the product of the concentrations of the reacting species raised to the power of their respective molecularities. Then the reaction rate v j (κ, x) of the j-th reaction is With this, the above reaction network defines the following system of ordinary differential equationṡ where N is the stoichiometric matrix. Here N is of dimension n × r and v(κ, x) of dimension r × 1. The i j-th entry of N is given by the difference of the stoichiometric coefficients: Let diag(κ) be the diagonal r × r matrix of rate constants where the κ i coordinate of an r-dimensional vector κ is the [i, i] entry of diag(κ). Let Y be the n × r matrix whose column vectors y j contain the stoichiometric coefficients α i j of the reactant of the j-th reaction. The matrix Y is sometimes called the kinetic order matrix. Given vectors x, y ∈ R n , we use the notation x y = n i=1 x y i i . Then the columns of Y define the monomial vector and v(κ, x) can be written as the product v(κ, x) = diag(κ)ψ(x).
For example the reaction network consists of the three reactions among the four species S 0 , K, S 1 and S 0 K taken in that order as species A 1 through A 4 . The stoichiometric coefficients of the first reaction are α 11 = α 21 = β 41 = 1 and α 31 = α 41 = β 11 = β 21 = β 31 = 0. The differences β i1 − α i1 define the first column of the following stoichiometric matrix (the remaining columns are defined in a similar way): With mass action kinetics, the reaction rate vector is The kinetic order matrix is

A theorem to preclude Hopf bifurcations
In this subsection we state a criterion in terms of the principal minors of the Hurwitz matrix (Proposition 1, Theorem 1) that we will use to exclude Hopf bifurcations. The criterion follows from the results in [17]. Related results for asserting Hopf bifurcations can be found in [18,19].
For ease of notation consider an ODE system parametrized by a single parameter µ ∈ R: where x ∈ R s , and g µ (x) varies smoothly in µ and x. Let x * ∈ R s be a steady state of the ODE system for some fixed value µ 0 , that is, g µ 0 (x * ) = 0. Furthermore, we assume that we have a smooth curve of steady states around µ 0 : that is, g µ (x(µ)) = 0 for all µ close enough to µ 0 such that x(µ 0 ) = x * . By the Implicit Function Theorem, this curves exists if the Jacobian of g µ 0 (x) evaluated at x * is non-singular. Let J(x(µ), µ) be the Jacobian of g µ (x) evaluated at x(µ). If, as µ varies, a single complex-conjugate pair of eigenvalues of J(x(µ), µ) crosses the imaginary axis while all other eigenvalues remain with nonzero real parts, then a simple Hopf bifurcation occurs at (x(µ), µ). In this case, a limit cycle arises. If the Hopf bifurcation is supercritical and all other eigenvalues have negative real part, then stable periodic solutions are generated for nearby parameter values.
Remark 1. As we pointed out in the informal discussion above, a simple Hopf bifurcation requires that exactly one pair of complex conjugate eigenvalues crosses the imaginary axis when the parameter µ varies. Thus, in particular, there must exist a value µ 0 with steady state x(µ 0 ), where the Jacobian J(x(µ 0 ), µ 0 ) has exactly one pair of purely imaginary eigenvalues. More generally, a Hopf bifurcation necessitates J(x(µ 0 ), µ 0 ), for some value µ 0 , to have a pair of purely imaginary eigenvalues.
The following proposition gives necessary and sufficient conditions for the existence of exactly one pair of purely imaginary eigenvalues in the scenario we encounter later on. It is based on Hurwitz matrices, which we define first: Definition 1. 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. Note H s = a s H s−1 .
We now state the criterion we will use to determine whether the characteristic polynomial has a pair of purely imaginary roots. Proposition 1. In the setup above, let s ≥ 2, µ 0 be fixed, and let H i (µ 0 ) denote the i-th Hurwitz matrix of p µ 0 (z). We assume that det H 1 (µ 0 ) > 0, . . . , det H s−2 (µ 0 ) > 0.
Then p µ 0 (z) has at most one pair of symmetric roots, and has exactly one if and only if det H s−1 (µ 0 ) = 0.
In this case, the pair consists of purely imaginary roots if and only if a s (µ 0 ) > 0.
Proof. The first statement follows from Corollary 3.2 in [17] using (4). So, assume p µ 0 (z) has a pair of symmetric roots z, −z such that det H s−1 (µ 0 ) = 0. By Lemma 3.3 in [17], the Routh-Hurwitz criterion for stable polynomials, and using (4), we conclude that all the other s − 2 roots of p µ 0 (z) have negative real part. Now a s (µ 0 ) is (−1) s times the product of the roots of p µ 0 (z), and hence its sign agrees with the sign of (−1) s (−1) s−2 z(−z) = −z 2 . We conclude that the pair of symmetric roots of p µ 0 (z) is real if a s (µ 0 ) ≤ 0 and purely imaginary if a s (µ 0 ) > 0.
As a consequence of Propostion 1, we obtain the following criterion to preclude Hopf bifurcations: Theorem 1. Consider the dynamical systemẋ = g(x) with s ≥ 2 and assume there exists a curve of steady states x(µ) as in the above setting. As before, let p µ (z) be the characteristic polynomial of the Jacobian J( for all µ and either then the system does not undergo a (simple) Hopf bifurcation.
Proof. Follows from Remark 1 and Proposition 1.
In the following sections we will use Theorem 1 to prove the non-existence of Hopf bifurcations in subnetworks of the mass action network derived from Figure 1a. For a reaction network with n species, if the rank s of the stoichiometric matrix N is not maximal, as it is the case for our networks, then the dynamics takes place in the invariant s-dimensional linear subspaces x 0 + im N. This implies that 0 is a root of the characteristic polynomial of Nv(κ, x) with multiplicity n − s and hence it factors as The polynomial a 0 (κ, x)z s + a 1 (κ, x)z s−1 + · · · + a s (κ, x) is the characteristic polynomial of the Jacobian of the restriction of system (2) to x + im N, and hence we apply Theorem 1 to this polynomial.

Convex parameters
By the previous subsection, in order to determine whether a Hopf bifurcation can arise in our systems, we need to analyze the Jacobian matrix of the right-hand side of (2) for all possible values of rate constants κ and positive steady states x * . Here we reparametrize the Jacobian matrices using so-called convex parameters. These parameters were introduced by Clarke in [20] to analyze the stability of mass action reaction systems (2), in the context of Stoichiometric Network Analysis (SNA) theory.
Let a positive steady state x * ∈ R n >0 of (2) satisfy the polynomial system N diag(κ)ψ(x * ) = 0. Then the rate functions v = v(κ, x * ) satisfy the linear problem The vector v is referred to as a flux vector in the SNA theory [21]. The solutions v of (5) define a convex polyhedral cone called the flux cone. Convex polyhedral cones have a finite number of extreme vectors up to a scalar positive multiplication [22]. Therefore, any flux vector v can be represented as a nonnegative linear combination of its extreme vectors where E is the matrix with columns E 1 , . . . , E l and λ = (λ 1 , . . . , λ l ).
Remark 2 (cf. [22]). The nonnegative coefficients λ 1 , . . . , λ l in (6) are often referred to as convex parameters in the literature. However, they do not account alone for all new parameters -the other group of parameters used in SNA theory are reciprocals of each positive steady state coordinate x * k > 0. They are denoted by Definition 2. A vector of convex parameters is a vector of the form (h, λ) = (h 1 , . . . , h n , λ 1 , . . . , λ l ) ∈ R n >0 × R l ≥0 such that Eλ ∈ R r >0 . The convex parameters are convenient for parameterizing the Jacobian J(κ, x) evaluated at a positive steady state x = x * . To see this, note that the ( j, i)-th entry of the Jacobian of v(κ, x) evaluated at x * is Hence, the Jacobian of Nv(κ, x) evaluated at x * is where we use the vector notation Using now (6) and (7) to write v(κ, x * ) = Eλ, the Jacobian of Nv(κ, x) evaluated at x * can be written as with (h, λ) a vector of convex parameters. Therefore, given a vector of rate constants κ and a corresponding positive steady state x * , there exist convex parameters (h, λ) such that equality (8) holds. Conversely, given convex parameters (h, λ), we define x * = 1/h and let which is a positive vector since all entries of Eλ are positive. Then, using that (8) holds as well. This proves the following proposition: Proposition 2. The set of Jacobian matrices J(κ, x * ) for all κ and corresponding positive steady states x * agrees with the set of matrices defined by the right-hand side of (8), for all h ∈ R n >0 and λ ∈ R l ≥0 such that Eλ ∈ R r >0 . The computation of the Jacobian in convex parameters (8) appears in great detail in previous works [23,24]. In [25] it is used to detect saddle-node bifurcations. In Section 5, we use the Jacobian in convex coordinates given in (8) and apply Theorem 1 to conclude that there does not exist a point (h,λ) where a Hopf bifurcation occurs. Using the equality between the two sets of matrices in Proposition 2, this will imply that there do not exist κ and x * where a Hopf bifurcation occurs.   At the level of mass action kinetics each of these phosphorylation events can be described by the following reactions (with i = 0, 1): Consequently, if all phosphorylation events depicted in Figure 1a are described at the mass action level one obtains the following reaction network: To apply Theorem 1, one has to compute Hurwitz determinants (see Section 3.1 above). These are determinants of matrices that are composed of the coefficients of the characteristic polynomial of the Jacobian of the ODEs defined by N. In order to show whether or not there exist some values of the rate constants where a Hopf bifurcation occurs, we have to treat all rate constants as fixed but unknown. The coefficients of the characteristic polynomial may contain several hundred terms (cf. the supporting information of [7]). To facilitate the analysis we consider the following simplifications of N: (i) We consider only the forward reaction of the reversible reactions This is a reasonable assumption if the rate constants for the reversible reactions are small.
(ii) We consider only two of the four enzyme-substrate complexes KS 0 , KS 1 , FS 2 and FS 1 .
There are six ways to choose two complexes out of four. Due to the symmetry of the ODE system obtained by interchanging K and F, S 0 and S 2 , and relabeling the rate constants as appropriate, it suffices to consider the following four simplified networks derived from N: • The network containing only KS 0 and FS 2 : • The network containing only KS 0 and FS 1 (mathematically equivalent to the network containing only FS 2 and KS 1 ): • The network containing only KS 0 and KS 1 (mathematically equivalent to the network containing only FS 2 and FS 1 ): (N 3 ) • The network containing only KS 1 and FS 1 : Specifically, under this assumption, the simplified system corresponds to the slow system arising from Tikhonov-Fenichel singular perturbation reduction of the original system, which additionally agrees with its quasi-steady-state reduction [26]. Upon removal of a reverse unbinding reaction, the assumption is that this reaction occurs at a much slower time scale than the other reactions in the network, and hence is negligible.

Remark 5.
A biochemical interpretation of the simplification leading to the networks (N 1 ) -(N 4 ) can be obtained by comparing it to the well-known Michaelis-Menten approximation: we view our simplification as similar in spirit but with different focus and hence concurrent. To see this recall that to simplify a reaction network based on the Michaelis-Menten approximation, a mass action network of the form is replaced by a network of the form where v MM is the familiar Michaelis-Menten kinetics with catalytic constant k cat , total enzyme concentration K 0 , Michaelis-Menten constant K m , and where [·] denotes the concentration. This is a standard practice of model simplification, for example if the condition formulated by Briggs and Haldane holds [27]; see [28,29] for the mathematical study of this reduction. In [30] it is shown that the dynamical system that arises from the Michaelis-Menten approximation of the full network N does not exhibit periodic orbits. In our simplification we replace a mass action network of the form (9) by a mass action network of the form where, according to the law of mass action the reaction rate is In summary, the simplification based on the Michaelis-Menten approximation on the one hand is well suited to describe the saturation of the enzyme with substrate but does not account for the influence of varying concentrations of free enzyme. Our simplification on the other hand cannot account for enzyme saturation, but for the influence of varying concentrations of free enzyme. Moreover, for small values of the substrate concentration and for concentrations of free enzyme close to its total amount both simplifications behave similar. Hence we view our simplification as biochemically concurrent to the one based on the Michaelis-Menten approximation. Both simplifications fail to accommodate the complete behavior of the distributive and sequential double phosphorylation cycle of Figure 1a. But both cover different but equally important aspects of its behavior and hence are well worth studying. Remark 6. If any of the networks (N 1 ) -(N 4 ) had a Hopf bifurcation giving rise to periodic solutions, then by [31], so would the full mechanism in N. In particular, by [31] the same is true if any of the irreversible reactions is made reversible and in particular, if unbinding reactions are considered in the formation the complexes KS 0 , KS 1 , FS 1 or FS 2 .

Absence of Hopf bifurcations
In this section we apply Theorem 1, using the discussion after it, to the networks (N 1 ) -(N 4 ). To this end we use Eq (8) to determine the Jacobian matrices J 1 (h, λ), . . . , J 4 (h, λ) of networks (N 1 ) -(N 4 ), correspondingly. The computations are done symbolically and can be found in the supporting file 'SupplMat1.mw'.
We first comment on some common features of the networks: Observe that a scaling αλ of the vector λ = (λ 1 , λ 2 ) with α > 0 translates into a scaling ακ of the vector of rate constants κ under the correspondence of parameterization in Proposition 2. Further x * is a steady state for κ if and only if it is for ακ and αJ(κ, x * ) = J(ακ, x * ). The latter implies that any eigenvalue of J(ακ, x * ) is in fact an eigenvalue of J(κ, x * ) multiplied by α > 0. Thus, J(ακ, x * ) has a pair of purely imaginary eigenvalues if and only if J(κ, x * ) has a pair of purely imaginary eigenvalues.
Hence, it is enough to take one of λ 1 , λ 2 to be one (since both are positive), and we let the elements in the kernel of the stoichiometric matrices N i be of the form Therefore, we consider every Jacobian matrix J i (h, λ) according to Eq (8) parametrized by 8 parameters: the parameter λ and h 1 , . . . , h 7 . By Remark 7 (i) and (ii) it follows that the characteristic polynomial of every Jacobian J i (h, λ) is a degree 7 polynomial of the form where each a i depends on the 8 parameters λ, h 1 , . . . , h 7 (cf. the file 'SupplMat1.mw'). Following the discussion after Theorem 1, for each network we compute a 0 (h, λ), . . . , a 4 (h, λ), det H 1 (h, λ), det H 2 (h, λ) and det H 3 (h, λ) (cf. Definition 1) and show the following proposition: In particular, this proposition (which is proven below) tells us that in all four networks, det H 3 (h, λ) 0 whenever a 4 (h, λ) > 0. As a consequence we obtain the following theorem.
Theorem 2. For the networks (N 1 ) -(N 4 ) there do not exist rate constants κ and a corresponding positive steady state x * such that the Jacobian matrix J i (κ, x * ) has a pair of purely imaginary eigenvalues. Thus, in particular, there do not exist κ and x * where a Hopf bifurcation occurs.
Proof. By Proposition 3 and the correspondence between the two parameterization of the Jacobian given in Proposition 2, there do not exist κ and a corresponding positive steady state x * such that the corresponding Hurwitz determinant det H 3 (κ, x * ) vanishes and the coefficient of lowest degree of the characteristic polynomial a 4 (κ, x * ) is positive. Further, for all networks det H 1 (κ, x * ) > 0 and det H 2 (κ, x * ) > 0 for all κ, x * . Hence, by Theorem 1, the Jacobian matrix does not have a pair of purely imaginary eigenvalues.
Remark 8. By [33], the networks N 1 , N 2 and N 3 are multistationary, while N 4 is not. Hence, by [32], the coefficient of lowest degree of the Jacobian must vanish for some κ and positive steady state x * , and consistently a 4 (h, λ) must contain monomials of both signs.
Remark 9. We have focused on networks with two intermediates, as network ((N 1 )) provides the first non-trivial case and might shed light on how to approach the full network. All simplifications of the full network (N) obtained after removing three intermediates, that is, with only one intermediate left, satisfy that det H s−1 (h, λ) is a sum of positive terms and hence does not vanish. Consequently, Hopf bifurcations do not arise.
All that remains is to show Proposition 3. This is done through mathematical reasoning aided by symbolic computations performed in Maple and Mathematica. In the following subsections we present for each network the stoichiometric matrix N i , the kinetic order matrix Y i , the matrix E i whose columns are vectors w 0 and w 1 that generate the cone (5) and the Jacobian matrix J i (h, λ). For the networks ((N 1 )) -((N 4 )) we found these vectors simply by finding a basis of the kernel of N, which has dimension 2, and deriving extreme vectors. For larger networks software like CellNetAnalyzer [34] can be used. Where appropriate we then discuss the coefficient a 4 (h, λ) and the determinants of the Hurwitz matrices det H i (h, λ). The computations are given in the supplementary file 'SupplMat1.mw'.

Network (N 1 )
We denote by x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 the concentrations of K, F, S 0 , S 1 , S 2 , KS 0 , FS 2 respectively. Then the stoichiometric matrix, the kinetic order matrix and the matrix E 1 are Using convex parameters, the Jacobian matrix in terms h 1 , . . . , h 7 , λ as given in (8) is We observe that the coefficient a 4 (h, λ) of the characteristic polynomial contains monomials of both signs. We compute the associated Hurwitz determinants det H 1 (h, λ), det H 2 (h, λ) and det H 3 (h, λ) and obtain that det H 1 (h, λ) and det H 2 (h, λ) are sums of positive monomials and that det H 3 (h, λ) contains monomials of both signs as well. Hence both a 4 (h, λ) and det H 3 (h, λ) contain monomials of both signs and can potentially be zero.
In the remainder of this section we prove the following: Proof. The coefficient a 4 (k, λ) of the Jacobian is Since λ 2 factors out and it does not affect the sign, we consider We show that c 0 > 0 implies det H 3 > 0, omitting the argument of H 3 . The computations are performed in Maple, but we explain here the computational procedure for the proof.
We start by noting that c 0 can be written as: We see immediately that if h 1 > h 6 and h 2 > h 7 , then c 0 < 0. We do not need to study this case.
We consider the case h 6 ≥ h 1 and h 7 ≥ h 2 , such that one of the two inequalities is strict, otherwise c 0 = 0. In this case c 0 ≥ 0. We introduce new nonnegative parameters v 1 , v 2 and substitute h 6 = h 1 + v 1 and h 7 = h 2 + v 2 . This encodes the inequalities. We perform this substitution into det H 3 using Maple, expand the new polynomial, and check the sign of the coefficients. All coefficients in det H 3 are positive, meaning that det H 3 will be positive in this case. This holds if v 1 , v 2 > 0 or if one of the two parameters is set equal to zero.
The latter means that we need to study the case h 1 ≥ h 6 and h 2 ≤ h 7 . We now perform the substitution h 1 = h 6 + v 1 and h 7 = h 2 + v 2 , again with one of v 1 or v 2 nonzero. We observe that if h 5 ≥ h 3 , then again det H 3 is positive. So we restrict to h 3 > h 5 and perform the substitution When we do that, det H 3 has coefficients of both signs, and therefore the sign is not clear. We have still to impose c 0 > 0 for this scenario. We perform the substitutions into c 0 and obtain: For c 0 to be positive, we need v 1 to be smaller than the root of c 0 seen as a polynomial in v 1 , which is: . Now, we need to check whether det H 3 can be negative when v 1 is smaller than z 1 . To check that, we make the substitution v 1 = µ µ + 1 z 1 .
Then any number in the interval [0, z 1 ) is of this form for some nonnegative µ. We perform this substitution in Maple, gather the numerator of the resulting det H 3 , and confirm that all signs are positive. Further det H 3 is positive even if some of µ, v 2 , v 3 are zero. The other case h 2 ≥ h 7 and h 1 ≤ h 6 is analogous by the symmetry of the system. This finishes the argument, since we have explored all possibilities for c 0 > 0, and they all give that det H 3 is positive.

Network N 2
We use the following ordering: x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 for the concentration of K, F, S 0 , S 1 , S 2 , KS 0 , FS 1 respectively. Under this ordering the stoichiometric matrix, the kinetic order matrix and the matrix E 2 are With this parametrization, the Jacobian of the system evaluated at a steady state defined by (h 1 , . . . , h 7 , λ) is: det H 3 (h, λ) contains only positive monomials and in particular it does not vanish for any positive h, λ.

Network N 3
We use the following ordering: x 7 for the concentration of K, F, S 0 , S 1 , S 2 , KS 0 , KS 1 respectively. Under this ordering the stoichiometric matrix, the kinetic order matrix and the matrix E 3 are With this parametrization, the Jacobian of the system evaluated at a steady state defined by (h 1 , . . . , h 7 , λ) is: det H 3 (h, λ) contains only positive monomials and in particular it does not vanish for any positive h, λ.

Network (N 4 )
We use the following ordering: x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 for the concentration of K, F, S 0 , S 1 , S 2 , KS 1 , FS 1 respectively. Under this ordering the stoichiometric matrix, the kinetic order matrix and the matrix E 4 are With this parametrization, the Jacobian of the system evaluated at a steady state defined by (h 1 , . . . , h 7 , λ) is: det H 3 (h, λ) contains only positive monomials and in particular it does not vanish for any positive h, λ.
6. Discussion and Outlook: Does the full network admit Hopf bifurcations?
As discussed in the Introduction, in [16, Section 5.36] a simplification of the mass action model of Figure 1a is examined and the authors provide steady state concentration values and rate constants of a candidate Hopf bifurcation point. Here we want to explain how this point fails to give a pair of purely imaginary eigenvalues (Proposition 1).
While det H 4 (h * , λ * ), det H 5 (h * , λ * ) and a 6 (h * , λ * ) are all very small numbers, det H 5 (h * , λ * ) is by orders of magnitude smaller. And if 10 −42 ≈ 0, then det H 4 (h * , λ * ) and a 6 (h * , λ * ) are nonzero. In particular a 6 (h * , λ * ) < 0 and hence this point gives rise to a pair of symmetric real eigenvalues. We have additionally verified our claim by considering the set of corresponding eigenvalues (computed with Mathematica): All eigenvalues are real numbers, three eigenvalues are ≈ 0 (as expected as the system has three conservation relations) and the real values 0.0000339721 and −0.0000339721 give a pair of symmetric eigenvalues (up to 10 digits). Hence, this justifies our claim that the question of whether or not a Hopf bifurcation can occur in the full network (N) is still open. For the full network, det H 1 , . . . , det H s−2 are all positive, and hence the problem is reduced to making det H s−1 vanish while having a s positive. Network (N 1 ) is the smallest simplified network where both det H s−1 and a s attain both signs and hence poses the same challenge.
For networks admitting pairs of purely imaginary eigenvalues, a typical approach in this situation is to find values for which det H s−1 vanishes, and then verify that a s is positive when evaluated at these values. But if the two conditions happen to be incompatible, as it might be the case for the full network (N), then it is unclear how to proceed. Automated approaches [16] cannot handle the analysis of the signs of these two large polynomials at the same time. Here, by combining guided heavy symbolic computations with manual intervention based on the visual inspection of a s , we have managed to prove that the two sign conditions are incompatible for network (N 1 ).
The idea is to reduce the problem into checking the sign of the coefficients of a polynomial. We first break the condition a s > 0 into subcases involving simple inequalities between the parameters. Then, one uses the inequalities to substitute one parameter in H s−1 by a new parameter that is allowed to take any positive value (e.g. k 1 > k 2 leads to the substitution k 1 = k 2 + a, such that the new constraints are k 2 , a > 0). Finally, H s−1 is transformed into a polynomial (or rational function) that only needs to be evaluated at positive values to guarantee a s > 0. If all coefficients of H s−1 are positive, then trivially the polynomial is positive. If some coefficients are negative, then further exploration is required. This strategy has worked nicely for network (N 1 ). However, substitutions of the type k 1 = k 2 + a typically increase dramatically the number of terms as a polynomial of the form k d 1 p(κ, x) becomes (k 2 +a) d p(κ, x). We also performed substitutions of the form k 1 = µ µ+1 p(κ,x) q(κ,x) with µ > 0, which introduced large denominators.
All this illustrates the difficulty in the analysis of the full network: the computation of the Hurwitz determinants is very demanding, but a posterior analysis, in par with the analysis of network (N 1 ), is prohibitive. Nevertheless, we think ideas introduced here in the analysis of network (N 1 ) might be applied to other networks posing similar challenges, and maybe similar ideas end up being successful for the full network either after increasing computational power or by depicting new strategies to simplify computations.
As explained in Remark 6, the existence of periodic orbits in simplified models can be lifted to more complex models under certain network modifications [31]. However, the non-existence of periodic orbits or Hopf bifurcations does not, in principle, tell us anything about complex models including the simplified models in some form. Considering the demanding computational cost involved in describing the dynamics of a reaction network for arbitrary parameter values, it would be of great help to have a better understanding of what type of network operations or parameter regimes guarantee that the nonexistence of a behavior in a simplified model is preserved in the more complex model. For example, we are not aware of any result nor reasoning we could employ to ensure that our conclusions on networks (N 1 )-(N 4 ) are preserved after introducing unbinding reactions with a small rate constant.