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

1.   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 ^[10,9].

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 5.1). 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].

2.   Notation

We consider systems of $n$ chemical species $A_1 , \ldots, A_n$ and $r$ chemical reactions of the form

where the integer numbers $\alpha_{ij} \geq 0$, $\beta_{ij}\geq 0$ are the stoichiometric coefficients and ${\kappa}_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({\kappa}, x)$ of the $j$-th reaction is

With this, the above reaction network defines the following system of ordinary differential equations

where $N$ is the stoichiometric matrix. Here $N$ is of dimension $n \times r$ and $v({\kappa}, x)$ of dimension $r \times 1$. The $ij$-th entry of $N$ is given by the difference of the stoichiometric coefficients:

Let ${\rm{diag}}(\kappa)$ be the diagonal $r \times r$ matrix of rate constants

where the ${\kappa}_i$ coordinate of an $r$-dimensional vector ${\kappa}$ is the $[i, i]$ entry of ${\rm{diag}}({\kappa})$. Let $Y$ be the $n \times r$ matrix whose column vectors $y_j$ contain the stoichiometric coefficients $\alpha_{ij}$ of the reactant of the $j$-th reaction. The matrix $Y$ is sometimes called the kinetic order matrix. Given vectors $x, y\in \mathbb{R}^n$, we use the notation $x^y = \prod_{i = 1}^n x_i^{y_i}$. Then the columns of $Y$ define the monomial vector

and $v({\kappa}, x)$ can be written as the product

For example the reaction network

consists of the three reactions

among the four species $S_0$, $K$, $S_1$ and $S_0K$ taken in that order as species $A_1$ through $A_4$. The stoichiometric coefficients of the first reaction are $\alpha_{11} = \alpha_{21} = \beta_{41} = 1$ and $\alpha_{31} = \alpha_{41} = \beta_{11} = \beta_{21} = \beta_{31} = 0$. The differences $\beta_{i1}-\alpha_{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

3.   Background

3.1.   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 $\mu \in \mathbb{R}$:

where $x \in \mathbb{R}^s$, and $g_{\mu}(x)$ varies smoothly in $\mu$ and $x$. Let $x^* \in \mathbb{R}^s$ be a steady state of the ODE system for some fixed value $\mu_0$, that is, $g_{\mu_0}(x^*) = 0$. Furthermore, we assume that we have a smooth curve of steady states around $\mu_0$:

that is, $g_{\mu}\left(x(\mu) \right) = 0$ for all $\mu$ close enough to $\mu_0$ such that $x(\mu_0) = x^*$. By the Implicit Function Theorem, this curves exists if the Jacobian of $g_{\mu_0}(x)$ evaluated at $x^*$ is non-singular.

Let $J(x(\mu), \mu)$ be the Jacobian of $g_\mu(x)$ evaluated at $x(\mu)$. If, as $\mu$ varies, a single complex-conjugate pair of eigenvalues of $J(x(\mu), \mu)$ crosses the imaginary axis while all other eigenvalues remain with nonzero real parts, then a simple Hopf bifurcation occurs at $(x(\mu), \mu)$. 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 $\mu$ varies. Thus, in particular, there must exist a value $\mu_0$ with steady state $x(\mu_0)$, where the Jacobian $J(x(\mu_0), \mu_0)$ has exactly one pair of purely imaginary eigenvalues. More generally, a Hopf bifurcation necessitates $J(x(\mu_0), \mu_0)$, for some value $\mu_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} + \cdots + a_s$ is the following $i \times i$ matrix:

in which the $(k, l)$-th entry is $a_{2k-l}$ as long as $0\leq 2 k - l \leq s$, and $0$ otherwise. Note $H_s = a_s H_{s-1}$.

Returning to the ODE system $\dot{x} = g_\mu(x)$, consider the characteristic polynomial of the Jacobian matrix $J(x(\mu), \mu)$

where the determinant of a matrix $A$ is denoted by $\det A$. For $i = 1, \dots, s$, we let $H_i(\mu)$ be the $i$-th Hurwitz matrix of $p_{\mu}(z)$.

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\geq 2$, $\mu_0$ be fixed, and let $H_i(\mu_0)$ denote the $i$-th Hurwitz matrix of $p_{\mu_0}(z)$. We assume that

Then $p_{\mu_0}(z)$ has at most one pair of symmetric roots, and has exactly one if and only if $\det H_{s-1}(\mu_0) = 0$. In this case, the pair consists of purely imaginary roots if and only if $a_s(\mu_0) > 0$.

Proof. The first statement follows from Corollary 3.2 in ^[17] using (4). So, assume $p_{\mu_0}(z)$ has a pair of symmetric roots $z, -z$ such that $\det H_{s-1}(\mu_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_{\mu_0}(z)$ have negative real part. Now $a_s(\mu_0)$ is $(-1)^s$ times the product of the roots of $p_{\mu_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_{\mu_0}(z)$ is real if $a_s(\mu_0)\leq 0$ and purely imaginary if $a_s(\mu_0) > 0$.

As a consequence of Propostion 1, we obtain the following criterion to preclude Hopf bifurcations:

Theorem 1. Consider the dynamical system $\dot x = g(x)$ with $s\geq 2$ and assume there exists a curve of steady states $x(\mu)$ as in the above setting. As before, let $p_{\mu}(z)$ be the characteristic polynomial of the Jacobian $J(x(\mu), \mu)$ of $g_{\mu}(x)$ evaluated at $x(\mu)$. Further let $H_i(\mu)$ denote the $i$-th Hurwitz matrix of $p_{\mu}(z)$.

If

for all $\mu$ and either

or

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+\text{im} N$. This implies that $0$ is a root of the characteristic polynomial of $N v({\kappa}, x)$ with multiplicity $n-s$ and hence it factors as

The polynomial $a_{0}({\kappa}, x) z^s + a_{1}({\kappa}, x) z^{s-1} + \cdots + a_s({\kappa}, x)$ is the characteristic polynomial of the Jacobian of the restriction of system (2) to $x+\text{im} N$, and hence we apply Theorem 1 to this polynomial.

3.2.   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 ${\kappa}$ 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^*\in \mathbb{R}^n_{ > 0}$ of (2) satisfy the polynomial system $N~{\rm{diag}} ({\kappa}) \psi (x^*) = 0$. Then the rate functions $v = v({\kappa}, 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 $\{E_1, \ldots, E_l \}$

where $E$ is the matrix with columns $E_1, \dots, E_l$ and $\lambda = (\lambda_1, \dots, \lambda_l)$.

Remark 2. (cf. ^[22]). (a) The vectors $E_1, E_2, \ldots, E_l$ need not be linearly independent.

(b) If all extreme vectors $E_1, \ldots, E_l$ are unit vectors, then their choice is unique.

(c) When $v = v({\kappa}, x^*)$ and $x^*\in \mathbb{R}^n_{ > 0}$, then all components of $v$ are positive. This might impose some restrictions on the possible values of $\lambda$ in (6).

The nonnegative coefficients $\lambda_1, \ldots, \lambda_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, \lambda) = (h_1, \dots, h_n, \lambda_1, \ldots, \lambda_l)\in \mathbb{R}^n_{ > 0} \times \mathbb{R}^l_{\geq 0}$ such that $E\lambda\in \mathbb{R}^r_{ > 0}$.

The convex parameters are convenient for parameterizing the Jacobian $J({\kappa}, x)$ evaluated at a positive steady state $x = x^*$. To see this, note that the $(j, i)$-th entry of the Jacobian of $v({\kappa}, x)$ evaluated at $x^*$ is

Hence, the Jacobian of $N v({\kappa}, x)$ evaluated at $x^*$ is

where we use the vector notation

Using now (6) and (7) to write $v({\kappa}, x^*) = E\lambda$, the Jacobian of $Nv({\kappa}, x)$ evaluated at $x^*$ can be written as

with $(h, \lambda)$ a vector of convex parameters. Therefore, given a vector of rate constants ${\kappa}$ and a corresponding positive steady state $x^*$, there exist convex parameters $(h, \lambda)$ such that equality (8) holds.

Conversely, given convex parameters $(h, \lambda)$, we define $x^* = 1/h$ and let

which is a positive vector since all entries of $E\lambda$ are positive. Then, using that $\psi_j(x)^{-1} = \psi_j(x^{-1})$, we obtain $v({\kappa}, x^*) = {\rm{diag}}({\kappa}) \psi(x^*) = E\lambda$, and equality (8) holds as well. This proves the following proposition:

Proposition 2. The set of Jacobian matrices $J({\kappa}, x^*)$ for all ${\kappa}$ and corresponding positive steady states $x^*$ agrees with the set of matrices defined by the right-hand side of (8), for all $h\in \mathbb{R}^n_{ > 0}$ and $\lambda\in \mathbb{R}^l_{\geq 0}$ such that $E\lambda\in \mathbb{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$, $\lambda$) 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 ${\kappa}$ and $x^*$ where a Hopf bifurcation occurs.

Remark 3. The coefficients of the characteristic polynomial of $J({\kappa}, x^*)$ are homogeneous polynomials in the convex parameters if the Jacobian matrix is parametrized as in (8).

4.   The mass action model derived from Figure 1a and its simplifications

Figure 1a contains four phosphorylation events: the phosphorylation of $S_0$ and $S_1$ and the dephosphorylation of $S_2$ and $S_1$. 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 9. 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 ($\mathcal{N}$):

(ⅰ) 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.

(ⅱ) 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 $\mathcal{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$):

● The network containing only $KS_1$ and $FS_1$:

Remark 4. The dynamics of the network obtained upon removal of an intermediate resembles that of the original network when the binding reaction occurs at a slower time scale than unbinding reactions. 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 ($\mathcal{N}_1$) – ($\mathcal{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_{\text{MM}}$ is the familiar Michaelis-Menten kinetics

with catalytic constant $k_{cat}$, total enzyme concentration $K_0$, Michaelis-Menten constant $\mathit{K}_m$, and where $[\cdot]$ 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 9 does not exhibit periodic orbits.

A particular feature of the Michaelis-Menten kinetic $v_{\text{MM}}$ is that it is approximately linear in $[S_i]$ for small values of $[S_i]$ and approximately constant for very large values. Thus the simplification based on the Michaelis-Menten approximation covers, for example, the saturation of available enzyme with substrate very well. However, it does not capture the influence of varying concentrations of the free enzyme $K$.

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

Clearly, $v_{\text{MA}}$ is linear in $[S_i]$, hence for small values of $S_i$ and values of $[K]$ close to its total concentration $K_0$ the two reaction rates $v_{\text{MA}}$ and $v_{\text{MM}}$ behave similar. However, $v_{\text{MA}}$ is not able to reproduce the saturation of enzyme by the substrate. But it does capture the influence of varying concentrations of free enzyme $[K]$. Hence for some values of $[K]$ and $[S_i]$ our simplification behaves similar to the one based on the Michaelis-Menten approximation, but not everywhere. Moreover, our simplification covers the influence of varying concentrations of $K$.

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 ($\mathcal{N}_1$) – ($\mathcal{N}_4$) had a Hopf bifurcation giving rise to periodic solutions, then by ^[31], so would the full mechanism in 9. 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 $K S_0$, $K S_1$, $F S_1$ or $F S_2$.

5.   Absence of Hopf bifurcations

In this section we apply Theorem 1, using the discussion after it, to the networks ($\mathcal{N}_1$) – ($\mathcal{N}_4$). To this end we use Eq (8) to determine the Jacobian matrices $J_1(h, \lambda)$, $\ldots$, $J_4(h, \lambda)$ of networks ($\mathcal{N}_1$) – ($\mathcal{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:

Remark 7. (ⅰ) Every network ($\mathcal{N}_1$) – ($\mathcal{N}_4$) consists of $7$ species and $6$ reactions.

(ⅱ) The stoichiometric matrix of every network has rank $s = 4$, as every network has three conserved quantities (the total amount of substrate, kinase and phosphatase).

(ⅲ) For every network the cone (5) is spanned by two nonnegative vectors $w_0$ and $w_1$ such that $\lambda_1 w_0 + \lambda_2 w_1$ is a positive vector if and only if $\lambda_1, \lambda_2 > 0$.

Observe that a scaling $\alpha \lambda$ of the vector $\lambda = (\lambda_1, \lambda_2)$ with $\alpha > 0$ translates into a scaling $\alpha {\kappa}$ of the vector of rate constants ${\kappa}$ under the correspondence of parameterization in Proposition 2. Further $x^*$ is a steady state for ${\kappa}$ if and only if it is for $\alpha {\kappa}$ and $\alpha J({\kappa}, x^*) = J(\alpha {\kappa}, x^*)$. The latter implies that any eigenvalue of $J(\alpha {\kappa}, x^*)$ is in fact an eigenvalue of $J({\kappa}, x^*)$ multiplied by $\alpha > 0$. Thus, $J(\alpha {\kappa}, x^*)$ has a pair of purely imaginary eigenvalues if and only if $J({\kappa}, x^*)$ has a pair of purely imaginary eigenvalues.

Hence, it is enough to take one of $\lambda_1, \lambda_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, \lambda)$ according to Eq (8) parametrized by $8$ parameters: the parameter $\lambda$ and $h_1, \dots, h_7$.

By Remark 7 (ⅰ) and (ⅱ) it follows that the characteristic polynomial of every Jacobian $J_i(h, \lambda)$ is a degree $7$ polynomial of the form

where each $a_i$ depends on the $8$ parameters $\lambda, h_1, \ldots, h_7$ (cf. the file 'SupplMat1.mw'). Following the discussion after Theorem 1, for each network we compute $a_0(h, \lambda), \ldots, a_4(h, \lambda)$, $\det H_1(h, \lambda)$, $\det H_2(h, \lambda)$ and $\det H_3(h, \lambda)$ (cf. Definition 1) and show the following proposition:

Proposition 3. With the notation above, we have

(i) Network ($\mathcal{N}_1$): $\det H_1(h, \lambda)$ and $\det H_2(h, \lambda)$ contain only positive monomials, $a_4(h, \lambda)$ and $\det H_3(h, \lambda)$ contain monomials of both signs. But $\det H_3(h, \lambda)$ is positive whenever $a_4(h, \lambda) > 0$ (Proposition 4 in Section 5.1).

(ii) Network ($\mathcal{N}_2$) and ($\mathcal{N}_3$): $\det H_1(h, \lambda)$, $\det H_2(h, \lambda)$ and $\det H_3(h, \lambda)$ contain only positive monomials, $a_4(h, \lambda)$ contains monomials of both signs. Thus $\det H_i(h, \lambda) > 0$, $i = 1, 2, 3$ for positive $h$ and $\lambda$ and in particular, $\det H_3(h, \lambda)\neq 0$.

(iii) Network ($\mathcal{N}_4$): $\det H_1(h, \lambda)$, $\det H_2(h, \lambda)$ and $\det H_3(h, \lambda)$ and $a_4(h, \lambda)$ contain only positive monomials. As in case of $\mathcal{N}_2$ and $\mathcal{N}_3$, one has $\det H_i(h, \lambda) > 0$, $i = 1, 2, 3$ for positive $h$ and $\lambda$ and in particular, $\det H_3(h, \lambda)\neq 0$.

In particular, this proposition (which is proven below) tells us that in all four networks, $\det H_3(h, \lambda)\neq 0$ whenever $a_4(h, \lambda) > 0$. As a consequence we obtain the following theorem.

Theorem 2. For the networks ($\mathcal{N}_1$) – ($\mathcal{N}_4$) there do not exist rate constants ${\kappa}$ and a corresponding positive steady state $x^*$ such that the Jacobian matrix $J_i({\kappa}, x^*)$ has a pair of purely imaginary eigenvalues. Thus, in particular, there do not exist ${\kappa}$ 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 ${\kappa}$ and a corresponding positive steady state $x^*$ such that the corresponding Hurwitz determinant $\det H_3({\kappa}, x^*)$ vanishes and the coefficient of lowest degree of the characteristic polynomial $a_4({\kappa}, x^*)$ is positive. Further, for all networks $\det H_1({\kappa}, x^*) > 0$ and $\det H_2({\kappa}, x^*) > 0$ for all ${\kappa}, x^*$. Hence, by Theorem 3.1, the Jacobian matrix does not have a pair of purely imaginary eigenvalues.

Remark 8. By ^[33], the networks $\mathcal{N}_1$, $\mathcal{N}_2$ and $\mathcal{N}_3$ are multistationary, while $\mathcal{N}_4$ is not. Hence, by ^[32], the coefficient of lowest degree of the Jacobian must vanish for some ${\kappa}$ and positive steady state $x^*$, and consistently $a_4(h, \lambda)$ must contain monomials of both signs.

Remark 9. We have focused on networks with two intermediates, as network (($\mathcal{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 ($\mathcal{N}$) obtained after removing three intermediates, that is, with only one intermediate left, satisfy that $\det H_{s-1}(h, \lambda)$ 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, \lambda)$. For the networks (($\mathcal{N}_1$)) – (($\mathcal{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, \lambda)$ and the determinants of the Hurwitz matrices $\det H_i(h, \lambda)$. The computations are given in the supplementary file 'SupplMat1.mw'.

5.1.   Network ($\mathcal{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, \dots, h_7, \lambda$ as given in (8) is

We observe that the coefficient $a_4(h, \lambda)$ of the characteristic polynomial contains monomials of both signs. We compute the associated Hurwitz determinants $\det H_1(h, \lambda)$, $\det H_2(h, \lambda)$ and $\det H_3 (h, \lambda)$ and obtain that $\det H_1(h, \lambda)$ and $\det H_2(h, \lambda)$ are sums of positive monomials and that $\det H_3(h, \lambda)$ contains monomials of both signs as well. Hence both $a_4(h, \lambda)$ and $\det H_3(h, \lambda)$ contain monomials of both signs and can potentially be zero.

In the remainder of this section we prove the following:

Proposition 4. Consider the coefficient $a_4(h, \lambda)$ and the Hurwitz determinant $\det H_3(h, \lambda)$ of network ($\mathcal{N}_1$) and given in the file 'SupplMat1.mw'. Then

Proof. The coefficient $a_4(k, \lambda)$ of the Jacobian is

Since $\lambda^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\geq h_1$ and $h_7\geq h_2$, such that one of the two inequalities is strict, otherwise $c_0 = 0$. In this case $c_0\geq 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\geq h_6$ and $h_2\leq 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\geq h_3$, then again $\det H_3$ is positive. So we restrict to $h_3 > h_5$ and perform the substitution $h_3 = h_5+v_3$ into $\det H_3$.

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

Then any number in the interval $[0, z_1)$ is of this form for some nonnegative $\mu$. 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 $\mu, v_2, v_3$ are zero.

The other case $h_2\geq h_7$ and $h_1\leq 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.

5.2.   Network $\mathcal{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, K S_0, F S_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, \dots, h_7, \lambda)$ is:

$\det H_3(h, \lambda)$ contains only positive monomials and in particular it does not vanish for any positive $h, \lambda$.

5.3.   Network $\mathcal{N}_3$

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, K S_0, K S_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, \dots, h_7, \lambda)$ is:

$\det H_3(h, \lambda)$ contains only positive monomials and in particular it does not vanish for any positive $h, \lambda$.

5.4.   Network ($\mathcal{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, K S_1, F S_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, \dots, h_7, \lambda)$ is:

$\det H_3(h, \lambda)$ contains only positive monomials and in particular it does not vanish for any positive $h, \lambda$.

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).

In ^[16], the authors consider the following irreversible version of the network ($\mathcal{N}$):

In network ${{\mathcal{N}}_{f}}$ the constants ${\kappa}_2$, ${\kappa}_5$, ${\kappa}_8$ and ${\kappa}_{11}$ describing the 'backward' reactions in the full network ($\mathcal{N}$) have been assigned the value zero. Hence the matrix $E_f$ whose columns generate the cone (5) is

Using $x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9$ for the concentrations of $S_0, K, KS_0, S_1, KS_1, S_2, F, FS_2, FS_1$ respectively one obtains the Jacobian matrix

The authors of ^[16] provide the candidate point where $\det H_{5}(h^*, \lambda^*) = 0$:

We find that

While $\det H_4(h^*, \lambda^*)$, $\det H_5(h^*, \lambda^*)$ and $a_6(h^*, \lambda^*)$ are all very small numbers, $\det H_5(h^*, \lambda^*)$ is by orders of magnitude smaller. And if $10^{-42} \approx 0$, then $\det H_4(h^*, \lambda^*)$ and $a_6(h^*, \lambda^*)$ are nonzero. In particular $a_6(h^*, \lambda^*) < 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 $\approx 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 ($\mathcal{N}$) is still open. For the full network, $\det H_1, \dots, \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 ($\mathcal{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 ($\mathcal{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 ($\mathcal{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 ($\mathcal{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_1^d\, p({\kappa}, x)$ becomes $(k_2+a)^d\, p({\kappa}, x)$. We also performed substitutions of the form $k_1 = \tfrac{\mu}{\mu+1} \tfrac{p({\kappa}, x)}{q({\kappa}, x)}$ with $\mu > 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 ($\mathcal{N}_1$), is prohibitive. Nevertheless, we think ideas introduced here in the analysis of network ($\mathcal{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 non-existence 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 ($\mathcal{N}_1$)-($\mathcal{N}_4$) are preserved after introducing unbinding reactions with a small rate constant.

Acknowledgements

CC acknowledges funding from Deutsche Forschungsgemeinschaft, 284057449. EF acknowledges funding from the Danish Research Council for Independent Research.

The authors thank Angélica Torres for providing source code to compute Hurwitz determinants and Alan Rendall for helpful discussions about (simple) Hopf bifurcations. The authors thank the hospitality Erwin Schrödinger International Institute for Mathematics and Physics in Vienna, where this project was started during the workshop 'Advances in Chemical Reaction Network Theory'.

Conflict of interest

The authors declare there is no conflicts of interest.

Supplementary

The Maple computations for the proof of Theorem 5.1 are provided in the file 'SupplMat1.mw'. For the convenience of the readers without access to Maple, we provide as well as pdf version of the file.