Dynamics of ERK regulation in the processive limit

Abstract

We consider a model of extracellular signal-regulated kinase regulation by dual-site phosphorylation and dephosphorylation, which exhibits bistability and oscillations, but loses these properties in the limit in which the mechanisms underlying phosphorylation and dephosphorylation become processive. Our results suggest that anywhere along the way to becoming processive, the model remains bistable and oscillatory. More precisely, in simplified versions of the model, precursors to bistability and oscillations (specifically, multistationarity and Hopf bifurcations, respectively) exist at all “processivity levels”. Finally, we investigate whether bistability and oscillations can exist together.

Appendices

A files in the supporting information

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

Table 5 Supporting Information files and the results they support. Here, Text* indicates an output file from using the Julia package HomotopyContinuation.jl (Breiding and Timme 2018)

B Procedure to study multistationarity numerically

Here we describe the procedure we used in Sect. 4.3 for numerically studying multistationarity in the minimally bistable ERK network at various processivity levels $p_k$ and $p_{\ell }$.

We begin by mirroring the analysis of Sect. 4.2. Specifically, we use the parameters given in (14) to study the critical function $C(\kappa ,{{\hat{x}}})$ for $x_1=x_2=T$ and $x_3=1$. Due to this choice of $\kappa$ and ${{\hat{c}}}$, the critical function is a (rational) function of $p_k$, $p_\ell$, and T only, i.e., $C(\kappa ,{{\hat{x}}}) \equiv C(p_k,p_\ell ,T)$. The numerator is the following polynomial:

$$\begin{aligned} \begin{aligned}&q(p_k,p_\ell ,T) = -p_k \left( (3-2 p_\ell ) p_\ell +p_k \left( 1-3 p_\ell +2 p_\ell ^2\right) \right) \\&\quad +\left( -5 p_\ell ^2+p_k p_\ell (-9+11 p_\ell )+p_k^3 \left( -1+3 p_\ell -2 p_\ell ^2\right) -p_k^2 \left( 3-3 p_\ell +p_\ell ^2\right) \right) T \\&\quad +\left( -8 p_\ell ^2+p_k p_\ell (-13+9 p_\ell )+p_k^3 \left( -1+p_\ell -p_\ell ^2\right) +p_k^2 \left( -6+2 p_\ell +7 p_\ell ^2\right) \right) T^2 \\&\quad +\left( -3 p_\ell ^2-p_k p_\ell (5+8 p_\ell )+p_k^3 \left( -4+3 p_\ell +p_\ell ^2\right) +p_k^2 \left( -3-10 p_\ell +13 p_\ell ^2\right) \right) T^3 \\&\quad +p_k \left( (5-8 p_\ell ) p_\ell +3 p_k \left( 1-5 p_\ell +p_\ell ^2\right) +p_k^2 \left( -7+4 p_\ell +2 p_\ell ^2\right) \right) T^4 \\&\quad +p_k \left( 3 p_\ell +p_k^2 \left( -3-3 p_\ell +2 p_\ell ^2\right) -p_k \left( -2+p_\ell +4 p_\ell ^2\right) \right) T^5 -2 (-1+p_k) p_k^2 p_\ell T^6 \end{aligned} \end{aligned}$$

(25)

As $0< p_k, p_\ell < 1$, the leading coefficient of $q(p_k,p_\ell ,T)$ as a polynomial in T is positive. Next, the steady-state parametrization $\phi$ from Proposition 3.1 is as follows (cf. eq. (6)):

$$\begin{aligned} x_1 = T,\; x_2 = T,\; x_3 =1,\; x_{4} = \frac{p_k T^2 (1+T)}{p_\ell +p_k p_\ell T},\; x_{5} = -\frac{p_k (-1+p_\ell ) T (1+T)}{p_\ell +p_k p_\ell T}, \nonumber \\ x_{6} = \frac{T-p_k T}{1+p_k T},\; x_{7} = -\frac{(-1+p_k) T (1+T)}{1+p_kT},\; x_{8} = -\frac{p_k (-1+p_\ell ) T (1+T)}{p_\ell +p_k p_\ell T},\nonumber \\ x_{9} = \frac{T-p_k T}{1+p_k T},\; x_{10} = -\frac{p_k (-1+p_\ell ) (1+T)}{p_\ell +p_k p_\ell T},\; x_{11} = T^2,\; x_{12} = \frac{p_k (1+T)}{p_\ell +p_k p_\ell T}\nonumber \\ \end{aligned}$$

(26)

Table 6 Values of $p_k$, $p_\ell$, and T used in Figs. 3 and 4

Table 7 Values of $p_k$, $p_\ell$, and T used in Fig. 5

To numerically study multistationarity for $p_k$, $p_\ell \rightarrow 1$, we proceed as follows:

• (i) Pick values of $0< {{\tilde{p}}}_k, {{\tilde{p}}}_\ell < 1$ and ${{\tilde{T}}}>0$ such that $q({{\tilde{p}}}_k, {{\tilde{p}}}_\ell , {{\tilde{T}}})>0$ (recall eq. (25)).

• (ii) Substitute into (26) the values of ${{\tilde{p}}}_k$, ${{\tilde{p}}}_\ell$, and ${{\tilde{T}}}$ from the previous step to obtain a steady state ${{\tilde{x}}}$.

• (iii) Compute, using (5), the total amounts $\tilde{c}_1$, ${{\tilde{c}}}_2$, and ${{\tilde{c}}}_3$ at ${{\tilde{x}}}$.

• (iv) Use Matcont with initial condition near ${{\tilde{x}}}$ and bifurcation parameter $c_2$, to obtain a bifurcation curve.

• (v) To compare curves corresponding to distinct ${{\tilde{p}}}_k$ and ${{\tilde{p}}}_\ell$, compute relative concentrations $\frac{x_i}{{{\tilde{c}}}_1}$ and $\frac{c_2}{{{\tilde{c}}}_2}$ that relate $x_i$ and $c_2$ to the ${{\tilde{x}}}$ and ${{\tilde{c}}}_2$ computed in steps (ii) and (iii) .

Step (v) is crucial for interpreting the numerical results obtained by the above procedure, because certain total amounts differ by orders of magnitude as $p_k$, $p_\ell \rightarrow 1$, and so it is more meaningful to compare values relative to the reference point ${{\tilde{x}}}$ obtained in step (ii) Figs. 6 and 7.