Oscillations and bistability in a model of ERK regulation

TL;DR

This paper analyzes how oscillations and bistability arise in a model of ERK signaling, showing oscillations are robust while bistability depends on specific reversible reactions, with implications for understanding cellular signaling dynamics.

Contribution

The study provides a detailed mathematical analysis of ERK network dynamics, revealing conditions for oscillations and bistability, and how these properties depend on reaction reversibility.

Findings

Oscillations persist even with simplified, irreversible reactions.

Bistability is lost when intermediates are removed or reactions are made irreversible.

Reversible catalytic reactions are crucial for bistability.

Abstract

This work concerns the question of how two important dynamical properties, oscillations and bistability, emerge in an important biological signaling network. Specifically, we consider a model for dual-site phosphorylation and dephosphorylation of extracellular signal-regulated kinase (ERK). We prove that oscillations persist even as the model is greatly simplified (reactions are made irreversible and intermediates are removed). Bistability, however, is much less robust -- this property is lost when intermediates are removed or even when all reactions are made irreversible. Moreover, bistability is characterized by the presence of two reversible, catalytic reactions: as other reactions are made irreversible, bistability persists as long as one or both of the specified reactions is preserved. Finally, we investigate the maximum number of steady states, aided by a network's "mixed volume" (a concept from convex geometry). Taken together, our results shed light on the question of how oscillations and bistability emerge from a limiting network of the ERK network -- namely, the fully processive dual-site network -- which is known to be globally stable and therefore lack both oscillations and bistability. Our proofs are enabled by a Hopf bifurcation criterion due to Yang, analyses of Newton polytopes arising from Hurwitz determinants, and recent characterizations of multistationarity for networks having a steady-state parametrization.

Tables

Table 1. Table 1 : Summary of results. Yes* indicates that the fully irreversible ERK network exhibits oscillations (see Figure  3 ), and 5* indicates that 5 is the maximum number of complex-number steady states for the network obtained from the full ERK network by setting k on = 0 subscript 𝑘 on 0 k_{\mathrm{on}}=0 . For details on results, see Propositions  4.1 , 4.5 , and  5.9 , and Theorem  4.6 .

ERK			Maximum #	Maximum #
network	Oscillations?	Bistability?	steady states	over $ℂ$
Full	Yes [40]	Yes [40]	Conjecture: 3	7
Irreversible	Yes*	If $k_{\mathrm{on}}>0$ or ${\ell }_{\mathrm{on}}>0$	1	5*
Reduced	Yes	No	1	3

Table 2. Table 2 : Assignment of variables to species for the ERK network in Figure  1 .

$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$	$x_8$	$x_9$	$x_{10}$	$x_{11}$	$x_{12}$
$S_{00}$	$E$	$F$	$S_{11}F$	$S_{10}F$	$S_{01}F$	$S_{01}E$	$S_{10}E$	$S_{01}$	$S_{10}$	$S_{00}E$	$S_{11}$

Table 3. Table 3 : Assignment of variables to species for the reduced ERK network in Figure  2 . (Many of the variables that are also in the full ERK, in Table  2 , have been relabeled.)

$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$x_7$	$x_8$	$x_9$	$x_{10}$
$S_{00}$	$E$	$S_{00}E$	$S_{01}E$	$S_{11}$	$S_{01}$	$S_{10}$	$F$	$S_{11}F$	$S_{10}F$

Table 4. Table 4 : Results on ERK networks.

ERK	Max #	Max #	Mixed
network	positive steady states	over ${ℂ}^{\ast }$	volume
Full	$\ge 3$	7	7
Full with $k_{\mathrm{on}}=0$	$\ge 3$	5	5
Fully irreversible	1	3	3
Reduced	1	3	3

Table 5. Table 5 : Supporting Information files and the results they support.

Name	File type	Result
ERK-Matcont.txt	text file with MATCONT instructions	Figures 3 and 4
irreversibleERK.mw	Maple	Theorem 4.6
reducedERK-noMSS.mw	Maple	Proposition 4.5
reducedERK-hopf.mw	Maple	Theorem 4.3
reducedERK-cones.sws	Sage	Theorem 4.3
ERK-mixedVol.m2	PHCPack	Proposition 5.9
ERK-MaxComplexNumber.nb	Mathematica	Proposition 5.9

Equations

$${\alpha }_{1j}{X}_1+{\alpha }_{2j}{X}_2+\cdots +{\alpha }_{sj}{X}_s\to {\beta }_{1j}{X}_1+{\beta }_{2j}{X}_2+\cdots +{\beta }_{sj}{X}_s,\mathrm{for}j=1,2,\dots ,m,$$

$$\dot{x}=f(x):=N\cdot \left(\begin{array}{c}{\kappa }_1{x}_1^{{\alpha }_{11}}{x}_2^{{\alpha }_{21}}\cdots x_s^{{\alpha }_{s1}}\\ {\kappa }_2{x}_1^{{\alpha }_{12}}{x}_2^{{\alpha }_{22}}\cdots x_s^{{\alpha }_{s2}}\\ ⋮\\ {\kappa }_m{x}_1^{{\alpha }_{1m}}{x}_2^{{\alpha }_{2m}}\cdots x_s^{{\alpha }_{sm}}\end{array}\right),$$

$${\mathcal{S}}_c:=\{x\in {\mathbb{R}}_{\ge 0}^s\mid Wx=c\}.$$

$$f_{c,\kappa ,i}=f_{c,\kappa }(x{)}_i:=\{\begin{array}{ll}{f}_i(x)& \text{if}i\notin I,\\ (Wx-c{)}_k& \text{if}i=i_k\in I.\end{array}$$

$$\bar{a}:{\mathbb{R}}_{>0}^m$$

$$\kappa$$

$$({\bar{h}}_{j_1},{\bar{h}}_{j_2},\dots ,{\bar{h}}_{j_{s-d}}{)}^{\top }:=M(\kappa )(f_{j_1},f_{j_2},\dots ,f_{j_{s-d}}{)}^{\top },$$

$${\bar{h}}_{j_{\ell }}=h_{j_{\ell }}{|}_{a_1={\bar{a}}_1(\kappa ),\dots ,a_{\bar{m}}={\bar{a}}_{\bar{m}}(\kappa )}.$$

$$h_{c,a}(x{)}_i:=\{\begin{array}{ll}{h}_i(a;x)& \text{if}i\notin I\\ (Wx-c{)}_k& \text{if}i=i_k\in I.\end{array}$$

$$C(\hat{a};x)=\left(\det \mathrm{Jac}{h}_{c,a}\right){|}_{(a;x)=\varphi (\hat{a};x)},$$

$$\mathrm{sign}(C({\hat{a}}^{\ast };x^{\ast }))=(-1{)}^{\mathrm{rank}(N)+1}.$$

$$\mathrm{sign}(C(\hat{a};x))=(-1{)}^{\mathrm{rank}(N)}.$$

$$H_i=\left(\begin{array}{ccccccc}{b}_1& b_0& 0& 0& 0& \cdots & 0\\ b_3& b_2& b_1& b_0& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& & ⋮\\ b_{2i-1}& b_{2i-2}& b_{2i-3}& b_{2i-4}& b_{2i-5}& \cdots & b_i\end{array}\right),$$

$$\dot{x}=g_{\mu }(x),$$

$$\mu \mapsto x(\mu )$$

$$p_{\mu }(\lambda ):=\det \left(\lambda I-\mathrm{Jac}{g}_{\mu }\right){|}_{x=x(\mu )}={\lambda }^n+b_1(\mu ){\lambda }^{n-1}+\cdots +b_n(\mu ),$$

$$q(\lambda ):=\frac{1}{{\lambda }^d}\det \left(\lambda I-\mathrm{Jac}f\right){|}_{(\kappa ;x)=\varphi (\hat{\kappa };x)}.$$

$$\frac{1}{{\lambda }^d}\det (\lambda I-\mathrm{Jac}f){|}_{x=x^{\ast },{\kappa }_i={\kappa }_i^{\ast }foralli\ne j}.$$

$$\begin{array}{rl}\dot{x_1}& =-k_1{x}_1{x}_2+k_2{x}_{11}+{\ell }_{\mathrm{cat}}{x}_5+n_3{x}_6\\ \dot{x_2}& =-k_1{x}_1{x}_2-k_{\mathrm{on}}{x}_2{x}_9-m_2{x}_2{x}_{10}+k_2{x}_{11}+k_{\mathrm{cat}}{x}_7+k_{\mathrm{off}}{x}_7+m_1{x}_8+m_3{x}_8\\ \dot{x_3}& =-{\ell }_1{x}_3{x}_{12}-{\ell }_{\mathrm{on}}{x}_3{x}_{10}-n_1{x}_3{x}_9+{\ell }_2{x}_4+{\ell }_{\mathrm{cat}}{x}_5+{\ell }_{\mathrm{off}}{x}_5+n_2{x}_6+n_3{x}_6\\ \dot{x_4}& ={\ell }_1{x}_3{x}_{12}-{\ell }_2{x}_4-{\ell }_3{x}_4\\ \dot{x_5}& ={\ell }_{\mathrm{on}}{x}_3{x}_{10}+{\ell }_3{x}_4-{\ell }_{\mathrm{cat}}{x}_5-{\ell }_{\mathrm{off}}{x}_5\\ \dot{x_6}& =n_1{x}_3{x}_9-n_2{x}_6-n_3{x}_6\\ \dot{x_7}& =k_{\mathrm{on}}{x}_2{x}_9+k_3{x}_{11}-k_{\mathrm{cat}}{x}_7-k_{\mathrm{off}}{x}_7\\ \dot{x_8}& =m_2{x}_2{x}_{10}-m_1{x}_8-m_3{x}_8\\ \dot{x_9}& =-k_{\mathrm{on}}{x}_2{x}_9-n_1{x}_3{x}_9+k_{\mathrm{off}}{x}_7+n_2{x}_6\\ \dot{x_{10}}& =-{\ell }_{\mathrm{on}}{x}_3{x}_{10}-m_2{x}_2{x}_{10}+{\ell }_{\mathrm{off}}{x}_5+m_1{x}_8\\ \dot{x_{11}}& =k_1{x}_1{x}_2-k_2{x}_{11}-k_3{x}_{11}\\ \dot{x_{12}}& =-{\ell }_1{x}_3{x}_{12}+k_{\mathrm{cat}}{x}_7+{\ell }_2{x}_4+m_3{x}_8\end{array}$$

$$\begin{array}{rl}{x}_1+x_4+x_5+x_6+x_7+x_8+x_9+x_{10}+x_{11}+x_{12}& ={\text{S}}_{\text{tot}}=:c_1\\ x_2+x_7+x_8+x_{11}& ={\text{E}}_{\text{tot}}=:c_2\\ x_3+x_4+x_5+x_6& ={\text{F}}_{\text{tot}}=:c_3.\end{array}$$

$$h_{c,a,1}$$

$$h_{c,a,2}$$

$$h_{c,a,3}$$

$$h_{c,a,4}$$

$$h_{c,a,5}$$

$$h_{c,a,6}$$

$$h_{c,a,7}$$

$$h_{c,a,8}$$

$$h_{c,a,9}$$

$$h_{c,a,10}$$

$$h_{c,a,11}$$

$$h_{c,a,12}$$

$$\varphi :{\mathbb{R}}_{>0}^{2+\text{mathds}{1}_{k_{\mathrm{on}}}+\text{mathds}{1}_{{\ell }_{\mathrm{on}}}+12}$$

$$(\hat{a};x_1,x_2,\dots ,x_{12})$$

$$a_1$$

$$a_7$$

$$a_{11}$$

$${\bar{a}}_1=\frac{{\ell }_{cat}}{k_{\mathrm{cat}}}$$

Full text

Keywords: chemical reaction network, Hopf bifurcation, oscillation, bistable, Newton polytope, mixed volume

1 Introduction

In recent years, significant attention has been devoted to the question of how bistability and oscillations emerge in biological networks involving multisite phosphorylation [14]. Such networks are of great biological importance [9]. The one we consider is the network, depicted in Figure 1, comprising extracellular signal-regulated kinase (ERK) regulation by dual-site phosphorylation by the kinase MEK (denoted by $E$) and dephosphorylation by the phosphatase MKP3 ($F$) [40]. This network, which we call the ERK network, has an important role in regulating many cellular activities, with dysregulation implicated in many cancers [42]. Accordingly, an important problem is to understand the dynamical properties of the ERK network, with the goal of predicting effects arising from mutations or drug treatments [21].

The ERK network was shown by Rubinstein, Mattingly, Berezhkovskii, and Shvartsman [40] to be bistable and exhibit oscillations (for some choices of rate constants). Rubinstein et al. also observed that the ERK network “limits” to a network without bistability or oscillations. Namely, when the rate constants $k_{\mathrm{cat}}$ and ${\ell}_{\mathrm{cat}}$ are much larger than $k_{\mathrm{off}}$ and ${\ell}_{\mathrm{off}}$, respectively, this yields the “fully processive” network obtained by deleting all vertical arrows in Figure 1, which is globally convergent to a unique steady state [13, 18]. Accordingly, Rubinstein et al. asked, How do bistability and oscillations in the ERK network emerge from the processive limit? This question was subsequently articulated as follows [14]:

Question 1.1**.**

When the processivity levels $p_{k}:={k_{\mathrm{cat}}}/{(k_{\mathrm{cat}}+k_{\mathrm{off}})}$ and $p_{\ell}:={{\ell}_{\mathrm{cat}}}/{({\ell}_{\mathrm{cat}}+{\ell}_{\mathrm{off}})}$ are arbitrarily close to 1, is the ERK network still bistable and oscillatory?

One of our main contributions is to lay foundation toward answering Question 1.1. Specifically, we answer a related question, How do bistability and oscillations emerge from simpler versions of the ERK network? Our main results, summarized in Table 1, are that oscillations are surprisingly robust to operations that simplify the network, while bistability is lost more easily. Specifically, oscillations persist even as reactions are made irreversible and intermediates are removed (see Section 4.1), while bistability is lost more quickly, when only a few reactions are made irreversible (Section 4.2). Taken together, our results form a case study for the problem of model choice – an investigation into the simplifications of a model that preserve important dynamical properties.

Our focus here – on determining which operations on the ERK network preserve oscillations and bistability – is similar in spirit to the recent approach of Sadeghimanesh and Feliu [41]. Indeed, there has been significant interest in understanding which operations on networks preserve oscillations [3], bistability [4, 20, 32], and other properties [24].

A related topic – mentioned earlier – is the question of how dynamical properties arise in phosphorylation systems. Several works have examined this problem at the level of parameters, focusing on the question of which rate constants and/or initial conditions give rise to oscillations [12] or bistability [10, 11]. Our perspective is slightly different; instead of allowing parameter values to change, we modify the network itself. Accordingly, our work is similar in spirit to recent investigations into minimal oscillatory or bistable networks [3, 4, 27, 28, 34].

A key tool we use is a parametrization of the steady states. Such parametrizations have been shown in recent years to be indispensable for analyzing multistationarity (multiple steady states, which are necessary for bistability) and oscillations [22, 31, 45]. Indeed, here we build on results in [10, 12, 17].

Specifically, following [12], we investigate oscillations by employing a steady-state parametrization together with a criterion of Yang [47] that characterizes Hopf bifurcations in terms of determinants of Hurwitz matrices. In [12], this approach showed that the Hopf bifurcations of a mixed-mechanism phosphorylation network lie on a hypersurface defined by the vanishing of a single Hurwitz determinant. For our ERK networks, however, the problem does not reduce to the analysis of a single polynomial, and the size of these polynomials makes the system difficult to solve. To this end, we introduce an algorithm for analyzing these polynomials, through their Newton polytopes, by using techniques from polyhedral geometry. Using this algorithm, we succeed in finding, for the reduced ERK network, a Hopf bifurcation giving rise to oscillations.

Finally, we investigate the precise number of steady states in ERK networks. For general networks, much has been done for determining which networks admit multiple steady states – see e.g. [10, 15, 17, 19, 22, 33, 37] – but there are few techniques for determining a network’s maximum number of steady states. To this end, we introduce two related measures of a network, the maximum number of complex-number steady states and the “mixed volume”. In general, the mixed volume is an upper bound on the number of complex-number steady states, but we show that these numbers are equal for ERK networks (Section 5).

The outline of our work is as follows. Section 2 contains background on chemical reaction systems, steady-state parametrizations, and Hopf bifurcations. We present steady-state parametrizations for the ERK network and the reduced ERK network in Section 3. Section 4 contains our main results on oscillations and bistability. Section 5 investigates the number of steady states and the relationship to mixed volumes. We end with a Discussion in Section 6.

2 Background

Here we introduce chemical reaction systems (Section 2.1), their steady-state parametrizations (Section 2.2), and Hopf bifurcations (Section 2.3).

2.1 Chemical reaction systems

As in [17], our notation closely matches that of Conradi, Feliu, Mincheva, and Wiuf [10]. A reaction network $G$ (or network for short) comprises a set of $s$ species $\{X_{1},X_{2},\ldots,X_{s}\}$ and a set of $m$ reactions:

[TABLE]

where each $\alpha_{ij}$ and $\beta_{ij}$ is a non-negative integer. The stoichiometric matrix of $G$, denoted by $N$, is the $s\times m$ matrix with $(i,j)$-entry equal to $\beta_{ij}-\alpha_{ij}$. Let $d=s-{\rm rank}(N)$. The stoichiometric subspace, denoted by $S$, is the image of $N$. A conservation-law matrix of $G$, denoted by $W$, is a row-reduced $d\times s$-matrix whose rows form a basis of the orthogonal complement of $S$. If there exists a choice of $W$ for which every entry is nonnegative and each column contains at least one nonzero entry (equivalently, each species occurs in at least one nonnegative conservation law), then $G$ is conservative.

We denote the concentrations of the species $X_{1},X_{2},\ldots,X_{s}$ by $x_{1},x_{2},\ldots,x_{s}$, respectively. These concentrations, under the assumption of mass-action kinetics, evolve in time according to the following system of ODEs:

[TABLE]

where $x=(x_{1},x_{2},\ldots,x_{s})$, and each $\kappa_{j}\in\mathbb{R}_{>0}$ is called a reaction rate constant. By considering the rate constants as a vector of parameters $\kappa=(\kappa_{1},\kappa_{2},\dots,\kappa_{m})$, we have polynomials $f_{\kappa,i}\in\mathbb{Q}[\kappa,x]$, for $i=1,2,\dots,s$. For ease of notation, we often write $f_{i}$ rather than $f_{\kappa,i}$.

A trajectory $x(t)$ beginning at a nonnegative vector $x(0)=x^{0}\in\mathbb{R}^{s}_{\geq 0}$ remains, for all positive time, in the following stoichiometric compatibility class with respect to the total-constant vector $c\coloneqq Wx^{0}\in{\mathbb{R}}^{d}$:

[TABLE]

A steady state of (1) is a nonnegative concentration vector $x^{*}\in\mathbb{R}_{\geq 0}^{s}$ at which the right-hand sides of the ODEs (1) vanish: $f(x^{*})=0$. We distinguish between positive steady states $x^{*}\in\mathbb{R}^{s}_{>0}$ and boundary steady states $x^{*}\in{\mathbb{R}}_{\geq 0}^{s}\backslash{\mathbb{R}}_{>0}^{s}$. Also, a steady state $x^{*}$ is nondegenerate if ${\rm Im}\left({\rm Jac}(f)(x^{*})|_{S}\right)$ is the stoichiometric subspace ${S}$. (Here, ${\rm Jac}(f)(x^{*})$ is the Jacobian matrix of $f$, with respect to $x$, at $x^{*}$.) A nondegenerate steady state is exponentially stable if each of the $\sigma\coloneqq\dim({S})$ nonzero eigenvalues of ${\rm Jac}(f)(x^{*})$ has negative real part.

A network $G$ is multistationary (respectively, bistable) if, for some choice of positive rate-constant vector $\kappa\in\mathbb{R}^{m}_{>0}$, there exists a stoichiometric compatibility class (2) with two or more positive steady states (respectively, exponentially stable positive steady states) of (1). A network is monostationary111Some authors define monostationary to be non-multistationary; the two definitions are equivalent for the ERK networks in this work. if, for every choice of positive rate constants, there is exactly one positive steady state in every stoichiometric compatibility class.

To analyze steady states within a stoichiometric compatibility class, we will use conservation laws in place of linearly dependent steady-state equations, as follows. Let $I=\{i_{1}<i_{2}<\dots<i_{d}\}$ denote the indices of the first nonzero coordinate of the rows of conservation-law matrix $W$. Consider the function $f_{c,\kappa}:{\mathbb{R}}_{\geq 0}^{s}\rightarrow{\mathbb{R}}^{s}$ defined by

[TABLE]

We call system (3), the system augmented by conservation laws. By construction, positive roots of the system of polynomial equations $f_{c,\kappa}=0$ are precisely the positive steady states of (1) in the stoichiometric compatibility class (2) defined by the total-constant vector $c$.

2.2 Steady-state parametrizations

Here we introduce steady-state parametrizations (Definition 2.2) and recall from [17] how to use them to determine whether a network is multistationary (Proposition 2.4). Later we will see how to use parametrizations to detect Hopf bifurcations (Proposition 4.1).

Definition 2.1**.**

Let $G$ be a network with $m$ reactions and $s$ species, and let $\dot{x}=f(x)$ denote the resulting mass-action system. Denote by $W$ a $d\times s$ row-reduced conservation-law matrix and by $I$ the set of indices of the first nonzero coordinates of its rows. Enumerate the complement of $I$ as follows: $[s]\setminus I=\{j_{1}<j_{2}<\dots<j_{s-d}\}$. A set of effective parameters for $G$ is formed by polynomials $\bar{a}_{1}(\kappa),\bar{a}_{2}(\kappa),\dots,\bar{a}_{\bar{m}}(\kappa)\in\mathbb{Q}(\kappa)$ for which the following hold:

1. (i)

$\bar{a}_{i}(\kappa^{*})$ is defined and, moreover, $\bar{a}_{i}(\kappa^{*})>0$ for every $i=1,2,\dots,\bar{m}$ and for all $\kappa^{*}\in\mathbb{R}^{m}_{>0}$, 2. (ii)

the reparametrization map below is surjective:

[TABLE] 3. (iii)

there exists an $(s-d)\times(s-d)$ matrix $M(\kappa)$ with entries in $\mathbb{Q}(\kappa):=\mathbb{Q}(\kappa_{1},\kappa_{2},\dots,\kappa_{m})$ such that:

1. (a)

for all $\kappa^{*}\in\mathbb{R}^{m}_{>0}$, the matrix $M(\kappa^{*})$ is defined and, moreover, $\det M(\kappa^{*})>0$, and 2. (b)

letting $(\bar{h}_{j_{\ell}})$ denote the functions obtained from $(f_{j_{\ell}})$ as follows:

[TABLE]

every nonconstant coefficient in every $\bar{h}_{j_{\ell}}$ is equal to a rational-number multiple of some $\bar{a}_{i}(\kappa)$.

Given such a set of effective parameters, we consider for $\ell=1,2,\dots,s-d,$ polynomials $h_{j_{\ell}}=h_{j_{\ell}}(a;x)\in\mathbb{Q}[a_{1},a_{2},\dots,a_{\bar{m}}][x]$ (here, the $a_{i}$’s are indeterminates) such that:

[TABLE]

For $i=1,2,\dots,s$ and any choice of $c\in\mathbb{R}^{d}_{>0}$ and $a\in\mathbb{R}^{\bar{m}}_{>0}$, set

[TABLE]

We call the function $h_{c,a}:\mathbb{R}^{s}_{>0}\to\mathbb{R}^{s}$ an effective steady-state function of $G$.

The “steady-state parametrizations” that we will use in this work belong to a subclass of the ones introduced in [17]. Thus, for simplicity, Definition 2.2 below is more restrictive than [17, Definition 3.6]. Specifically, our parametrizations have the form $\phi(\hat{a};x)$, while those in [17] are of the form $\phi(\hat{a};\hat{x})$.

Definition 2.2**.**

Let $G$ be a network with $m$ reactions, $s$ species, and conservation-law matrix $W$. Let $f_{c,\kappa}$ arise from $G$ and $W$ as in (3). Suppose that $h_{c,a}$ is an effective steady-state function of $G$, as in (7), arising from a matrix $M(\kappa)$, as in (5), a reparametrization map $\bar{a}$, as in (4), and polynomials $h_{j_{\ell}}$’s as in (6). The positive steady states of $G$ admit a positive parametrization with respect to $h_{c,a}$ if there exists a function $\phi:\mathbb{R}^{\hat{m}}_{>0}\times\mathbb{R}^{s}_{>0}\rightarrow\mathbb{R}_{>0}^{\bar{m}}\times\mathbb{R}_{>0}^{s}$, for some $\hat{m}\leq\bar{m}$, which we denote by $(\hat{a};x)\mapsto\phi(\hat{a};x)$, such that:

1. (i)

$\phi(\hat{a};x)$ extends the vector $(\hat{a};x)$. More precisely, there exists a natural projection $\pi:\mathbb{R}^{\bar{m}}_{>0}\times\mathbb{R}_{>0}^{s}\to\mathbb{R}^{\hat{m}}_{>0}\times\mathbb{R}^{s}_{>0}$ such that $\pi\circ\phi$ is equal to the identity map. 2. (ii)

Consider any $(a;x)\in\mathbb{R}^{\bar{m}}_{>0}\times\mathbb{R}_{>0}^{s}$. Then, the equality $h_{i}(a;x)=0$ holds for every $i\notin I$ if and only if there exists $\hat{a}^{*}\in\mathbb{R}^{\hat{m}}_{>0}$ such that $(a;x)=\phi(\hat{a}^{*};x)$.

We call $\phi$ a positive parametrization or a steady-state parametrization.

Definition 2.3**.**

Under the notation and hypotheses of Definition 2.2, assume that the steady states of $G$ admit a positive parametrization with respect to $h_{c,a}$. For such a positive parametrization $\phi$, the critical function $C:\mathbb{R}^{\hat{m}}_{>0}\times\mathbb{R}^{s}_{>0}\rightarrow\mathbb{R}$ is given by:

[TABLE]

where ${\rm Jac}(h_{c,a})$ denotes the Jacobian matrix of $h_{c,a}$ with respect to $x$.

The following result is a specialization222As noted earlier, here we consider parametrizations of the form $\phi(\hat{a};x)$, while [17] allowed those of the form $\phi(\hat{a};\hat{x})$. Also, “conservative” in Proposition 2.4 can be generalized to “dissipative” [17]. of [17, Theorem 3.12]:

Proposition 2.4**.**

Under the notation and hypotheses of Definitions 2.1–2.3, assume also that $G$ is a conservative network without boundary steady states in any compatibility class. Let $N$ denote the stoichiometric matrix of $G$.

1. (A)

Multistationarity.* $G$ is multistationary if there exists $(\hat{a}^{*};x^{*})\in\mathbb{R}^{\hat{m}}_{>0}\times\mathbb{R}^{s}_{>0}$ such that*

[TABLE] 2. (B)

Monostationarity.* $G$ is monostationary if for all $(\hat{a};x)\in\mathbb{R}^{\hat{m}}_{>0}\times\mathbb{R}^{s}_{>0}$,*

[TABLE]

2.3 Hopf bifurcations

A simple Hopf bifurcation is a bifurcation in which a single complex-conjugate pair of eigenvalues of the Jacobian matrix crosses the imaginary axis, while all other eigenvalues remain with negative real parts. Such a bifurcation, if it is supercritical, generates nearby oscillations or periodic orbits [35].

To detect simple Hopf bifurcations, we will use a criterion of Yang that characterizes Hopf bifurcations in terms of Hurwitz-matrix determinants (Proposition 2.6).

Definition 2.5**.**

The $i$-th Hurwitz matrix of a univariate polynomial $p(\lambda)=b_{0}\lambda^{n}+b_{1}\lambda^{n-1}+\cdots+b_{n}$ is the following $i\times i$ matrix:

[TABLE]

in which the $(k,l)$-th entry is $b_{2k-l}$ as long as $n\geq 2k-l\geq 0$, and [math] otherwise.

Consider an ODE system parametrized by $\mu\in\mathbb{R}$:

[TABLE]

where $x\in\mathbb{R}^{n}$, and $g_{\mu}(x)$ varies smoothly in $\mu$ and $x$. Assume that $x^{0}\in\mathbb{R}^{n}$ is a steady state of the system defined by $\mu_{0}$, that is, $g_{\mu_{0}}(x^{0})=0$. Assume, furthermore, that we have a smooth curve of steady states:

[TABLE]

(that is, $g_{\mu}\left(x(\mu)\right)=0$ for all $\mu$) and that $x(\mu_{0})=x^{0}$. Denote the characteristic polynomial of the Jacobian matrix of $g_{\mu}$, evaluated at $x(\mu)$, as follows:

[TABLE]

and, for $i=1,\dots,n$, define $H_{i}(\mu)$ to be the $i$-th Hurwitz matrix of $p_{\mu}(\lambda)$.

Proposition 2.6** (Yang’s criterion [47]).**

Assume the above setup. Then, there is a simple Hopf bifurcation at $x_{0}$ with respect to $\mu$ if and only if the following hold:

1. (i)

$b_{n}(\mu_{0})>0$, 2. (ii)

$\det H_{1}(\mu_{0})>0$, $\det H_{2}(\mu_{0})>0$, …, $\det H_{n-2}(\mu_{0})>0$, and 3. (iii)

$\det H_{n-1}(\mu_{0})=0$* and $\frac{d(\det H_{n-1}(\mu))}{d\mu}|_{\mu=\mu_{0}}\neq 0$.*

2.4 Using parametrizations to detect Hopf bifurcations

Here we prove a new result on how to use steady-state parametrizations to detect Hopf bifurcations (Theorem 2.8). The result, which uses Yang’s criterion, is a straightforward generalization of the approach used in [12]. We include it here to use later in Section 4, and we also expect it to be useful in future work.

Lemma 2.7**.**

Let $G$ be a network with $s$ species, $m$ reactions, and $d$ conservation laws. Denote the ODEs by $\dot{x}=f(x)$, as in (1). Assume that the positive steady states of $G$ admit a positive parametrization $\phi$ with respect to an effective steady-state function for which the reparametrization map (4) is just the identity map. In other words, the effective parameters $\bar{a}_{i}$ are the original rate constants $\kappa_{i}$, and so we write $\phi:\mathbb{R}^{\hat{m}}_{>0}\times\mathbb{R}^{s}_{>0}\rightarrow\mathbb{R}_{>0}^{m}\times\mathbb{R}_{>0}^{s}$ as $(\hat{\kappa};x)\mapsto\phi(\hat{\kappa};x)$. Assume moreover that each coordinate of $\phi_{i}$ is a rational function: $\phi_{i}(\hat{\kappa};x)\in\mathbb{Q}(\hat{\kappa};x)$ for $i=1,2,\dots,\hat{m}+s$. Then the following is a univariate, degree-$(s-d)$ polynomial in $\lambda$, with coefficients in $\mathbb{Q}(\hat{\kappa};x)$:

[TABLE]

Proof.

This result is straightforward from the fact that the characteristic polynomial of ${\rm Jac}(f)$ is a polynomial of degree $s$ and has zero as a root with multiplicity $d$ (because of the $d$ conservation laws). ∎

Theorem 2.8** (Hopf-bifurcation criterion).**

Assume the hypotheses of Lemma 2.7. Let $\mathfrak{h}_{i}$ (for $i=1,2,\dots,s-d$) be the determinant of the $i$-th Hurwitz matrix of $q(\lambda)$ in (9). Let $\kappa_{j}$ be one of the rate constants in the vector $\hat{\kappa}$. Then the following are equivalent:

1. (1)

there exists a rate-constant vector $\kappa^{*}\in\mathbb{R}^{m}_{>0}$ such that the resulting system (1) exhibits a simple Hopf bifurcation with respect to $\kappa_{j}$ at some $x^{*}\in\mathbb{R}^{s}_{>0}$, and 2. (2)

there exist $\hat{\kappa}^{*}\in\mathbb{R}^{\hat{m}}_{>0}$ and $x^{*}\in\mathbb{R}^{s}_{>0}$ such that

1. (i)

the constant term of the polynomial $q(\lambda)$, when evaluated at $(\hat{\kappa};x)=(\hat{\kappa}^{*};x^{*})$, is positive, 2. (ii)

$\mathfrak{h}_{1}(\hat{\kappa}^{*};x^{*})>0$, $\mathfrak{h}_{2}(\hat{\kappa}^{*};x^{*})>0$, …, $\mathfrak{h}_{s-d-2}(\hat{\kappa}^{*};x^{*})>0$ , and 3. (iii)

$\mathfrak{h}_{s-d-1}(\hat{\kappa}^{*};x^{*})=0$* and $\frac{\partial\mathfrak{h}_{s-d-1}}{\partial\kappa_{j}}|_{(\hat{\kappa};x)=(\hat{\kappa}^{*};x^{*})}\neq 0$.*

Moreover, given $\hat{\kappa}^{*}$ and $x^{*}$ as in (2), a simple Hopf bifurcation with respect to $\kappa_{j}$ occurs at $x^{*}$ when the vector of rate constants is taken to be $\kappa^{*}:=\widetilde{\pi}(\phi(\hat{\kappa}^{*};x^{*}))$. Here, $\widetilde{\pi}:\mathbb{R}_{>0}^{m}\times\mathbb{R}_{>0}^{s}\to\mathbb{R}_{>0}^{m}$ is the natural projection.

Proof.

Due to the $d$ conservation laws, we apply Yang’s criterion (Proposition 2.6) to:

[TABLE]

Now our result follows directly from Proposition 2.6 and Definition 2.2. ∎

Remark 2.9**.**

Theorem 2.8 easily generalizes beyond parametrizations of the form $\phi(\hat{\kappa};{x})$ to those of the form $\phi(\hat{\kappa};\hat{x})$ or $\phi(\kappa;\hat{x})$. Indeed, one of the form $\phi(\kappa;\hat{x})$ was used in [12] to establish Hopf bifurcations in a mixed-mechanism phosphorylation system.

3 ERK networks and steady-state parametrizations

Here we introduce steady-state parametrizations for the full ERK network and also irreducible and reduced versions of the network (Propositions 3.1 and 3.3).

3.1 The (full) ERK network

For the full ERK network shown earlier in Figure 1, we let $x_{1},x_{2},\ldots,x_{12}$ denote the concentrations of the species in the order given in Table 2. The resulting ODE system (1) is as follows:

[TABLE]

There are 18 rate constants $k_{i},\ell_{i},m_{i},n_{i}$. The 3 conservation laws correspond to the total amounts of substrate $S$, kinase $E$, and phosphatase $F$, respectively:

[TABLE]

A steady-state parametrization for the full ERK network was given in [17, Examples 3.1 and 3.7]. That parametrization, however, can not specialize to accommodate irreversible versions of the network (in the effective parameters given in [17], two of the denominators are $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$, so we can not set those rate constants to 0). So, in the next subsection, we give an alternate steady-state parametrization that, although quite similar to the one in [17], specializes when considering irreversible versions of the network (see Proposition 3.1).

3.2 Irreversible versions of the ERK network

Here we consider networks obtained from the full ERK network (Figure 1) by making some reversible reactions irreversible. Specifically, we delete one or more of the reactions marked in blue in Figure 1. Our motivation for removing those specific reactions (the ones with rate constants $k_{2},k_{\mathrm{on}},m_{1},\ell_{2},{\ell}_{\mathrm{on}},n_{2}$) rather than any of their opposite reactions is to preserve the main reaction pathways (from $S_{00}$ to $S_{11}$, as well as $S_{10}$ to $S_{11}$, $S_{11}$ to $S_{00}$, and $S_{01}$ to $S_{00}$). At the same time, we do not remove the reactions for $k_{\mathrm{off}}$ or ${\ell}_{\mathrm{off}}$, so that we can still pursue Question 1.1 (which involves $k_{\mathrm{off}}$ and ${\ell}_{\mathrm{off}}$) in a model with fewer reactions. We instead allow the removal of reactions $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$.

Proposition 3.1** (Steady-state parametrization for full and irreversible ERK networks).**

Let $\mathcal{N}$ be the full ERK network or any network obtained from the full ERK network by deleting one or more the reactions corresponding to rate constants $k_{2},k_{\mathrm{on}},m_{1},\ell_{2},{\ell}_{\mathrm{on}},n_{2}$ (marked in blue in Figure 1 ). Let $\mathds{1}_{k_{2}}$ denote the indicator function that is 1 if the reaction labeled by $k_{2}$ is in $\mathcal{N}$ and 0 otherwise; analogously, we also define $\mathds{1}_{k_{\mathrm{on}}}$, $\mathds{1}_{m_{1}}$, $\mathds{1}_{\ell_{2}}$, $\mathds{1}_{{\ell}_{\mathrm{on}}}$, and $\mathds{1}_{n_{2}}$. Then $\mathcal{N}$ admits an effective steady-state function $h_{c,a}:\mathbb{R}^{12}_{>0}\to\mathbb{R}^{12}$ given by:

[TABLE]

Moreover, with respect to this effective steady-state function, the positive steady states of $\mathcal{N}$ admit the following positive parametrization:

[TABLE]

given by

[TABLE]

Here, $\hat{a}=(a_{2},a_{4},a_{6},a_{8})$ if $\mathcal{N}$ contains the reactions labeled by $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$, and $\hat{a}=(a_{2},a_{4},a_{6})$ if $\mathcal{N}$ contains the reaction ${\ell}_{\mathrm{on}}$ but not $k_{\mathrm{on}}$, and so on.

Proof.

We will show that the map $\bar{a}:\mathbb{R}^{12+\mathds{1}_{k_{2}}+\mathds{1}_{k_{\mathrm{on}}}+\mathds{1}_{m_{1}}+\mathds{1}_{\ell_{2}}+\mathds{1}_{{\ell}_{\mathrm{on}}}+\mathds{1}_{n_{2}}}_{>0}\to\mathbb{R}^{11+\mathds{1}_{k_{\mathrm{on}}}+\mathds{1}_{{\ell}_{\mathrm{on}}}}_{>0}$, defined as follows, is a reparametrization map as in (4):

[TABLE]

In particular, we remove the effective parameter $\bar{a}_{6}$ (respectively, $\bar{a}_{8}$) if $\mathds{1}_{{\ell}_{\mathrm{on}}}=0$ (respectively, $\mathds{1}_{k_{\mathrm{on}}}=0$). Notice that each $\bar{a}_{i}$ (if it is not removed) is defined and positive for all $\kappa=(k_{1},\dots,n_{3})\in\mathbb{R}^{12+\mathds{1}_{k_{2}}+\mathds{1}_{k_{\mathrm{on}}}+\mathds{1}_{m_{1}}+\mathds{1}_{\ell_{2}}+\mathds{1}_{{\ell}_{\mathrm{on}}}+\mathds{1}_{n_{2}}}_{>0}$.

We must show that the map $\bar{a}$ is surjective. Indeed, given $a\in\mathbb{R}^{11+\mathds{1}_{k_{\mathrm{on}}}+\mathds{1}_{{\ell}_{\mathrm{on}}}}_{>0}$, it is easy to check that $a$ is the image under $\bar{a}$ of the vector obtained by removing every 0 coordinate from the following vector:

[TABLE]

Next, consider the following $9\times 9$ matrix:

[TABLE]

It is straightforward to check that ${\rm det}M(\kappa)$ is the product of all diagonal terms, and hence is positive for all $\kappa\in\mathbb{R}^{12+\mathds{1}_{k_{2}}+\mathds{1}_{k_{\mathrm{on}}}+\mathds{1}_{m_{1}}+\mathds{1}_{\ell_{2}}+\mathds{1}_{{\ell}_{\mathrm{on}}}+\mathds{1}_{n_{2}}}_{>0}$.

The mass-action ODEs of $\mathcal{N}$ are obtained from those (10) of the full ERK network by replacing the rate constants $k_{2},k_{\mathrm{on}},m_{1},\ell_{2},{\ell}_{\mathrm{on}},n_{2}$, respectively, by $\mathds{1}_{k_{2}}k_{2}$, $\mathds{1}_{k_{\mathrm{on}}}k_{\mathrm{on}}$, $\mathds{1}_{m_{1}}m_{1}$, $\mathds{1}_{\ell_{2}}\ell_{2}$, $\mathds{1}_{{\ell}_{\mathrm{on}}}{\ell}_{\mathrm{on}}$, and $\mathds{1}_{n_{2}}n_{2}$, respectively. To the right-hand sides of these ODEs, we apply the recipe given in equations (5)–(7), using the effective parameters $\bar{a}_{i}$ in (14), the matrix $M(\kappa)$ in (24), and the conservation-law matrix $W$ arising from the conservation laws (11). It is straightforward to check that the result is the function $h_{c,a}(x)$ given in (3.1).

Observe that, for the non-conservation-law equations $h_{c,a,4},\dots,h_{c,a,12}$ in (3.1), each non-constant coefficient is, up to sign, one of the $a_{i}$’s. Hence, the $\bar{a}_{i}$’s in (14) are effective parameters, and the function in (3.1) is an effective steady-state function. Finally, the fact that $\phi$ is a positive parametrization with respect to (3.1) (as in Definition 2.2) follows directly from comparing equations (3.1) and (3.1). ∎

Remark 3.2** (Multistationarity depends on only $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$).**

Proposition 3.1 considers any network obtained by deleting any (or none) of the six reactions labeled by $k_{2}$, $k_{\mathrm{on}}$, $m_{1}$, $\ell_{2}$, ${\ell}_{\mathrm{on}}$, $n_{2}$. Nonetheless, the resulting steady-state parametrization (3.1) depends on $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ but not any of the other rate constants. Thus, multistationarity for these irreversible networks depends only on whether the network contains $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ (see Theorem 4.6).

3.3 The reduced ERK network

In the previous subsection, we consider irreversible versions of the ERK network. Now we further reduce the network by additionally removing some “intermediate complexes” (namely, $S_{10}E$ and $S_{01}F$). These operations yield the reduced ERK network in Figure 2. Note that in the process of removing intermediates, the reactions $m_{2}$ and $m_{3}$ (similarly, $n_{1}$ and $n_{3}$) are collapsed into a single reaction labeled $m$ (respectively, $n$). A biological motivation for collapsing these reactions is the fact that intermediates are usually short-lived, so the simpler model may approximate the dynamics well.

Our notion of removing intermediates matches that of Feliu and Wiuf [20], who initiated the recent interest in the question of when dynamical properties are preserved when intermediates are added or removed (e.g., $S_{10}+E\leftrightarrows S_{10}E\to S_{11}+E$ versus $S_{10}\to S_{11}$). Our work, therefore, fits into this circle of ideas [7, 36, 41].

In the reduced ERK network, the remaining 10 rate constants are as follows: $k_{1},k_{3},k_{\mathrm{cat}},k_{\mathrm{off}},m,\ell_{1},\ell_{3},{\ell}_{\mathrm{cat}},{\ell}_{\mathrm{off}},n$. Letting $x_{1},x_{2},\ldots,x_{10}$ denote the species concentrations in the order given in Table 3, the resulting mass-action kinetics ODEs are as follows:

[TABLE]

[TABLE]

The 3 conservation equations are:

[TABLE]

Proposition 3.3** (Steady-state parametrization for reduced ERK network).**

The reduced ERK network (Figure 2) admits an effective steady-state function $h_{c,a}:\mathbb{R}^{10}_{>0}\to\mathbb{R}^{10}$ given by:

[TABLE]

Moreover, with respect to this effective steady-state function, the positive steady states admit the following positive parametrization:

[TABLE]

given by

[TABLE]

In particular, the image of $\phi$ is the following set of pairs of positive steady states and rate constants:

[TABLE]

Here, $k$ denotes the vector $(k_{1},k_{3},k_{\mathrm{cat}},k_{\mathrm{off}},m,\ell_{1},\ell_{3},{\ell}_{\mathrm{cat}},{\ell}_{\mathrm{off}},n)$.

Proof.

Let $W$ denote the conservation-law matrix arising from the conservation laws (3.3) for the reduced ERK network. Then $I=\{1,2,8\}$ is the set of indices of the first nonzero coordinates of the rows of $W$. We take $\mathbb{Q}(k_{\mathrm{cat}},k_{\mathrm{off}})$-linear combinations of the $f_{i}$’s in (25), where $i\notin I$, to obtain the following binomials in the $x_{i}$’s:

[TABLE]

Consider the (above) linear transformation from $f_{i}$ to $h_{i}$ $(i\not\in I)$. Let $M$ denote the corresponding matrix representation ($M$ plays the role of the matrix denoted by $M(\kappa)$ in Definition 2.1). It is straightforward to check that $\det M=k_{\mathrm{cat}}^{2}$, which is positive when $k_{\mathrm{cat}}>0$.

Consider the reparametrization map $\bar{a}:\mathbb{R}^{10}\to\mathbb{R}^{10}$ defined by the identity map (and so is surjective). Then $\bar{a}$, together with the conservation-law matrix $W$ and the matrix $M$, yield (as in Definition 2.1333In this case, Definition 2.1(iii)(b) requires every nonconstant coefficient in the effective steady-state function (32) to be a rational-number multiple of one of the rate constants. However, for the non-conservation-law equations in (32), many of the non-constant coefficients – such as $k_{\mathrm{off}}{\ell}_{\mathrm{cat}}$ – are not rational-number multiples of one of the rate constants. Nonetheless, these coefficients are all polynomials in the rate constants, and the relevant results in [17] hold in that generality.) the effective steady-state function $h_{c,a}(x)$ given in (32).

To show that $\phi$ is a positive steady-state parametrization with respect to (32), as in Definition 2.2, it suffices to show the following claim:

Claim: For every $(k^{*};x^{*})\in\mathbb{R}^{10+10}_{>0}$, the steady-state condition holds – namely, $h_{i}(k^{*};x^{*})=0$ for all $i\in\{3,4,5,6,7,9,10\}$ – if and only if $\phi(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})=(k^{*};x^{*})$.

For the “$\Rightarrow$” direction, assume $h_{i}(k^{*};x^{*})=0$ for all $i$. Then $h_{9}(k^{*};x^{*})=0$ implies that

[TABLE]

In other words, $\lambda_{\rm cat}$ – when evaluated at $(k_{\mathrm{cat}},k_{\mathrm{off}},{\ell}_{\mathrm{off}};x)=(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})$ – equals ${\ell}_{\mathrm{cat}}^{*}$. Next, the equality $h_{3}(k^{*};x^{*})=0$ implies that

[TABLE]

where the final equality follows from equation (35). Thus, the expression for $k_{1}$ given after (3.3) – when evaluated at $(k_{\mathrm{cat}},k_{\mathrm{off}},{\ell}_{\mathrm{off}};x)=(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})$ – equals $k_{1}^{*}$.

Similarly, the equality $h_{4}(k^{*};x^{*})=0$ (respectively, $h_{5}(k^{*};x^{*})=0$, $h_{6}(k^{*};x^{*})=0$, $h_{7}(k^{*};x^{*})=0$, or $h_{10}(k^{*};x^{*})=0$) implies that $\kappa_{3}$ (respectively, $m$, $\ell_{1}$, $\ell_{3}$, or $n$) – when evaluated at $(k_{\mathrm{cat}},k_{\mathrm{off}},{\ell}_{\mathrm{off}};x)=(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})$ – equals $k_{3}^{*}$ (respectively, $m^{*}$, $\ell^{*}_{1}$, $\ell^{*}_{3}$, or $n^{*}$). Thus, $\phi(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})=(k^{*};x^{*})$.

The “$\Leftarrow$” direction is similar. Assume $\phi(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})=(k^{*};x^{*})$. That is, the expressions for $k_{1}$, $k_{3}$, $m$, $\ell_{1}$, $\ell_{3}$, ${\ell}_{\mathrm{cat}}$, and $n$ evaluate to, respectively, $k_{1}^{*}$, $k_{3}^{*}$, $m^{*}$, $\ell_{1}^{*}$, $\ell_{3}^{*}$, ${\ell}_{\mathrm{cat}}^{*}$, and $n^{*}$, when $(k_{\mathrm{cat}},k_{\mathrm{off}},{\ell}_{\mathrm{off}};x)=(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*};x^{*})$. In particular, equation (35) holds, and so $h_{9}(k^{*};x^{*})=0$. Similarly, $h_{i}(k^{*};x^{*})=0$ for all other $i$ (here we also use equation (35)). ∎

Remark 3.4**.**

The proof of Proposition 3.3 proceeds by performing linear operations on the steady-state polynomials to yield binomials $g_{i}$, and then solving for one $k_{j}$ from each binomial to obtain the parametrization (33). This is similar in spirit to – but more general than – the approach prescribed in [17, §4] for “linearly binomial” networks. Also, our linear operations were found “by hand”, and so an interesting future direction is to develop efficient and systematic approaches to finding such operations leading to binomials.

Remark 3.5**.**

The proof of Proposition 3.3 shows that the “steady-state ideal” (the ideal generated by the right-hand sides of the ODEs) of the reduced ERK network is generated by the binomials $g_{i}$. This network, therefore has, “toric steady states” [39]. In contrast, the steady-state ideal of the full ERK network is not a binomial ideal (it is straightforward to check this computationally, e.g., using the Binomials package in Macaulay2 [23]). As for the irreversible versions of the ERK network, when the reactions with rate constants $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ are deleted, we see from (3.1) that the steady-state ideal becomes binomial. Hence, irreversible ERK networks that are missing both $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ are “linearly binomial” as in [17, §4].

Remark 3.6**.**

All networks considered in this section are conservative, which can be seen from the conservation laws (11) for the full and irreversible ERK networks, and (3.3) for the reduced ERK network. Also for these networks, there are no boundary steady states in any compatibility class (it is straightforward to check this using results from [1] or [43]).

4 Main Results

Each ERK network we investigated admits oscillations via a Hopf bifurcation (Section 4.1). Bistability, however, is more subtle (Section 4.2).

4.1 Oscillations

The full ERK system (Figure 1) exhibits oscillations for some values of the rate constants [40]. We now investigate oscillations in the fully irreversible and reduced ERK networks.

4.1.1 Fully irreversible ERK network

As shown in Figure 3, the fully irreversible ERK network admits oscillations. That figure was generated using the following rate constants:

[TABLE]

These rate constants (37) come from the ones that Rubinstein et al. showed generate oscillations for the full ERK network [40, Table 2] (we simply ignore their rate constants for the six deleted reactions). The approximate initial species concentrations used to generate Figure 3 are as follows (see supplementary file ERK-Matcont.txt):

[TABLE]

In Figure 3, we notice some peculiarities in the graphs $x_{i}(t)$ of the species concentrations. The species concentrations $x_{1}$ and $x_{2}$ (corresponding to \ceS00 and \ceE, respectively) peak dramatically, while $x_{3}$ and $x_{6}$ (\ceF and \ceS01F) stabilize momentarily at each peak. Also, each of $x_{1},x_{2},x_{3},x_{4},x_{5},x_{10},x_{11}$ deplete for some time in each period, whereas $x_{12}$ (\ceS11) never depletes. Finally, the graphs of the pairs $x_{1}$ and $x_{2}$ are qualitatively similar, and also the pair $x_{3}$ and $x_{6}$, the pair $x_{4}$ and $x_{5}$, and the pair $x_{10}$ and $x_{11}$.

Going beyond the fully irreversible ERK network, all other irreversible ERK networks – those obtained from the full ERK network by deleting one or more the reactions $k_{2},k_{\mathrm{on}},m_{1},\ell_{2},{\ell}_{\mathrm{on}},n_{2}$ – also admit oscillations. This claim follows from a result of Banaji that “lifts” oscillations when one or more reactions are made reversible [3, Proposition 4.1].

4.1.2 Reduced ERK network

We saw in the previous subsection that the fully irreversible ERK network exhibits oscillations. We now show that a simpler network - the reduced ERK network - also undergoes oscillations via a Hopf bifurcation. These oscillations are shown in Figure 4, and the rate constants that yield the corresponding Hopf bifurcation are specified in Theorem 4.3.

Compared to the oscillations for the irreversible ERK network (Figure 3), the oscillations in the reduced ERK network (Figure 4) are more uniform. Also, the period of oscillation is much shorter, and the amplitudes for species $x_{3}$, $x_{8}$, and $x_{10}$ are small (this may be due to the choice of rate constants). Finally, three of the six species shown do not deplete completely, whereas nearly all the species of the fully irreversible ERK do deplete in each period.

We discovered oscillations by finding a Hopf bifurcation. How we found this bifurcation – via the Hopf-bifurcation criterion in Section 2.4 – is the focus of the rest of this subsection.

Proposition 4.1** (Hopf criterion for reduced ERK).**

Consider the reduced ERK network, and let the polynomials $f_{i}$ denote the right-hand sides of the resulting ODEs, as in (25). Let $\hat{\kappa}:=(k_{\mathrm{cat}},k_{\mathrm{off}},{\ell}_{\mathrm{off}})$ and $x:=(x_{1},x_{2},\dots,x_{10})$, and let $\phi$ be the steady-state parametrization (33). Then the following is a univariate, degree-7 polynomial in $\lambda$, with coefficients in $\mathbb{Q}(x)[\hat{\kappa}]$:

[TABLE]

Now let $\mathfrak{h}_{i}$, for $i=4,5,6$, denote the determinant of the $i$-th Hurwitz matrix of the polynomial $q(\lambda)$ in (39). Then the following are equivalent:

1. (1)

there exists a rate-constant vector $\kappa^{*}\in\mathbb{R}^{10}_{>0}$ such that the resulting system (25) exhibits a simple Hopf bifurcation, with respect to $k_{\mathrm{cat}}$, at some $x^{*}\in\mathbb{R}^{10}_{>0}$, and 2. (2)

there exist $x^{*}\in\mathbb{R}^{10}_{>0}$ and $\hat{\kappa}^{*}\in\mathbb{R}^{3}_{>0}$ such that

[TABLE]

Moreover, given $\hat{\kappa}^{*}$ and $x^{*}$ as in (2), a simple Hopf bifurcation with respect to $k_{\mathrm{cat}}$ occurs at $x^{*}$ when the rate constants are taken to be $\kappa^{*}:=\widetilde{\pi}(\phi(\hat{\kappa}^{*};x^{*}))$. Here, $\widetilde{\pi}:\mathbb{R}_{>0}^{10}\times\mathbb{R}_{>0}^{10}\to\mathbb{R}_{>0}^{10}$ is the natural projection to the first 10 coordinates.

Proof.

The fact that $q(\lambda)$ is a degree-7 polynomial follows from Lemma 2.7, and the fact that its coefficients are in $\mathbb{Q}(x)[\hat{\kappa}]$ follows from inspecting equations (25) and (33). The rest of the result will follow immediately from Theorem 2.8 and Proposition 3.3, once we prove that $\mathfrak{h}_{1}$, $\mathfrak{h}_{2}$, $\mathfrak{h}_{3}$, and the constant term of $q(\lambda)$ are all positive when evaluated at any $(\hat{\kappa};x)\in\mathbb{R}^{3}_{>0}\times\mathbb{R}^{10}_{>0}$. Indeed, this is shown in the supplementary file reducedERK-hopf.mw. (In fact, even before substituting the parametrization $(\kappa;x)=\phi(\hat{\kappa};x)$, the corresponding Hurwitz determinants are already positive polynomials.) ∎

Remark 4.2**.**

Note that $k_{\mathrm{cat}}$ is the only free parameter, so it is the natural bifurcation parameter.

We now prove that the reduced ERK network gives rise to a Hopf bifurcation.

Theorem 4.3** (Hopf bifurcation in reduced ERK).**

The reduced ERK network exhibits a simple Hopf bifurcation with respect to the bifurcation parameter $k_{\mathrm{cat}}$ at the following point:

[TABLE]

when the rate constants are as follows:

[TABLE]

*Here, $\phi$ is the parametrization (33), and $\widetilde{\pi}$ is the projection to the first 10 coordinates. *

Proof.

By Proposition 4.1, we need only show that the inequalities and equality in (40) are satisfied at $x=x^{*}$ (with $x^{*}$ given in the statement of the theorem) and $\hat{\kappa}=\hat{\kappa}^{*}=(9,1,1)$. These are verified in the supplementary file reducedERK-hopf.mw. ∎

Remark 4.4**.**

The Hopf bifurcation given in Theorem 4.3 was found by analyzing the Newton polytopes of $\mathfrak{h_{4}}$, $\mathfrak{h_{5}}$, and $\mathfrak{h_{6}}$. The theory behind this approach is presented in Appendix B, and the steps we took to find the Hopf bifurcation are listed in Appendix C. We include these appendices for readers who wish to apply similar approaches to other systems.

4.2 Bistability

Although the full ERK network is bistable [40], we now prove that the reduced ERK network is not bistable (Proposition 4.5). As for irreversible ERK networks, some of them are bistable, and we show that bistability is controlled by the two reactions $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ (Theorem 4.6).

Proposition 4.5**.**

*The reduced ERK network is not multistationary, and hence not bistable. *

Proof.

Let $\mathcal{N}$ denote the reduced ERK network. By definition and Proposition 3.3, we obtain the following critical function for $\mathcal{N}$:

[TABLE]

where $\hat{a}=(k_{\mathrm{cat}},k_{\mathrm{off}},{\ell}_{\mathrm{off}})$, the function $h_{c,a}$ is as in (32), and $\phi(\hat{a};x)$ is as in (33).

This critical function $C(\hat{a};x)$ (see the supplementary file reducedERK-noMSS.mw) is a rational function, where the denominator is the following monomial: $x_{1}x_{2}x_{3}x_{5}x_{6}x_{7}x_{8}x_{9}$. The numerator of $C(\hat{a};x)$ is the following polynomial, which is negative when evaluated at any $(\hat{a};x)\in\mathbb{R}_{>0}^{3}\times\mathbb{R}_{>0}^{10}$:

[TABLE]

Thus, the following holds for all $(\hat{a};x)\in\mathbb{R}_{>0}^{3}\times\mathbb{R}_{>0}^{10}$:

[TABLE]

where the final equality uses the fact that the stoichiometric matrix $N$ has rank $10-3=7$.

So, by Proposition 2.4 and the fact that $\mathcal{N}$ is conservative with no boundary steady states in any stoichiometric compatibility class (Remark 3.6), $\mathcal{N}$ is monostationary. Thus, $\mathcal{N}$ is not multistationary and so, by definition, is not bistable. ∎

Although the reduced ERK network is not bistable (Proposition 4.5), the next result shows that irreversible versions of the full ERK network are bistable, as long as one of the reactions labeled by $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ is present. That is, this result tells us which reactions can be safely deleted (in contrast to standard results concerning reactions that can be added) while preserving bistability.

Theorem 4.6** (Bistability in irreversible ERK networks).**

*Consider any network $\mathcal{N}$ obtained from the full ERK network by deleting one or more of the reactions corresponding to rate constants $k_{2},k_{\mathrm{on}},m_{1},\ell_{2},{\ell}_{\mathrm{on}},n_{2}$ (blue in Figure 1). Then the following are equivalent: *

1. (1)

$\mathcal{N}$* is multistationary,* 2. (2)

$\mathcal{N}$* is bistable, and* 3. (3)

$\mathcal{N}$* contains at least one of the reactions labeled by $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$.*

Proof.

By definition, every bistable network is multistationary, so $(2)\Rightarrow(1)$. We therefore need only show $(1)\Rightarrow(3)\Rightarrow(2)$. (All computations below are found in our supplementary file irreversibleERK.mw).

For $(1)\Rightarrow(3)$, we will prove $\neg(3)\Rightarrow\neg(1)$: Assume that $\mathcal{N}$ contains neither the reaction labeled by $k_{\mathrm{on}}$ nor the reaction ${\ell}_{\mathrm{on}}$. Our proof here is analogous to that of Proposition 4.5. By Proposition 3.1, we obtain a critical function, $C(\hat{a};x)$, for $\mathcal{N}$ of the form (41), where now $h_{c,a}$ is as in (3.1) (with $\mathds{1}_{k_{\mathrm{on}}}=\mathds{1}_{{\ell}_{\mathrm{on}}}=0$) and $\phi(\hat{a};x)$ is as in (3.1) (with $\hat{a}=(a_{2},a_{4})$).

Here, $\det{\rm Jac}(h_{c,a})$ is a rational function with denominator equal to $k_{\mathrm{off}}x_{2}(n_{2}+n_{3}){\ell}_{\mathrm{cat}}l_{3}k_{3}m_{3}$, which is always positive. The numerator is a polynomial of degree 5 in the variables $x_{2},x_{3}$, and $x_{9}$ with coefficients that are always negative (see the supplementary file). The critical function $C(\hat{a};x)$ is obtained by substituting the positive parametrization into $\det{\rm Jac}(h_{c,a})$. Hence, for all $(\hat{a};x)\in\mathbb{R}_{>0}^{2}\times\mathbb{R}_{>0}^{12}$, the equality ${\rm sign}(C(\hat{a};x))=-1=(-1)^{{\rm rank}(N)}$ holds, because the stoichiometric matrix $N$ has rank $12-3=9$. So, by Proposition 2.4 (recall from Remark 3.6 that $\mathcal{N}$ is conservative with no boundary steady states in any stoichiometric compatibility class), $\mathcal{N}$ is not multistationary.

Now we show $(3)\Rightarrow(2)$, that is, if $\mathcal{N}$ contains at least one of the reactions labeled by $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$ then $\mathcal{N}$ is bistable. By symmetry (from exchanging in the network $E$, $S_{00}$, and $S_{01}$ with, respectively, $F$, $S_{11}$, and $S_{10}$), we may assume that $\mathcal{N}$ contains $k_{\mathrm{on}}$.

Consider the network $\mathcal{N}^{\prime}$ obtained from the full ERK network by deleting all reactions marked in blue in Figure 1, except for $k_{\mathrm{on}}$ (equivalently, we set $k_{2}=m_{1}=\ell_{2}={\ell}_{\mathrm{on}}=n_{2}=0$). We will show that the following total constants and rate constants yield bistability:

[TABLE]

Among the resulting three steady states (see the supplementary file), one of them is approximately:

[TABLE]

At the above steady state, the Jacobian matrix (of the system obtained from (10) by making the substitutions (4.2) and $k_{2}=m_{1}=\ell_{2}={\ell}_{\mathrm{on}}=n_{2}=0$) has three zero eigenvalues (due to the three conservation laws). For the remaining eigenvalues, the real parts are approximately:

[TABLE]

Thus, the nonzero eigenvalues have strictly negative real part, so the steady state is exponentially stable.

Another steady state is approximately

[TABLE]

At this steady state, the real part of the eigenvalues of the Jacobian matrix of the system are, in addition to the three zero eigenvalues, approximately as follows:

[TABLE]

This steady state is also exponentially stable. (A third steady state, not shown, is unstable.) Hence, $\mathcal{N}^{\prime}$ is bistable. Finally, as $\mathcal{N}^{\prime}$ is a subnetwork obtained from $\mathcal{N}$ by making some reversible reactions irreversible, then by [32, Theorem 3.1], bistability “lifts” from $\mathcal{N}^{\prime}$ to $\mathcal{N}$. Thus, $\mathcal{N}$ is bistable. ∎

We obtain the following immediate consequence of Theorem 4.6.

Corollary 4.7**.**

The fully irreversible ERK network is monostationary.

5 Maximum number of steady states

In the previous section, we saw that the full ERK network and some irreversible ERK networks (those with $k_{\mathrm{on}}$ or ${\ell}_{\mathrm{on}}$) are bistable, admitting two stable steady states in a stoichiometric compatibility class. The question arises, Do these networks admit three or more such steady states? We suspect not (Conjecture 5.10).

As a step toward resolving this problem, here we investigate the maximum number of positive steady states in ERK networks, together with some related measures we introduce, the maximum number of (non-boundary) complex-number steady states and the “mixed volume”. The mixed volume is always an upper bound on the number of complex steady states (Proposition 5.8), but we show these numbers are equal for ERK networks (Proposition 5.9).

5.1 Background and new definitions

Here we recall from [34] a network’s maximum number of positive steady states, and then extend the definition to allow for complex-number steady states.

Definition 5.1**.**

A network admits $k$ positive steady states (for some $k\in\mathbb{Z}_{\geq 0}$) if there exists a choice of positive rate constants so that the resulting mass-action system (1) has exactly $k$ positive steady states in some stoichiometric compatibility class (2).

In [34], $k=\infty$ was allowed when there are infinitely many steady states in a stoichiometric compatibility class. Here we do not allow $k=\infty$ so that we consider isolated roots only (as in Proposition 5.5 below).

Definition 5.2**.**

Let $G$ be a network with $s$ species, $m$ reactions, and a $d\times s$ conservation-law matrix $W$, which results in the system augmented by conservation laws $f_{c,\kappa}$, as in (3). The network $G$ admits $k$ steady states over $\mathbb{C}^{*}$ if there exists a choice of positive rate constants $\kappa\in\mathbb{R}^{m}_{>0}$ and a total-constant vector $c\in\mathbb{R}^{d}$ such that the system $f_{c,\kappa}=0$ has exactly $k$ solutions in $(\mathbb{C}^{*})^{s}=(\mathbb{C}\setminus\{0\})^{s}$.

It is straightforward to check that Definition 5.2 does not depend on the choice of $W$.

Definition 5.3**.**

The maximum number of positive steady states (respectively, maximum number of steady states over $\mathbb{C}^{*}$) of a network $G$ is the maximum value of $k$ for which $G$ admits $k$ positive steady states (respectively, $k$ steady states over $\mathbb{C}^{*}$).

Next we recall, from convex geometry, the concept of mixed volume, which we will apply to reaction networks. For background on convex and polyhedral geometry (such as polytopes and Minkowski sums), we direct the reader to Ziegler’s book [48]. In particular, for a polynomial $f=b_{1}x^{\sigma_{1}}+b_{2}x^{\sigma_{2}}+\dots+b_{\ell}x^{\sigma_{\ell}}\in\mathbb{R}[x_{1},x_{2},\dots,x_{s}]~{}$, where the exponent vectors $\sigma_{i}\in\mathbb{Z}^{s}$ are distinct and $b_{i}\neq 0$ for all $i$, the Newton polytope of $f$ is the convex hull of its exponent vectors: ${\rm Newt}(f)\coloneqq{\rm conv}\{\sigma_{1},~{}\sigma_{2},~{}\dots~{},~{}\sigma_{\ell}\}~{}\subseteq~{}\mathbb{R}^{s}.$

Definition 5.4**.**

Let $P_{1},P_{2},\ldots,P_{s}\subseteq\mathbb{R}^{s}$ be polytopes. The volume of the Minkowski sum $\lambda_{1}P_{1}+\lambda_{2}P_{2}+\ldots+\lambda_{s}P_{s}$ is a homogeneous polynomial of degree $s$ in nonnegative variables $\lambda_{1},\lambda_{2},\ldots,\lambda_{s}$. In this polynomial, the coefficient of $\lambda_{1}\lambda_{2}\cdots\lambda_{s}$, denoted by $\mathrm{Vol}(P_{1},P_{2},\ldots,P_{s})$, is the mixed volume of $P_{1},P_{2},...,P_{s}$.

The mixed volume counts the number of solutions in $(\mathbb{C}^{*})^{s}$ of a generic polynomial system.

Proposition 5.5** (Bernstein’s theorem [5]).**

Consider $s$ real polynomials $g_{1},g_{2},\dots,g_{s}\in\mathbb{R}[x_{1},x_{2},\dots,x_{s}]$. Then the number of isolated solutions in $(\mathbb{C}^{*})^{s}$, counted with multiplicity, of the system $g_{1}(x)=g_{2}(x)=\cdots=g_{s}(x)=0$ is at most $\mathrm{Vol}(\mathrm{Newt}(g_{1}),\ldots,\mathrm{Newt}(g_{s}))$.

Definition 5.6**.**

Let $G$ be a network with $s$ species, $m$ reactions, and a $d\times s$ conservation-law matrix $W$, which results in the system augmented by conservation laws $f_{c,\kappa}$, as in (3). Let $c^{*}\in\mathbb{R}^{d}_{\neq 0}$, and let $\kappa^{*}\in\mathbb{R}^{m}_{>0}$ be generic. Let $P_{1},P_{2},\ldots,P_{s}\subset\mathbb{R}^{s}$ be the Newton polytopes of $f_{c^{*},\kappa^{*},1},f_{c^{*},\kappa^{*},2},\ldots,f_{c^{*},\kappa^{*},s}$, respectively. The mixed volume of $G$ (with respect to $W$) is the mixed volume of $P_{1},P_{2},\ldots,P_{s}$.

A closely related definition is introduced and analyzed by Gross and Hill [25].

Remark 5.7**.**

The mixed volume (Definition 5.6) is well defined. Indeed, it is straightforward to check that the exponents appearing in $f_{c^{*},\kappa^{*}}$ are the same as long as $c^{*}\in\mathbb{R}^{d}_{\neq 0}$ and $\kappa^{*}$ is chosen generically (so that no coefficients of $f_{c^{*},\kappa^{*}}$ vanish, or equivalently certain linear combinations of the $\kappa_{j}$’s do not vanish).

5.2 Results

Every positive steady state is a steady state over $\mathbb{C}^{*}$. Also, the mixed volume pertains to polynomial systems with the same supports (i.e., the exponents that appear in each polynomial) as the augmented system $f_{c,\kappa}=0$ (but without constraining the coefficients to come from a reaction network). We obtain, therefore, the bounds in the following result:

Proposition 5.8**.**

For every network, the following inequalities hold among the maximum number of positive steady states, the maximum number of steady states over $\mathbb{C}^{*}$, and the mixed volume of the network (with respect to any conservation-law matrix):

[TABLE]

Proof.

This result follows from Proposition 5.5 and Definitions 5.1–5.3. ∎

We investigate the numbers in Proposition 5.8 for ERK networks in the following result.

Proposition 5.9**.**

Consider four ERK networks: the full ERK network, the full ERK network with the reaction $k_{\mathrm{on}}$ removed, the fully irreversible network, and the reduced network. For these networks, the following numbers (or bounds on them) are given in Table 4: the maximum number of positive steady states, the maximum number of steady states over $\mathbb{C}^{*}$, and the mixed volume of the network (with respect to the consveration laws (11) or (3.3)).

Proof.

The results on the mixed volume were computed using the PHCpack [26] package in Macaulay2 [23]. See the supplementary file ERK-mixedVol.m2.

The mixed volume is an upper bound on the maximum number of steady states over $\mathbb{C}^{*}$ (Proposition 5.8), so we need only show that each network admits the number shown in Table 4 for steady states over $\mathbb{C}^{*}$.

The full ERK network admits 7 steady states over $\mathbb{C}^{*}$ (including 3 positive steady states) [17, Example 3.18]. Next, we consider the remaining three networks (see the supplementary file ERK-MaxComplexNumber.nb).

For the full ERK network with $k_{on}=0$, when $(c_{1},c_{2},c_{3})=(1,2,3)$ and $(k_{1},k_{2},$ $k_{3},k_{\mathrm{cat}},k_{\mathrm{on}},k_{\mathrm{off}},{\ell}_{1},{\ell}_{2},{\ell}_{3},{\ell}_{\mathrm{cat}},{\ell}_{\mathrm{on}},{\ell}_{\mathrm{off}},m_{1},m_{2},m_{3},n_{1},n_{2},n_{3})=(3,25,1,5,0,6,5,23,11,13,$ $43,41,12,7,8,12,31,21),$ we obtain $5$ steady states over $\mathbb{C}^{*}$, three real and one complex-conjugate pair, which are approximately as follows:

[TABLE]

For the fully irreversible ERK network, when $(c_{1},c_{2},c_{3})=(1,2,3)$ and $(k_{1},k_{3},k_{\mathrm{cat}},$ $k_{\mathrm{off}},{\ell}_{1},{\ell}_{3},{\ell}_{\mathrm{cat}},{\ell}_{\mathrm{off}},m_{2},m_{3},n_{1},n_{3})=(3,1,5,6,5,11,13,41,7,8,12,21)$, there are 3 steady states over $\mathbb{C}^{*}$, all real, with approximate values:

[TABLE]

For the reduced ERK network, let $(c_{1},c_{2},c_{3})=(1,2,3)$ and $(k_{1},k_{3},k_{\mathrm{cat}},k_{\mathrm{off}},m,n,{\ell}_{1},$ ${\ell}_{3},{\ell}_{\mathrm{cat}},{\ell}_{\mathrm{off}})=(3,4,1,5,6,8,7,11,12,5)$. We obtain $3$ steady states over $\mathbb{C}^{*}$, all real, which are approximately:

[TABLE]

Finally, we examine the maximum number of positive steady states. We already saw that the fully irreversible and reversible networks are monostationary (Corollary 4.7 and Proposition 4.5, respectively). For the “partially irreversible” network, we saw in the proof of Theorem 4.6 that it admits 3 positive steady states. As for the full network, as noted above, 3 positive steady states were shown in [17, Example 3.18]. ∎

Table 4 suggests that the mixed volume is a measure of the complexity of a network. The full ERK network is multistationary, and its mixed volume is 7. The mixed volume drops to 5 when $k_{\mathrm{on}}=0$. When the network is further simplified to the fully irreversible, or even to the reduced ERK network, the mixed volume becomes 3, and bistability is lost as well.

Finally, we conjecture that the bounds in Table 4 are strict, and ask about stability.

Conjecture 5.10**.**

For the full ERK network and the full ERK network with $k_{\mathrm{on}}=0$, the maximum number of positive (respectively, positive stable) steady states is 3 (respectively, 2).

6 Discussion

Phosphorylation plays a key role in cellular signaling networks, such as mitogen-activated protein kinase (MAPK) cascades, which enable cells to make decisions (to differentiate, proliferate, die, and so on) [8]. This decision-making role of MAPK cascades suggests that they exhibit switch-like behavior, i.e., bistability. Indeed, bistability in such cascades has been seen in experiments [2, 6]. Oscillations also have been observed [29, 30], hinting at a role in timekeeping. Indeed, multisite phosphorylation is the main mechanism for establishing the 24-hour period in eukaryotic circadian clocks [38, 46].

These experimental findings motivated the questions we pursued. Specifically, we investigated robustness of oscillations and bistability in models of ERK regulation by dual-site phosphorylation. Bistability, we found, is quickly lost when reactions are made irreversible. Indeed, bistability is characterized by the presence of two specific reactions. Oscillations, in contrast, persist even as the network is greatly simplified. Indeed, we discovered oscillations in the reduced ERK network. Moreover, this network has the same number of reactions (ten) as the mixed-mechanism network which Suwanmajo and Krishnan surmised “could be the simplest enzymatic modification scheme that can intrinsically exhibit oscillation” [44, §3.1]. Our reduced ERK network, therefore, may also be such a minimal oscillatory network.

Returning to our bistability criterion (Theorem 4.6), recall that this result elucidates which reactions can be safely deleted while preserving bistability – in contrast to standard results concerning reactions that can be added [4, 20, 32]. We desire more results of this type, so we comment on how we proved our result. The key was the special form of the steady-state parametrization. In particular, following [17], our parametrizations allow both species concentrations and rate constants to be solved (at steady state) in terms of other variables. Additionally, a single parametrizations specialized (by setting rates to zero for deleted reactions) to obtain parametrizations for a whole family of networks. Together, these properties gave us access to new information on how bistability is controlled. We are interested, therefore, in the following question: Which networks admit a steady-state parametrization that specializes for irreversible versions of the network?

Our results on oscillations were enabled by new mathematical approaches to find Hopf bifurcations. Specifically, building on [12], we gave a Hopf-bifurcation criterion for networks admitting a steady-state parametrization. Additionally, we successfully applied this criterion to the reduced ERK network by analyzing the Newton polytopes of certain Hurwitz determinants. We expect these techniques to apply to more networks.

Finally, our work generated a number of open questions. First, what are the mixed volumes of irreversible versions of the ERK network (beyond those shown in Table 4)? In particular, is there a mixed-volume analogue of our bistability criterion, which is in terms of the reactions $k_{\mathrm{on}}$ and ${\ell}_{\mathrm{on}}$? And, what is the maximum number of (stable) steady states in the full ERK network (Conjecture 5.10)? Progress toward these questions will yield further insight into robustness of bistability and oscillations in biological signaling networks.

Acknowledgements

NO, AS, and XT were partially supported by the NSF (DMS-1752672). AT was partially supported by the Independent Research Fund Denmark. The authors thank Elisenda Feliu for insightful comments on an earlier draft, and thank Carsten Conradi, Elizabeth Gross, Cvetelina Hill, Maya Mincheva, Stanislav Shvartsman, Frank Sottile, Elise Walker, and Timo de Wolff for helpful discussions. This project was initiated while AT was a visiting scholar at Texas A&M University, and while XT was hosted by ICERM. We thank Texas A&M University and ICERM for their hospitality.

Appendix A Files in the Supporting Information

Table 5 lists the files in the Supporting Information, and the result/proof each file supports. All files can be found at the online repository: https://github.com/neeedz/ERK.

Appendix B Newton-polytope method

Here we show how analyzing the Newton polytopes of two polynomials can reveal whether there is a positive point at which one polynomial is positive and simultaneously the other is zero (Proposition B.2 and Algorithm 1). In Appendix C, we show how we used this approach, which we call the Newton-polytope method, to find a Hopf bifurcation leading to oscillations in the reduced ERK network (in Theorem 4.3).

Notation B.1**.**

Consider a polynomial $f=b_{1}x^{\sigma_{1}}+b_{2}x^{\sigma_{2}}+\dots+b_{\ell}x^{\sigma_{\ell}}\in\mathbb{R}[x_{1},x_{2},\dots,x_{s}]$, where the exponent vectors $\sigma_{i}\in\mathbb{Z}_{\geq 0}^{s}$ are distinct and $b_{i}\neq 0$ for all $i$. A vertex $\sigma_{i}$ of $\mathrm{Newt}(f)$, the Newton polytope of $f$, is a positive vertex (respectively, negative vertex) if the corresponding monomial of $f$ is positive, i.e., $b_{i}>0$ (respectively, $b_{i}<0$). Also, $N_{f}(\sigma)$ denotes the outer normal cone of the vertex $\sigma$ of $\mathrm{Newt}(f)$, i.e., the cone generated by the outer normal vectors to all supporting hyperplanes of $\mathrm{Newt}(f)$ containing the vertex $\sigma$. Finally, for a cone $C$, let $\mathrm{int}(C)$ denote the relative interior of the cone.

For an extensive discussion on polytopes and normal cones, see [48].

Proposition B.2**.**

Let $f,g\in\mathbb{R}[x_{1},x_{2},\ldots x_{s}]$. Assume that $\alpha$ is a positive vertex of ${\rm Newt}(f)$, $\beta_{+}$ is a positive vertex of ${\rm Newt}(g)$, and $\beta_{-}$ is a negative vertex of ${\rm Newt}(g)$. Then, if $\mathrm{int}(N_{f}(\alpha))\cap\mathrm{int}(N_{g}(\beta_{+}))$ and $\mathrm{int}(N_{f}(\alpha))\cap\mathrm{int}(N_{g}(\beta_{-}))$ are both nonempty, then there exists $x^{*}\in\mathbb{R}^{s}_{>0}$ such that $f(x^{*})>0$ and $g(x^{*})=0$.

To prove Proposition B.2 we use the following well-known lemma and its proof.

Lemma B.3**.**

*For a real, multivariate polynomial $f=b_{1}x^{\sigma_{1}}+b_{2}x^{\sigma_{2}}+\dots+b_{\ell}x^{\sigma_{\ell}}\in\mathbb{R}[x_{1},x_{2},\dots,x_{s}],$ if $\sigma_{i}$ is a positive vertex (respectively, negative vertex) of ${\rm Newt}(f)$, then there exists $x^{*}\in\mathbb{R}^{s}_{>0}$ such that $f(x^{*})>0$ (respectively, $f(x^{*})<0$). *

Proof.

Let $\sigma_{i}$ be a vertex of ${\rm Newt}(f)$. Pick $w=(w_{1},w_{2},\ldots,w_{s})$ in the relative interior of the outer normal cone $N_{f}(\sigma_{i})$, which exists because $\sigma_{i}$ is a vertex. Then, by construction, the linear functional $\langle w,-\rangle$ is maximized over the exponent-vectors $\sigma_{1},\sigma_{2},\dots,\sigma_{\ell}$ at $\sigma_{i}$. Thus, we have the following univariate “polynomial with real exponents” in $t$:

[TABLE]

So, for $t$ large, $\text{sign}(f(t^{w_{1}},t^{w_{2}},\ldots,t^{w_{s}}))=\text{sign}(b_{i})$. Note that $(t^{w_{1}},t^{w_{2}},\ldots,t^{w_{s}})\in\mathbb{R}^{s}_{>0}$. ∎

Our proof of Proposition B.2 is constructive, through the following algorithm, where we use the notation $f_{w}(t):=f(t^{w_{1}},t^{w_{2}},\ldots,t^{w_{s}})$, for $t\in\mathbb{R}$ and $w=(w_{1},w_{2},\dots,w_{s})\in\mathbb{R}^{s}$.

Proof of Proposition B.2.

Let $a_{+}x^{\alpha}$ be the term of $f$ corresponding to the vertex $\alpha$ of $\mathrm{Newt}(f)$, and similarly let $b_{+}x^{\beta_{+}}$ (respectively, $b_{-}x^{\beta_{-}}$) be the term of $g$ corresponding to the vertex $\beta_{+}$ (respectively, $\beta_{-}$) of $\mathrm{Newt}(g)$. Thus, $a_{+}>0$, $b_{+}>0$, and $b_{-}<0$. Let $\{a_{1},a_{2},\ldots,a_{d}\}\subseteq\mathbb{R}$ denote the remaining set of coefficients of $f$, so that $f=a_{+}x^{\alpha}+(a_{1}x^{\sigma_{1}}+a_{2}x^{\sigma_{2}}+\dots+a_{d}x^{\sigma_{d}})$, for some exponent vectors $\sigma_{i}\in\mathbb{Z}_{\geq 0}^{s}$.

Algorithm 1 terminates: First, $\ell$ and $m$ in line 2 exist by hypothesis. Also, $\tau_{\ell}$ and $\tau_{m}$ in lines 4–5 exist by the proof of Lemma B.3 and by construction. Next, $\min h(r)$ in line 8 exists because $h$ is a continuous univariate function defined on a compact interval.

By construction and because cones are convex, the vector $r\ell+(1-r)m$, which is a convex combination of $\ell$ and $m$, is in the relative interior of $N_{f}(\alpha)$ for all $r\in[0,1]$. Thus, $\langle r\ell+(1-r)m,~{}\alpha-\sigma_{i}\rangle>0$ for all $i=1,2,\dots,d$ and for all $r\in[0,1]$. This (together with a straightforward argument using continuity and compactness) implies the following:

[TABLE]

Next, let $\beta:=\inf_{r\in[0,1]}\langle r\ell+(1-r)m,~{}\alpha\rangle$. Then, for all $r\in[0,1]$ and $t>0$,

[TABLE]

In $\widetilde{f}(t)$, the term $a_{+}t^{\beta}$ dominates the other term, for $t$ large, so there exists $T^{*}>0$ such that $\widetilde{f}(t)\geq 0$ when $t\geq T^{*}$. So, by (B), the while loop in line 8 ends when $T\geq T^{*}$ (or earlier).

Algorithm 1 is correct: For $T$ fixed, the minimum of $\psi(r):=\left(g(T^{r\ell+(1-r)m})\right)^{2}$ over the compact set $[0,1]$ is attained, because $\psi$ is continuous. Next we show that this minimum value is 0, or equivalently that for $\chi(r):=g(T^{r\ell+(1-r)m})$ there exists some $r^{*}\in(0,1)$ such that $\chi(r^{*})=0$. Indeed, this follows from the Intermediate Value Theorem, because $\chi$ is continuous, $\chi(0)=g(T^{m})<0$ (because $T>\tau_{m}$), and $\chi(1)=g(T^{\ell})>0$ (because $T>\tau_{\ell})$.

Finally, the inequality $f(T^{r^{*}\ell+(1-r^{*})m})>0$ holds by construction of $T$, so defining $x^{*}:=T^{r^{*}\ell+(1-r^{*})m}\in\mathbb{R}^{s}_{>0}$ yields the desired vector satisfying $f(x^{*})>0$ and $g(x^{*})=0$. ∎

Appendix C Using the Newton-polytope method

Here we show how we used Algorithm 1 to find the Hopf bifurcation in Theorem 4.3. (For details, see the supplementary files reducedERK-hopf.mw and reducedERK-cones.sws). Recall from the proof of that theorem, that our goal was to find some $x^{*}\in\mathbb{R}^{10}_{>0}$ and $\hat{\kappa}^{*}=(k_{\mathrm{cat}}^{*},k_{\mathrm{off}}^{*},{\ell}_{\mathrm{off}}^{*})\in\mathbb{R}^{3}_{>0}$ satisfying the following conditions from Proposition 4.1:

[TABLE]

Step One. Specialize some of the parameters: set $k_{\rm off}=\ell_{\rm off}=1$ and $x_{3}=x_{4}=x_{6}=x_{7}=x_{9}=x_{10}=1$. (Otherwise, $\mathfrak{h}_{5}$ and $\mathfrak{h}_{6}$ are too large to be computed.)

Step Two. Do a change of variables: let $y_{i}=1/x_{i}$ for $i=1,2,5,8$. These variables $x_{i}$ were in the denominator, so switching to the variables $y_{i}$ yield polynomials.

Let $\mathcal{H}_{4}$, $\mathcal{H}_{5}$, and $\mathcal{H}_{6}$ denote the resulting polynomials in $\mathbb{Q}[k_{\mathrm{cat}},y_{1},y_{2},y_{5},y_{8}]$ after performing Steps One and Two. Accordingly, our updated goal is to find $(k_{\mathrm{cat}}^{*},y^{*}_{1},y^{*}_{2},y^{*}_{5},y^{*}_{8})\in\mathbb{R}^{5}_{>0}$ at which $\mathcal{H}_{4}$ and $\mathcal{H}_{5}$ are positive and $\mathcal{H}_{6}$ is zero. (In a later step, we must also check the partial-derivative condition in (44).)

Step Three. Apply (a straightforward generalization of) Algorithm 1 as follows.

1. (i)

Find a positive vertex of $\mathcal{H}_{4}$ and a positive vertex of $\mathcal{H}_{5}$ whose outer normal cones intersect (denote the intersection by $C$), and a positive vertex and a negative vertex of $\mathcal{H}_{6}$ (denote their outer normal cones by $D_{+}$ and $D_{-}$, respectively) for which:

1. (a)

the intersection $D_{+}\cap D_{-}$ is 4-dimensional, and 2. (b)

the intersections $C\cap D_{+}$ and $C\cap D_{-}$ are both 5-dimensional. 2. (ii)

By Proposition B.2, a vector $(k_{\mathrm{cat}}^{*},y^{*}_{1},y^{*}_{2},y^{*}_{5},y^{*}_{8})$ that accomplishes our updated goal, is guaranteed. To find such a point, we follow Algorithm 1 to obtain $k_{\mathrm{cat}}^{*}=729,y_{1}^{*}\approx 16.79978292,y_{2}^{*}\approx 453.5941389,y_{5}^{*}\approx 6.587368051$, and $y_{8}^{*}\approx 80675.77181$.

Recall the specializations in Step One and change of variables in Step Two, to obtain $\hat{\kappa}=(729,1,1)$ and

[TABLE]

Step Four. Verify that the conditions in (44) hold.