A Case Study on the Parametric Occurrence of Multiple Steady States

TL;DR

This paper applies advanced symbolic computation techniques to analyze the occurrence of multiple steady states in biological network models, specifically the MAPK network, revealing conditions for multistationarity with partial parameter knowledge.

Contribution

It introduces a novel combination of symbolic methods to determine multistationarity in high-dimensional biological models with partial parameter information.

Findings

Successfully determined multistationarity conditions for an 11-dimensional MAPK model.

Demonstrated the effectiveness of symbolic computation methods in complex biological systems.

Provided insights into parameter regions leading to multiple steady states.

Abstract

We consider the problem of determining multiple steady states for positive real values in models of biological networks. Investigating the potential for these in models of the mitogen-activated protein kinases (MAPK) network has consumed considerable effort using special insights into the structure of corresponding models. Here we apply combinations of symbolic computation methods for mixed equality/inequality systems, specifically virtual substitution, lazy real triangularization and cylindrical algebraic decomposition. We determine multistationarity of an 11-dimensional MAPK network when numeric values are known for all but potentially one parameter. More precisely, our considered model has 11 equations in 11 variables and 19 parameters, 3 of which are of interest for symbolic treatment, and furthermore positivity conditions on all variables and parameters.

Tables

Table 1. Table 1: The unique solution x ( 200 ) superscript 𝑥 200 x^{(200)} for k 19 = 200 subscript 𝑘 19 200 k_{19}=200 and the three solutions x 1 ( 500 ) subscript superscript 𝑥 500 1 x^{(500)}_{1} , …,  x 3 ( 500 ) subscript superscript 𝑥 500 3 x^{(500)}_{3} for k 19 = 500 subscript 𝑘 19 500 k_{19}=500 . We have actually computed real algebraic numbers, which are pairs of univariate polynomials and isolated intervals. For convenience we are giving machine float approximations here, which can be made arbitrarily precise.

	$x^{\left(200\right)}$	$=(90.6512,2.67311,10.4996,17.8545,35.9695,32.0501,0.0954536,15.5631,2.39331,0.641001,45.4331)$
	$x_1^{\left(500\right)}$	$=(17.6392,6.97675,367.57,36.6772,5.50874,12.811,0.511775,83.4416,8.06095,0.25622,2.73253)$
	$x_2^{\left(500\right)}$	$=(122.034,14.6721,234.974,14.5102,7.16952,35.064,0.42579,69.4223,7.43877,0.70128,15.2681)$
	$x_3^{\left(500\right)}$	$=(323.761,9.49621,37.1013,6.72938,13.6295,43.1428,0.127807,20.8381,3.21139,0.862856,61.4581)$

Table 2. Table 2: Coefficients c i subscript 𝑐 𝑖 c_{i} and d j subscript 𝑑 𝑗 d_{j} of polynomials occurring in Equations ( 6 ) and ( 7 ), respectively.

	$c_{10}$	$=351590934502740290936895033267017158736060313940693076650155371250411$
	$c_9$	$=-213699072852157674283997527746395583273033983170426080574800781989093156$
	$c_8$	$=25374851641220554774259605635053469432582109883965015804077119110958034090$
	$c_7$	$=12972493018300022707027639267804259251235991618029852880330004508564391594000$
	$c_6$	$=-8468945963692802414226427249726123493448372439778349029355636316929687020660000$
	$c_5$	$=2231098270337406450670301663172664333421440833875848621423683265663846533079600000$
	$c_4$	$=-376265008904112258290319173193792052014899485528994925965885895511831873444245100000$
	$c_3$	$=39262101548790869407057994985320156500968958361396178908180026842806643766783104000000$
	$c_2$	$=-2492623990743029234974354081270296106309603462451517057779877596842448287799337600000000$
	$c_1$	$=70978850735887473459176997186175978425873267246760023212940616924643171868478080000000000$
	$c_0$	$=-1062871192838985876948077114923898204990434138901495394834749613184670362810368000000000000$
	$d_6$	$=16838105723097694257603469$
	$d_5$	$=-24078605201553273505077988k_{19}+7723967969644977896148686580$
	$d_4$	$=8176202638735769127032169k_{19}^2-7723411665463544477701499460k_{19}+1232154357941338876156606812900$
	$d_3$	$=1465408757440589841803452380k_{19}^2-798169557586805582842481309800k_{19}+83152655240002767729550477640000$
	$d_2$	$=85462524901276846107251669400k_{19}^2-35266411401427656834572095140000k_{19}+2556805354853318332197489636000000$
	$d_1$	$=1631685649719702672282505500000k_{19}^2-721989571100461862477342320000000k_{19}+28843755938318780823218400000000000$
	$d_0$	$=-7013104139459910876520500000000000k_{19}.$

Table 3. Table 3: Triangular solution formulae valid for all positive k 19 subscript 𝑘 19 k_{19} excluding three isolated points

	$x_{11}$	$=-\frac{1}{60}{x}_2^2+\frac{1}{600}\left(10k_{19}-10x_1-37x_3+10x_4-2100\right)x_2-\frac{9}{200}{x}_3^2+\frac{1}{600}\left(-27x_1+27x_4+27k_{19}-4650\right)x_3$
		$-x_1+x_4+k_{19}-50$
	$x_{10}$	$=\frac{1}{150}{x}_2\left(x_2+x_3-x_4-k_{19}+x_1+150\right)$
	$x_9$	$=\frac{1}{18200}\left(69x_3+182x_2\right)\left(x_2+x_3-x_4-k_{19}+x_1+150\right)$
	$x_8$	$=\frac{15}{364}\left(x_2+x_3-x_4-k_{19}+x_1+150\right)x_3$
	$x_7$	$=50-\frac{2}{101}{x}_4{x}_1-x_4$
	$x_6$	$=\frac{2}{101}{x}_4{x}_1$
	$x_5$	$=x_2+x_3-x_4-k_{19}+x_1+150$
	$x_4$	$=\frac{2525000}{101x_2+1000x_1+50500}$
	$x_3$	$=\frac{-101x_2^3-(-101k_{19}+1101x_1+65650)x_2^2-(1000x_1^2+(-1000k_{19}+200500)x_1-50500k_{19}+5050000)x_2+150000x_1)}{101x_2^2+\left(1000x_1+50500\right)x_2}$
	$x_2$
	$x_1$	has at most 6 solutions for a given value of $k_{19}$, according to Equation (7).

Table 4. Table 4: Constraints on k 19 subscript 𝑘 19 k_{19} for solution formulae in Table 3 and Equation ( 7 ) to be valid

	$k_{19}$	$>0$
	polynomial in (6)	$\ne 0$
	$23197989433419579994929k_{19}^2-89407400615452409453098800k_{19}-4822419303419166525491149190000$	$\ne 0$
	$505465566622475867655547880786544637953790406059982726185509k_{19}^4$
	$-12725780456964391893178560515183873684222178969868366920505134120k_{19}^3$
	$+1175510330915205241831243213231417517003037315562884193657451445400k_{19}^2$
	$-281867359883676159811192082978541193600292804324596911878337972560000k_{19}$
	$-42434363570215587465668423701563932185051066892741207931879307200000000$	$\ne 0$

Equations

$$\dot{x_1}=$$

$$\dot{x_2}=$$

$$x_2{x}_5(k_{11}+k_{12})-k_4{x}_2{x}_4$$

$$\dot{x_3}=$$

$$\dot{x_4}=$$

$$\dot{x_5}=$$

$$x_2{x}_5(k_{11}+k_{12})-k_7{x}_3{x}_5-k_{16}{x}_1{x}_5$$

$$\dot{x_6}=$$

$$\dot{x_7}=$$

$$\dot{x_8}=$$

$$\dot{x_9}=$$

$$\dot{x_{10}}=$$

$$\dot{x_{11}}=$$

$$k_1$$

$$k_5$$

$$k_9$$

$$k_{13}$$

$$x_5+x_8+x_9+x_{10}+x_{11}$$

$$x_4+x_6+x_7$$

$$x_1+x_2+x_3+x_6+x_7+x_8+x_9+x_{10}+x_{11}$$

$$k_{17}$$

$$F\subset \mathbb{Z}[k_1,\dots ,k_{19}][x_1,\dots ,x_{11}].$$

$$P=\{k_1,\dots ,k_{19},x_1,\dots ,x_{11}\}\subset \mathbb{Z}[k_1,\dots ,k_{19}][x_1,\dots ,x_{11}]$$

$$\phi =\underset{f\in F}{\bigwedge }f=0\wedge \underset{v\in P}{\bigwedge }v>0.$$

$${\phi }_{500}^{(i)}=\mathrm{VS}(\exists x_1\dots \exists x_{i-1}\exists x_{i+1}\dots \exists x_{11}{\phi }_{500}).$$

$$3\le |S_{500}|\text{and}{S}_{500}\subseteq {\prod }_{i=1}^{11}\{a_i,b_i,c_i\}.$$

$${\phi }_{k_{19}}^{(i)}=\mathrm{VS}(\exists x_1\dots \exists x_{i-1}\exists x_{i+1}\dots \exists x_{11}{\phi }_{k_{19}}).$$

$${\sum }_{i=0}^{10}{c}_i{k}_{19}^i\text{with integer coefficients ci as in Table refTAB:coeffs}.$$

$$f(x_1,k_{19})={\sum }_{i=0}^6{d}_i{x}_1^i=0$$

Full text

1 Introduction

The occurrence of multiple steady states (multistationarity) has important consequences on the capacity of signaling pathways to process biological signals, even in its elementary form of two stable steady states (bistability). Bistable switches can act as memory circuits storing the information needed for later stages of processing [26]. The response of bistable signaling pathways shows hysteresis, namely dynamic and static lags between input and output. Because of hysteresis one can have, at the same time, a sharp binary response and protection against chatter noise.

Bistability of signaling usually occurs as a result of activation of upstream signaling proteins by downstream components [3]. A different mechanism for producing bistability in signaling pathways was proposed by Markevich et al. [19]. In this mechanism bistability can be caused by multiple phosphorylation/dephosphorylation cycles that share enzymes. A simple, two-step phosphorylation/dephosphorylation cycle is capable of ultrasensitivity, a form of all or nothing response with no hysteresis (Goldbeter–Koshland mechanism). In multiple phosphorylation/dephosphorylation cycles, enzyme sharing provides competitive interactions and positive feedback that ultimately leads to bistability.

Algorithmically the task is to find the positive real solutions of a parameterized system of polynomial or rational systems, since the dynamics of the network is given by polynomial systems (arising from mass action kinetics) or rational functions (arising in signaling networks when some intermediates of the reaction mechanisms are reduced). Due to the high computational complexity of this task [13] considerable work has been done to use specific properties of networks and to investigate the potential of multistationarity of a biological network out of the network structure. This only determines whether or not there exist rate constants allowing multiple steady states, instead of coming up with a semi-algebraic description of the range of parameters yielding this property. These approaches can be traced back to the origins of Feinberg’s Chemical Reaction Network Theory (CRNT) whose main result is that networks of deficiency 0 have a unique positive steady state for all rate constants [12, 8]. We refer to [6, 20, 16] for the use of CRNT and other graph theoretic methods to determine potential existence of multiple positive steady states, with [17] giving a survey.

Given a bistable mechanism it is also important to compute the bistability domains in parameter space, namely the parameter values for which there is more than one stable steady state. The size of bistability domains gives the spread of the hysteresis and quantifies the robustness of the switches. The work of Wang and Xia [24] is relevant here: they used symbolic computation tools, including cylindrical algebraic decomposition as we do below, to determine the number of steady states and their stability for several systems. They reported results up to a 5-dimensional system using specified parameter values, but their method is extensible to parametric questions. Higher-dimensional systems were studied using sign conditions on the coefficients of the characteristic polynomial of the Jacobian. In some cases these guarantee uniqueness of the steady state [7].

In this paper we use an 11-dimensional model of a mitogen-activated protein kinases (MAPK) cascade [19] as a case study to investigate properties of the system using algorithmic methods towards the goal of semi-algebraic descriptions of parameter regions for which multiple positive steady states exist.

2 The MAPK Network

The model of the MAPK cascade we are investigating can be found in the Biomodels database [18].111www.ebi.ac.uk/biomodels-main/BIOMD0000000026 We have renamed the species names to $x_{1}$, …, $x_{11}$ and the rate constants to $k_{1}$, …, $k_{16}$ to facilitate reading:

[TABLE]

The Biomodels database also gives us meaningful values for the rate constants:

[TABLE]

Some of these values are measured and some are well-educated guesses. For the purpose of our study we assume they are suitable.

We add three linear conservation constraints introducing three further constant parameters $k_{17}$, $k_{18}$, $k_{19}$:

[TABLE]

Computations to produce these in MathWorks SimBiology use the left-null space of the stoichiometric matrix under positivity conditions, see for example [21].

Meaningful values for $k_{17},k_{18},k_{19}$ are harder to obtain than the constants in (2). The following are some realistic values estimated by ourselves on the basis of our understanding of the biological model:

[TABLE]

The long-term goal of our research is to treat all three of these together parametrically, although in the present work we focus on situations with one free-parameter.

The steady state problem for the MAPK cascade can now be formulated as a real algebraic problem. That is, we replace the left hand sides of all equations in (1) with [math]. This together with the equations in (3) yields an algebraic system with polynomials in

[TABLE]

All entities in our model are strictly positive, which yields in addition a system

[TABLE]

establishing a side condition on the solutions of $F$ that all variables $x_{i}$ and parameters $k_{i}$ of $P$ be positive. In terms of first-order logic our specification of $F$ and $P$ yields a quantifier-free Tarski formula

[TABLE]

The estimations for the rate constants in (2) formally establish a substitution rule $\sigma=[0.02/k_{1},\dots,0.0011/k_{16}]$, which can be applied to $F$, $P$, or $\varphi$ in postfix notation.

2.1 Symbolic Determination of Occurrences of Multiple Steady States

In this section we are going to analyze the system for multiple positive steady states. As we will not include a priori information about the stability of the fixed points, we do not only have to consider (at least) two stable fixed points but also unstable fixed points, i.e., we investigate the existence of at least three different roots $\mathbf{x}\in\left]0,\infty\right[^{11}$ of $F$ for given choices $\mathbf{k}\in\left]0,\infty\right[^{19}$ of parameters.

We present two investigations: one using the Redlog package in Reduce and the other using the Regular Chains Library in Maple. Both will make use of Cylindrical Algebraic Decomposition (CAD) [1] to solve the problem. The worst-case time complexity of CAD is doubly exponential.222Traditionally, doubly exponential in the number of variables. However recent progress on CAD in the presence of equational constraints [10], such as (1) with [math] for left-hand side, allows us to conclude it is actually doubly-exponential the number of variables minus the number of equational constraints at different levels of the projection [11]. Our approaches admit, in principle, arbitrary numbers of indeterminates. However, for the sake of realistic computation times we must restrict ourselves to one free parameter. Even then, the number of variables present is too large for contemporary CAD implementations. We make progress by combining CAD with additional symbolic methods. Our first approach uses virtual substitution techniques and the second real triangularization. In both cases we have combined the corresponding methods by hand, but automation is clearly possible.

2.1.1 Real Quantifier Elimination in Redlog

Real Quantifier Elimination (QE) can directly handle the parametric existence of steady states, taking as input $\exists x_{1}\dots\exists x_{11}\varphi$, possibly with substitutions for some parameters. However, we are not only interested in the existence but also in the number of solutions. We are going to combine Virtual Substitution (VS) [25] with CAD. The former smoothly eliminates the majority of the quantifiers while the latter allows us to count numbers of solutions via decomposition of the remaining low-dimensional spaces. That combination of methods requires the solution of several QE runs with each problem and some combinatorial arguments. Throughout this subsection we are using Redlog [9].

Parameter-Free Computations

We consider $\varphi_{500}=\varphi\sigma[100/k_{17},50/k_{18},500/k_{19}]$ where all parameters have been substituted with rational numbers. The closed formula $\bar{\varphi}_{500}=\exists x_{1}\dots\exists x_{11}\varphi_{500}$ states the existence of a suitable real solution. In a first step, we solve for $i\in\{1,\dots,11\}$ the following eleven QE problems using VS:

[TABLE]

Each $\varphi_{500}^{(i)}$ is a univariate quantifier-free formula describing all possible real choices for $x_{i}$ for which there exist real choices for all other variables such that $\varphi_{500}$ holds. CAD can easily decompose the corresponding one-dimensional spaces. It turns out that for each $x_{i}$ there are exactly three zero-dimensional cells $a_{i}$, $b_{i}$, $c_{i}\in\mathbb{R}$ for which $\varphi_{500}^{(i)}$ holds. We extract all $a_{i}$, $b_{i}$, and $c_{i}$ as real algebraic numbers, i.e., univariate defining polynomials with integer coefficients plus isolating intervals. By combinatorial arguments it is not hard to see that the following holds for the set $S_{500}$ of real solutions of $\varphi_{500}$:

[TABLE]

Notice that at this point we have proven multistationarity for $k_{19}=500$. We can furthermore compute $S_{500}$ by plugging the $3^{11}$ candidates from the Cartesian product into $\varphi_{500}$. A straightforward approach requires arithmetic with real algebraic numbers followed by the determination of the signs of the results, which is quite inefficient in practice. We use instead a heuristic approach combining refinements of the isolating intervals of the real algebraic numbers with interval arithmetic. This excludes $3^{11}-3$ of the candidate solutions. The three remaining candidates require no further checking since we already know that $|S_{500}|\geq 3$. The overall CPU time is 71.3 seconds for 11 runs of VS plus 11 runs of CAD, followed by 16 hours for checking candidates.333All QE-related computations have been carried out on a 2.4 GHz Intel Core i7 with 3 GB RAM or cores on a compute server with similar speed and memory limitations. Our checking procedure is a file-based prototype starting a Reduce process for every single of the $3^{11}$ candidates; there is considerable room for optimization.

For $k_{19}=200$ instead of $500$ all eleven univariate CAD computations yield unique solutions which can be straightforwardly combined to one unique solution for the corresponding $\varphi_{200}$. The overall CPU time here is 66.4 seconds for 11 runs of VS plus 11 runs of CAD. Machine float approximations of all our solutions are given in Table 1.

Parametric Analysis for $k_{19}$

We now consider $\varphi_{k_{19}}=\varphi\sigma[100/k_{17},50/k_{18}]$ leaving $k_{19}$ as a parameter. Again, we solve for $i\in\{1,\dots,11\}$ eleven QE problems using VS:

[TABLE]

This time each $\varphi_{k_{19}}^{(i)}$ is a bivariate quantifier-free formula in $k_{19}$ and the corresponding $x_{i}$. This time we construct a two-dimensional CAD for each $\varphi_{k_{19}}^{(i)}$. The projection order is important: we first project $x_{i}$, then the CAD base phase decomposes the $k_{19}$-axis, followed by an extension phase that decomposes the $x_{i}$-space over the $k_{19}$-cells obtained in the base phase. This is feasible with one limitation: we do not extend over zero-dimensional $k_{19}$-cells. In other words, we accept finitely many blind spots in parameter space, which we can explicitly read off from the CAD so that in the end we know exactly what we are missing.

Figure 1 shows our CAD tree for $\varphi_{k_{19}}^{(2)}$. The first layer next to the root shows the decomposition of the $k_{19}$-axis. The five zero-dimensional (rectangular) cells are the previously mentioned blind spots, among which the smallest one with negative value of $k_{19}$ is not relevant. Those zero-dimensional cells also establish the limits of the full dimensional (oval) cells in between. The cylinders over those one-dimensional $k_{19}$-cells each contain either one or three zero-dimensional $x_{2}$-cells where $\varphi_{k_{19}}^{(2)}$ holds. We have deleted from the tree all $x_{2}$-cells where $\varphi_{k_{19}}^{(2)}$ does not hold. We make two observations, important for a qualitative analysis of our system:

1. (i)

For all positive choices of $k_{19}$—extending to infinity—there is at least one positive solution for $x_{2}$. 2. (ii)

There is a break point around $k_{19}=409.253$ where the system changes from unique solutions to exactly three solutions.

Recall that for all floating point numbers given here as approximations we in fact know exact real algebraic numbers. For instance, the exact break point is the only real zero in the interval $\mathopen{(}409,410\mathclose{)}$ of an irreducible defining polynomial

[TABLE]

Figure 2 depicts all eleven CAD trees for $\psi_{k_{19}}^{(1)}$, …, $\psi_{k_{19}}^{(11)}$.

They are quite similar to the one just discussed. Even the break point from one to three solutions for $x_{i}$ is identical for all $i\in\{1,\dots,11\}$ so that we can generalize our observations:

1. (i)

For all positive choices of $k_{19}$—extending to infinity—there is at least one positive solution for $(x_{1},\dots,x_{11})$. 2. (ii)

There is a break point $\beta$ around $k_{19}=409.253$ where the system changes its qualitative behavior. We have exactly given $\beta$ as a real algebraic number in Equation (6). For $k_{19}<\beta$ there is exactly one positive solution for $(x_{1},\dots,x_{11})$. For $k_{19}>\beta$ there are at least $3$ and at most $3^{11}$ positive solutions for $(x_{1},\dots,x_{11})$.

The overall computation time for our parametric analysis is 4.3 minutes. It is strongly dominated by 2.8 minutes for the computation of one particular CAD tree, for $\varphi_{k_{19}}^{(11)}$. It turns out that the suitable projection order with $x_{i}$ eliminated first is computationally considerably harder than projecting the other way round. As a preprocessing step we apply CAD-based simplification of the $\varphi_{k_{19}}^{(i)}$ with the opposite, faster, projection order. Here we use Qepcad B, which performs better than Redlog at simple solution formula construction.

2.1.2 Triangular Decomposition methods with the Regular Chains Library

We now describe an alternative approach to the solution using regular chains methods. Regular chains are the triangular decompositions of systems of polynomial equations (triangular in terms of the variables in each polynomial). Highly efficient methods for working in complex space have been developed based on these; see [23] for a survey.

Recent work by Chen et al. [4] proposes adaptations of these tools to the real analogue: semi-algebraic systems. They describe two algorithms to decompose any real polynomial system into finitely many regular semi-algebraic systems. The first does so directly while the second, Lazy Real Triangularize (LRT) produces the highest dimension solution component and unevaluated function calls, which if all evaluated would combine to give the full solution. These algorithms are implemented in the Regular Chains Library444www.regularchains.org in Maple which we use throughout this subsection.

We apply LRT on the quantifier-free formula (5) evaluated with the parameter estimates for $k_{1}$, …, $k_{18}$ given at the start of Section 2, so we have one free parameter as in the previous section. We need to choose a variable ordering: our analysis requires that $k_{19}$ be the indeterminate considered alone; the remaining variables are placed in lexicographical order (the in-built heuristics to make the choice could suggest nothing better). The solutions must hence contain constraints in $k_{19}$, constraints in ($x_{1},k_{19})$, in ($x_{2},x_{1},k_{19})$ and so on. We define the main variable of a constraint to be the highest one present in this ordering.

LRT produces one solution component and 6 unevaluated function calls in less than 3 seconds. In the evaluated component: for each of $x_{2}$, …, $x_{11}$ there is a single equation which had this as the main variable. Further, these are all linear in their main variable meaning they can be easily rearranged into the solution formulae in Table 3.

The constraints on $(x_{1},k_{19})$ are that $x_{1}>0$ and that a polynomial equation of degree $6$ be satisfied:

[TABLE]

where the coefficients $d_{i}$ are univariate polynomials in $k_{19}$ of maximum degree $2$ as given in Table 2.

Finally, the constraints on $k_{19}$ are that it be positive; it not be a root of the polynomial in Equation (6); nor two other polynomials as described in Table 4.

Thus this solution component is valid for all positive values of $k_{19}$ excluding three points. As before, we could give these as exact algebraic numbers but for brevity give float approximations: $409.253$, $16473.337$, and $25084.536$.

Three of the six unevaluated function calls define the solutions at these points, however evaluating these solutions is not possible in reasonable time. The other three define empty solution sets (evaluating to discover this is instantaneous). So, as with our previous approach, we proceed accepting a small number of blind spots.

The output of LRT has quickly given us the structure of the solution space valid at all but three isolated values of $k_{19}$. However, it does not identify where the number of real solutions change: although the break point identified earlier has been rediscovered there is no information from which we can infer its significance; and there is no significance in our application of the other two isolated points.

To finish the analysis we need to decompose $(x_{1},k_{19})$-space according to the real roots of $f(x_{1},k_{19})$; and also $x_{1}$ since the constraint $x_{1}>0$ was specified separately in the output (the case for this variable only). CAD is ideally suited for this task. Using the Regular Chains algorithm [5] in Maple a CAD for $f(x_{1},k_{19})$ divides the plane into 135 cells in a few seconds. This CAD decomposes the $k_{19}$ axis into 11 cells, i.e. identifying five points which approximate to: $-379.993$, $-87.776$, [math], $409.253$, and $25084.536$.

On the cell for $k_{19}\in\mathopen{]}0,409.253\mathclose{[}$, the cylinder above in the $(x_{1},k_{19})$ plane is divided into 11 cells: three of which cover $x_{1}>0$ (two 2d sectors and a 1d section). This indicates that $f(x_{1},k_{19})$ has a single positive real solution for such $k_{19}$. On the two cells for $k_{19}\in\mathopen{]}409.253,25084.536\mathclose{[}$ and $k_{19}\in\mathopen{]}25084.536,\infty\mathclose{[}$ the cylinders above are divided into 15 cells; seven of which cover $x_{1}>0$. This indicates that $f(x_{1},k_{19})$ has three positive real solutions for such $k_{19}$.

At the end of this analysis we have rediscovered the break point where the system moves from a single positive real solution to three. We also have explicit solutions valid for all except three isolated $k_{19}$ values. To obtain a solution select the $k_{19}$ value of interest then identify the real roots of $f(x_{1},k_{19})$ (we know in advance how many depending on the $k_{19}$ value chosen); then for each $x_{1}$ solution substitute recursively into the equations of Table 3; starting from the bottom and including the new variable solution discovered from each substitution into the next. The solutions in Table 1 may be easily rediscovered this way.

Repeating the Process for Different Choices of the Lone Free Parameter/Fixed Parameter Values

We may repeat the approach described above for different choices of free parameter and different choices of fixed parameter values. For example:

• •

With $k_{17}$ set to $95$ instead of $100$ we find that the break point between 1 and 3 real positive solutions moves to $k_{19}=369.917$. With $k_{17}$ set to $105$ it moves to $k_{19}=450.077$.

• •

Allowing $k_{17}$ to be free and fixing $k_{19}=200$ we find that there is only ever one positive real solution.

• •

Allowing $k_{17}$ to be free and fixing $k_{19}=500$ we find the number of positive real solutions moving from 1 to 3 to 1 breaking at $k_{17}=85.988$ and $k_{17}=110.869$.

• •

Similarly, allowing $k_{18}$ to be free and fixing $k_{19}=200$ we find there is only ever one positive real solution; but fixing $k_{19}=500$ instead we find 3 real solutions between $k_{18}=51.382$ and $58.329$ and 1 otherwise.

The results above hint that there is a shape approximating a narrow paraboloid in $(k_{17},k_{18},k_{19})$-space within which bistability may occur; with bistability available for any $k_{17}$ and $k_{18}$ value but bounded from below in the $k_{19}$ coordinate. We note that these additional experiments all produce results which, as with the one described in detail, are invalid at a handful of isolated values of the free parameter.

2.2 Stability of the Fixed Points

We use the three linear conservation constraint equations (3) to eliminate $x_{1}$, $x_{7}$, and $x_{11}$ from system (1) and symbolically compute the Jacobian $\tilde{J}$ of the obtained reduced system. We then numerically compute the eigenvalues of $\tilde{J}$ for the instances arising from the substitution of the different positive fixed points for the variables and the corresponding parameter values.

We have used the float approximations for the unique solution $x^{(200)}$ with $k_{19}=200$ and the three solutions $x^{(500)}_{1}$, …, $x^{(500)}_{3}$ for $k_{19}=500$ in Table 1. For the single positive fixed point $x^{(200)}$ the Jacobian $\tilde{J}(x^{(200)})$ has eigenvalues with negative real part only and hence can be shown to be stable. For $k_{19}=500$ one of the three positive fixed points $x^{(500)}_{2}$ can be shown to be unstable, as $\tilde{J}(x_{2}^{(500)})$ has one eigenvalue with positive real part; the other seven had negative real parts. In contrast $x_{1}^{(500)}$ and $x_{3}^{(500)}$ can be shown to be stable. Hence for $k_{19}=500$ the system is indeed bistable.

A verification of the stability of the fixed points using exact real algebraic numbers by the well-known Routh–Hurwitz criterion is possible algorithmically [15], but seems to be out of range of current methods for this example. Notice that also in other studies on multistationarity of signalling pathways [6, 14] the question of stability has not been addressed.

2.3 Numerical Homotopy Methods

Finally, we compare our symbolic results with numerical ones obtained using the homotopy solver Bertini [2]. Bertini computes complex roots of polynomial systems using methods from numerical algebraic geometry [22].

For the parameter values as above and $k_{19}=500$ we obtain six solutions, three of which are positive real solutions. For $k_{19}=200$ we obtain a single positive solution. In both cases the relevant solutions coincide with the ones obtained with our symbolic analyses up to the used numeric precision.

However, for larger values of $k_{19}$ Bertini produces incorrect results due to numerical instability. For instance, we falsely obtain exactly one positive real solution for $k_{19}=6000$ and no positive real solution for $k_{19}=10000$.

Figure 3 shows a Bertini-based grid sampling of parameter regions, varying $k_{19}$ between 200 and 1000 and fixing one of $k_{17}$ and $k_{18}$ while varying the other among the default values (4). While this suffers from the discrete nature of sampling and potentially unreliable results as discussed, it is nevertheless useful for the generation of hypothesis about the nature of the parameter regions. Figure 3 seems to identify a region of bistability (in blue) within the parameter space, as hypothesised at the end of Section 2.1.2. The results of Bertini indicate holes in this region (the green dots within the blue). However, computation at these particular points reveals these to be the result of numerical errors: where an insufficiently high precision causes what is actually a positive real solution to appear to have a negative component. It seems there is scope for fruitful interplay between symbolic and numeric methods here; with numerics postulating hypotheses for the symbolic methods to check and refine.

3 Conclusions and Future Work

We have shown that the determination of multistationarity of an 11-dimensional MAPK network can be achieved by combinations of currently available symbolic computation methods for mixed equality/inequality systems if, for all but potentially one parameter, numeric values are known. The aspiration of a semi-algebraic description of the ranges for all parameters in the conservation laws (3) yielding multistationarity will now be pursued, with the present results demonstrating that this aspiration may be within reach.

As there are many very relevant systems having dimensions between 10 and 20 it seems to be worth the effort to enhance and improve the present algorithmic methods, and in particular their combination, to solve such important application problems for symbolic computation.

Acknowledgments

For our QE-related computations we used two great free software tools: GNU Parallel for distributing computations on several processors, and yEd for visualization of CAD trees. D. Grigoriev is grateful to the grant RSF 16-11-10075 and to MCCME for wonderful working conditions and inspiring atmosphere. M. Košta has been supported by the DFG/ANR Project STU 483/2-1 SMArT. H. Errami, A. Weber, and O. Radulescu thank the French-German Procope-DAAD program for partial support of this research. J. Davenport, M. England and T. Sturm are grateful to EU CSA project 712689 SC2.