Parameter regions that give rise to 2[n/2]+1 positive steady states in the n-site phosphorylation system

The distributive sequential n-site phosphorylation/dephosphorylation system is an important building block in networks of chemical reactions arising in molecular biology, which has been intensively studied. In the nice paper of Wang and Sontag (2008) it is shown that for certain choices of the reaction rate constants and total conservation constants, the system can have 2[n/2]+1 positive steady states (that is, n+1 positive steady states for n even and n positive steady states for n odd). In this paper we give open parameter regions in the space of reaction rate constants and total conservation constants that ensure these number of positive steady states, while assuming in the modeling that roughly only 1/4 of the intermediates occur in the reaction mechanism. This result is based on the general framework developed by Bihan, Dickenstein, and Giaroli (2018), which can be applied to other networks. We also describe how to implement these tools to search for multistationarity regions in a computer algebra system and present some computer aided results.


Introduction
Multisite phosphorylation cycles are common in cell regulation systems [18]. The important distributive sequential n-site phosphorylation network is given by the following reaction scheme: Arrows correspond to chemical reactions. The labels in the arrows are positive numbers, known as reaction rate constants. This reaction scheme represents one substrate S 0 that can sequentially acquire up to n phosphate groups (giving rise to the phosphoforms S 1 , . . . , S n ) via the action of the kinase E, and which can be sequentially released via the action of the phosphatase F , in both cases via an intermediate species (denoted by Y 0 , Y 1 , . . . , Y n−1 , U 0 , U 1 , . . . , U n−1 ) formed by the union of the substrate and the enzyme . The kinetics of this network is deduced by applying the law of mass action to the labeled digraph of reactions (1.1), which yields an autonomous polynomial system of ordinary differential equations depending on the positive reaction rate constants. This system describes the evolution in time of the concentrations of the different chemical species. We denote with lower case letters the concentrations of the 3n + 3 species. The derivatives of the concentrations of the substrates satisfy (see e.g. [7] for the general definition): and de dt = − dy 0 dt − · · · − dy n−1 dt , df dt = − du 0 dt − · · · − dun dt , which give two linear conservation relations. Indeed, there are three linearly independent conservation laws: where the linear conservation constants S tot , E tot , F tot are positive for any trajectory of the system intersecting the positive orthant.
A steady state of the system is a common zero of the polynomials in the right-hand side of the differential equations. There are 6n + 3 parameters: the reaction rate constants and the linear conservation constants (which correspond to points in the positive orthant R 6n+3 >0 ). It is known that for any n ≥ 2, the sequential n-site system shows multistationarity. This means that there exists a choice of parameters for which there are at least two positive steady states satisfying moreover the linear conservation relations (1.2) for the same S tot , E tot and F tot . In this case, it is said that the steady states are stoichiometrically compatible, or that they lie in the same stoichiometric compatibility class. This is an important notion, because the occurrence of multistationarity allows for different cell responses with the same linear conservation constants (that is, the same total amounts of substrate and enzymes).
The possible number of positive steady states of the n-site phosphorylation system (for fixed linear conservation constants) has been studied in several articles. For n = 2, it is a well known fact that the number of nondegenerate positive steady states is one or three [20,24]. The existence of bistability is proved in [15]. In [3] and [5], the authors give conditions on the reaction rate constants to guarantee the existence of three positive steady states based on tools from degree theory, but this approach does not describe the linear conservation constants for which there is multistationarity. For an arbitrary number n of phosphorylation sites, it was shown in [24] that the system has at most 2n − 1 positive steady states. In the same article, the authors showed that there exist reaction rate constants and linear conservation constants such that the network has n (resp. n + 1) positive steady states for n odd (resp. even); that is, there are 2[ n 2 ] + 1 positive steady states for any value of n, where [.] denotes integer part.
In [13] (see also [16]) the authors showed parameter values such that for n = 3 the system has five positive steady states, and for n = 4 the system has seven steady states, obtaining the upper bound given in [24]. In the recent article [4] the authors show that if the n-site phosphorylation system is multistationary for a choice of rate constants and linear conservation constants (S tot , E tot , F tot ) then for any positive real number c there are rate constants for which the system is multistationary when the linear conservation constants are scaled by c. Concerning the stability, in [23] it is shown evidence that the system can have 2[ n 2 ] + 1 positive steady states with [ n 2 ] + 1 of them being stable. Recently, a proof of this unlimited multistability was presented in [10], where the authors find a choice of parameters that gives the result for a smaller system, and then extend this result using techniques from singular perturbation theory.
In the previous paper [1], parameter regions on all the parameters are given for the occurrence of multistationarity for the n-site sequential phosphorylation system, but no more than three positive steady states are ensured. These conditions are based on a general framework to obtain multistationary regions jointly in the reaction rate constants and the linear conservation constants. Our approach in this article uses the systematic technique in [1], which we briefly summarize in Section 3.
The removal of intermediates was studied in [12]. More specifically, the emergence of multistationarity of the n-site phosphorylation system with less intermediates was studied in [11]. It is known that the n-site phosphorylation network without any intermediates complexes has only one steady state for any choice of parameters. In [11], the authors showed which are the minimal sets of intermediates that give rise to a multistationarity system, but they do not give information about how many positive steady states can occur, and also, they do not describe the parameter regions for which these subnetworks are multistationary.
In this paper, we work with subnetworks of the sequential n-site phosphorylation system that only have intermediates in the E component (that is, in the connected component of the network where the kinase E reacts), see Definition 1.1. In case of the full mechanism on the E component or if we only assume that there are intermediate species that are formed between the phosphorylated substrates with parity equal to n (that is, roughly only 1 4 of the intermediates of the n sequential phosphorylation cycle), we obtain precise conditions on all the parameters to ensure as many positive steady states as possible. Indeed, we show in Proposition 2.2 that the maximum number of complex solutions to the steady state equations intersected with the linear conservation relations is always n + 1, the maximum number of real roots is also n + 1, that could be all positive when n is even, while only n of them can be positive when n is odd, so the maximum number of positive steady states equals 2[ n 2 ] + 1 for any n. Exact conditions on the parameters so that the associated phosphorylation/dephosphorylation system admits these number of positive steady states, are given in Theorem 1.2 and Corollary 1.5. The latter follows from Theorem 1.4. In order to state these results, we need to introduce some notations. Definition 1.1. For any natural number n, we write I n = {0, . . . , n − 1}. Given n ≥ 1, and a subset J ⊂ I n , we denote by G J the network whose only intermediate complexes are Y j for j ∈ J, and none of the U i . It is given by the following reactions: where the labels of the arrows are positive numbers. We will also denote by G J the associated differential system with mass-action kinetics.
For all these systems G J , there are always three linearly independent conservation laws for any value of n: where the total conservation constants S tot , E tot , F tot are positive for any trajectory of the differential system which intersects the positive orthant. Note that the concentration of the phosphatase F is constant, equal to F tot .
To get lower bounds on the number of positive steady states with fixed positive linear conservation constants, we first consider the network G In , that is, when all the intermediates in the E component appear. It has the following associated digraph: With the previous notations, consider the network G In in (1.5), and suppose that the reaction rate constants k cat i and ν i , i = 0, . . . , n − 1, satisfy the inequality For any positive values S tot , E tot and F tot of the linear conservation constants with verifying the inequalities: there exist positive constants B 1 , . . . , B n such that for any choice of positive constants λ 0 , . . . , λ n−1 satisfying rescaling of the given parameters k on j by λ j k on j , for each j = 0, . . . , n − 1, gives rise to a system with exactly 2[ n 2 ] + 1 nondegenerate positive steady states. Remark 1.3. We will also show in the proof of Theorem 1.2, that for any reaction rate constants and linear conservation constants satisfying (1.6), there exist t 0 > 0 such that for any value of t ∈ (0, t 0 ), the system G In has exactly 2[ n 2 ] + 1 nondegenerate positive steady states after modifying the constants k on j by t j−n k on j for each j = 0, . . . , n − 1.
We now consider subnetworks G J , with J ⊂ J n where (1.8) J n := {i ∈ I n : (−1) i+n = 1}, for n ≥ 1, that is, subsets J with indexes that have the same parity as n.
Theorem 1.4. Let n ≥ 1, and consider a subset J ⊂ J n . Let G J be its associated system as in (1.3). Assume moreover that A multistationarity region in the space of all parameters for which the system G J admits at least 1 + 2|J| positive steady states can be described as follows. Given any positive value of F tot , choose any positive real numbers k cat j , ν j , with j ∈ J satisfying Then, there exist positive constants B 1 , . . . , B n such that for any choice of positive constants λ 0 , . . . , λ n−1 satisfying rescaling of the given parameters k on j by λ j k on j , for j ∈ J and τ j by λ j τ j if j / ∈ J gives rise to a system with at least 1 + 2|J| positive steady states.
The following immediate Corollary of Theorem 1.4 implies that we can obtain a region in parameters space with [ n 2 ] intermediates, where the associated system has 2[ n 2 ] + 1 positive steady states. Assume moreover that Then, there is a multistationarity region in the space of all parameters for which the network G Jn admits 2[ n 2 ] + 1 steady states (with fixed linear conservation constants corresponding to the coordinates of a vector of parameters in this region), described in the statement of Theorem 1.4.
We explain in Section 5 the computational approach to find new regions of multistationarity. This is derived from the framework in [1] that we used to get the previous results. We find new regions of multistationarity for the cases n = 2, 3, 4 and 5, after some manual organization of the automatic computations. Our computer aided results are summarized in Propositions 5.2, 5.3, 5.4, 5.5, 5.6 and 5.7.
The paper is organized as follows. In Section 2 we prove Proposition 2.2 and we show how to lift regions of multistationarity from the reduced system G J to the complete sequential n-site phosphorylation system. We also show in Lemma 2.3 that even with a single intermediate Y 0 it is possible to make a choice of all parameters such that the system has 2[ n 2 ] + 1 positive steady states. This result has been independently found by Elisenda Feliu (personal communication). This says that while Corollary 1.5 is optimal, the regions obtained for any subset J with indexes of the same parity of n in Theorem 1.4 properly contained in J n , only ensure 2|J| + 1 positive steady states. However, note that we are able to describe open regions in parameter space and Lemma 2.3 only allows us to get choices of parameters.
In Section 3 we briefly recall the framework in [1], which is the basis of our approach. In Section 4 we give the proofs of Theorems 1.2 and 1.4. Finally, as we mentioned above, we explain in Section 5 how to implement the computational approach to find new regions of multistationarity. We end with a discussion where we further compare our detailed results with previous results in the literature.

Upper bounds and extension of multistationarity
In this section we collect three results that complete our approach to describe multistationarity regions giving lower bounds for the number of positive steady states with fixed linear conservation relations for the systems G J in Definition 1.1 for any J ⊂ I n . We first show a positive parametrization of the concentrations at steady state of the species of the systems G J , which allows us to translate the equations defining the steady states in all the concentrations plus the linear conservation relations into a system of two equations in two variables. In Proposition 2.2, we prove that any mass-action system G J , has at most 2[ n 2 ] + 1 positive steady states. Lemma 2.3 shows that having a single intermediate is enough to get that number of positive steady states, for particular choices of the reaction rate constants. However, this computation by reduction to the univariate case is not systematic as the general approach from [1] that we use to describe multistationarity regions in Theorems 1.2 and 1.4, which can be applied to study other quite different mechanisms. It is known that if a system G J has m nondegenerate positive steady states for a subset J ⊂ I n , then it is possible to find parameters for the whole n-site phosphorylation system that also give at least m positive steady states (see [12]). In Theorem 2.4, we give precise conditions on the rate constants to lift the regions of multistationarity for the reduced networks to regions of multitationarity with 2[ n 2 ] + 1 positive steady states (with fixed linear conservation constants) of the complete n-sequential phosporylation cycle.

2.1.
Parametrizing the steady states. The following lemma gives a positive parametrization of the concentration of the species at steady state for the systems G J , in terms of the concentrations of the unphosphorylated substrate S 0 and the kinase E. It is a direct application of the general procedure presented in Theorem 3.5 in [21], and generalizes the parametrization given in Section 4 of [1].
Lemma 2.1. Given n ≥ 1 and a subset J ⊂ I n , consider the system G J as in Definition 1.1. Denote for each j ∈ J set T −1 = 1, and for any i = 0, . . . , n − 1: Then, the parametrization of the concentrations of the species at a steady state in terms of s 0 and e is equal to: The inverses K −1 j of the constants K j in (2.1) in the statement of Lemma 2.1 are usually called Michaelis-Menten constants.
Let n ≥ 1 and a subset J ⊂ I n . If we substitute the monomial parametrization of the concentration of the species at steady state (2.3) into the conservation laws, we obtain a system of two equations in two variables s 0 and e. We have: We can write system (2.4) in matricial form: is the matrix of coefficients: Note that if we order the variables s 0 , e, the support of the system (that is, the exponents of the monomials that occur) is the following set A: independently of the choice of J ⊂ I n . We also recall the classical Descartes rule of signs. Descartes rule of signs: Let p(x) = c 0 + c 1 x + · · · + c m x m be a nonzero univariate real polynomial with r positive real roots counted with multiplicity. Denote by s the number of sign variations in the ordered sequence of the signs sign(c 0 ), . . . , sign(c m ) of the coefficients, i.e., discard the 0's in this sequence and then count the number of times two consecutive signs differ. Then r ≤ s and r and s have the same parity, which is even if c 0 c m > 0 and odd if c 0 c m < 0.

Upper bounds
We then have that 2[ n 2 ] + 1 is an upper bound for the number of positive real solutions of the system of equations defining the steady states of any system G J in Definition 1.1: For any choice of rate constants and total conservation constants, the dynamical system G J associated with any subset J ⊂ I n has at most 2[ n 2 ] + 1 isolated positive steady states. In fact, the polynomial system of equations defining the steady states of G J can have at most n + 1 isolated solutions in (C * ) d .
Proof. The number of positive steady states of the subnetwork G J is the number of positive solutions of the sparse system (2.4) of two equations and two variables. The support of the system is (2.8) whose convex hull is a triangle with Euclidean volume equal to n+1 2 . By Kushnirenko Theorem, the number of isolated solutions in (C * ) 2 is at most 2! (n+1) 2 = n+1. In particular, the number of isolated positive solutions is at most n + 1.
Moreover, when all positive solutions are nondegenerate, their number is necessarily odd by Corollary 2 in [3], which is based on the notion of Brouwer's degree. Indeed, in our case, we can bypass the condition of nondegeneracy because we can use Descartes rule of signs in one variable. In fact, from the first equation of system (2.4), we can write: where p(e) is the following polynomial of degree n on the variable e: . Note that for any e > 0 it holds that p(e) > 0, and so s 0 > 0. If we replace (2.9) in the second equation of (2.4), we have: The number of positive solutions of the equation (2.12) is the same if we multiply by p(e): This last polynomial q has degree n + 1, with leading coefficient equal to α n−1 + β n−1 > 0 and constant coefficient equal to −E tot < 0. The sign variation of the coefficients of q has the same parity as the sign variation of the leading coefficient and the constant coefficient, which is one. So, by Descartes rule of signs, as the sign variation is odd, the number of positive solutions is also odd.

2.3.
One intermediate is enough in order to obtain 2[ n 2 ] + 1 positive steady states. As we mentioned in the introduction, the following result has been independently found by Elisenda Feliu (personal communication). Lemma 2.3. If J = {0}, then there exists parameter values such that the system G J admits 2[ n 2 ] + 1 positive steady states. Proof. Assume n is even, then n + 1 is odd. As we did in the proof of Proposition 2.2, the positive solutions of the system (2.4) are in bijection with the positive solutions of the polynomial q(e) in (2.13). Here β 0 = K 0 and β i = 0 for i = 0. We will consider the polynomialq(e) := q(e) Etot , with constant coefficient equal to −1. Consider any polynomial of degree n + 1 (2.14) c n+1 e n+1 + c n e n + · · · + c 1 e − 1, with n + 1 distinct positive roots, and with constant term equal to −1. Then, we have that c i (−1) i+1 > 0, and in particular, c n+1 > 0. Our goal is to find reaction rate constants and total conservation constants such that the polynomial (2.14) coincides with the polynomialq(e). Comparing the coefficient of e i , for i = 1, . . . , n + 1 in both polynomials, we need to have: Keep in mind that the values of c i are given. We may solve (2.15) in terms of E tot and of the c i , i = 1, . . . , n + 1 : Note that if we take any value for E tot > 0, then the values of α i for i = 1, . . . , n − 1, α 0 + K 0 and S tot K 0 are completely determined. So, we find an appropriate value of E tot such that the previous values α i , K 0 and S tot are all positive. For that, we choose K 0 = 1, and we take E tot large enough such that This is possible since c n+1 > 0 and that imposing the first condition just on k odd is enough to ensure that it holds for all k ∈ I n−1 as well because of the signs of the c i . With these conditions, and using the equations (2.16), the values of α i for each i = 0, . . . , n − 1 and S tot are determined and are all positive. Now, it remains to show that we can choose reaction rate constants such that the values of α i , i = 0, . . . , n − 1 are the given ones. Recall that Take for example F tot = 1, k on 0 = 2, k off 0 = 1 and k cat 0 = 1 (to obtain K 0 = 1). Then, τ 0 = 1, so we take ν 0 = 1 α 0 . As α i+1 = α i τ i+1 ν i+1 , for i = 0, . . . , n − 2, we can choose any positive values of τ i+1 , ν i+1 such that τ i+1 ν i+1 = α i+1 αi , and we are done. When n is odd, with a similar argument, we can find reaction rate constants and total conservation constants such that the polynomialq(e) gives a polynomial like (2.14) (but with n distinct positive roots and one negative root).

2.4.
Lifting regions of multistationarity. Multistationarity of any of the subsystems G J can be extended to the full n-site phosphorylation system (for instance, by Theorem 4 in [12]). We give a precise statement of this result in Theorem 2.4.
Consider the full n-site phosphorylation network (1.1), with a given vector of reaction rate constants κ ∈ R 6n : We define the following rational functions of κ: where µ j (κ) and η j (κ) are in turn the following rational functions: We denote by ϕ : the function that takes κ and gives a vector of (positive) reaction rate constants with the following order: first, the constants k on j , k off j , k cat j , j ∈ J, then τ (κ), j / ∈ J, and then ν j (κ), j = 0, . . . , n − 1.
Given a subset J ⊂ I n and a vector of reaction rate constants κ ∈ R 6n >0 , we consider the subnetwork G ϕ(κ) J as in Definition 1.1, with rate constants ϕ(κ): Applying Theorem 6.4 from [8], which is built on Theorem 4 in [12], we get the following lifting result. This last result allows us to obtain multistationary regions for the complete n-site phosphorylation system, combining the conditions on the parameters given in Theorem 1.2 and Theorem 1.4, with conditions (2.20) of Theorem 2.4. In particular, let J n ⊂ I n as in (1.8). By lifting a multistationarity region for the system G Jn in Corollary 1.5, we get a multistationarity region of parameters of the n-site phosphorylation cycle with 2[ n 2 ] + 1 positive steady states in the same stoichiometric compatibility class.

Positive solutions of sparse polynomial systems
As we saw in Subsection 2.1, the steady states of the systems G J correspond to the positive solutions of the sparse polynomial system (2.4) in two variables. In this section, we briefly recall the general setting from [1,2] to find lower bounds for sparse polynomial systems, that we will use in the proof of Theorems 1.2 and 1.4 in Section 4 and also in Section 5. For detailed examples of this approach, we refer the reader to Section 2 in [14].
We consider a fixed finite point configuration A = {a 1 , . . . , a n } ⊂ Z d , with n ≥ d + 2. A sparse polynomial system of d Laurent polynomial equations with support A is a system where the exponents belong to A. We call C = (c ij ) ∈ R d×n the coefficient matrix of the system and we assume that no column of C is identically zero. Recall that a zero of (3.1) is nondegenerate when it is not a zero of the Jacobian of f 1 , . . . , f d . Our method to obtain a lower bound on the number of positive steady states, based on [1,2], is to restrict our polynomial system (3.1) to subsystems which have a positive solution and then extend these solutions to the total system, via a deformation of the coefficients. The first step is then to find conditions in the coefficient matrix C that guarantee a positive solution to each of the subsystems. The hypothesis of having subsystems supported on simplices of a regular subdivision of A is the key to the existence of a on open set in the space of coefficients where all these solutions can be extended. We recall below only the concepts that we need for our work. A comprehensive treatment of this subject can be found in [6]. Following Section 3 in [2], we define: A d-simplex with vertices in A is a subset of d + 1 points of A which is affinely independent. By Proposition 3.3 in [2], a system of d polynomial equations in d variables with a d-simplex as support, has one non-degenerate positive solution if and only if its d × (d + 1) matrix of coefficients is positively spanning. We further define, following [2]: Given a fixed finite point configuration A, take a height function h : A → R, h = (h(a 1 ), . . . , h(a n )). Consider the lower convex hull L of the n lifted points (a j , h(a j )) ∈ R d+1 , j = 1, . . . , n (see Figure 1). Project to R d the subsets of points in each of the faces of L (that it, points where the affine linear function which defines the lower face is minimized). These subsets define a regular subdivision of A induced by h. When the height vector h is generic, the regular subdivision is a regular triangulation, in which all the subsets are simplices. Note that the set of all height vectors inducing a regular subdivision of A that contains certain d-simplices ∆ 1 , . . . , ∆ p is defined by a finite number of linear inequalities. Thus, this set is a finitely generated convex cone C ∆ 1 ,...,∆p in R n with apex at the origin. In particular, the set of all height vectors inducing a regular subdivision Γ of A is a finitely generated convex cone in R n , which we call C Γ . In particular, if the simplices ∆ 1 , . . . , ∆ p completely determine the regular subdivision Γ, the cone C ∆ 1 ,...,∆p is equal to the cone C Γ .
Given A ⊂ Z d , consider a system of sparse real polynomials f 1 = f 2 = · · · = f d = 0 as in (3.1) with support A. Let Γ be a regular subdivision of A and h any height function that induces Γ. We define the following family of real polynomial systems parametrized by a positive real number t > 0: We also consider the following family of polynomial systems parametrized by γ ∈ R n >0 : c ij γ j x a j = 0, i = 1, . . . , d.

Theorem 3.3 (Theorem 3.4 of [2] and Theorem 2.11 of [1]
). Consider ∆ 1 , . . . , ∆ p dsimplices which occur in a regular subdivision Γ of a finite configuration A ⊂ Z d , and which are positively decorated by a matrix C ∈ R d×n .
(1) Let h be any height function that defines Γ. Then, there exists a positive real number t 0 such that for all 0 < t < t 0 , the number of nondegenerate positive solutions of (3.2) is at least p. (2) Let m 1 . . . , m ∈ R n be vectors that define a presentation of the cone C ∆ 1 ,...,∆p : Then, for any ε ∈ (0, 1) there exists t 0 (ε) > 0 such that for any γ in the open set We start this section with a lemma.  Proof. We can take h : A → R, with h(0, 0) = 0, h(0, 1) = n and h(1, j) = j(j−1) 2 , for j = 0, . . . , n − 1. It is easy to check h defines a regular triangulation that is equal to T .
The idea in the proofs of Theorem 1.2 and Theorem 1.4 is to detect positively decorated simplices in the regular triangulation T .
Proof of Theorem 1.2. By Proposition 2.2, the number of positive solutions of the system G In is at most 2[ n 2 ] + 1. So, it is enough to prove that this number is also a lower bound. The number of positive steady states of the system G In is the number of positive solutions of the system (2.4). As we saw before, the support of this last system is A = {(1, 0), (0, 1), (1, 1), (1, 2), . . . , (1, n), (0, 0)} ⊂ Z 2 , with coefficient matrix C (2.6). Note that if one multiplies a column of C by a positive number, then a simplex is positively decorated by C if and only if it is positively decorated by the new matrix. After multiplying the columns by convenient positive numbers, we obtain the following matrix from C: where M i = kcat i ν i Ftot + 1, for each i = 0, . . . , n − 1. We will work with this new matrix Csimple.
We consider the regular triangulation T in  We can include the last simplex if and only if n is even (because otherwise the inequalities are not compatible), and in this case we have at least n + 1 positively decorated simplices.
We can obtain 2[ n 2 ] + 1 positively decorated simplices if the inequalities (4.1) are satisfied. These inequalities are equivalent to the inequalities (1.6) in the statament.
Assume (1.6) holds. Given any height function h inducing the triangulation T , by item (1) in Theorem 3.3 there exists t 0 in R >0 such that for all 0 < t < t 0 , the number of positive nondegenerate solutions of the deformed system as in (3.2) with support A and coefficient matrix C t , with (C t ) ij = t h(α j ) c ij (with α j ∈ A, C = (c ij )), is at least 2[ n 2 ] + 1. In particular if we choose h as in the proof of Lemma 4.1, there exists t 0 in R >0 , such that for all 0 < t < t 0 , the system has at least 2[ n 2 ] + 1 positive solutions. If we change the variableē = t n e, we get the following system: It is straightforward to check that if we scale the constants K j by while keeping fixed the values of the constants k cat j , ν j for j = 0, . . . , n−1 and the values of the linear conservation constants E tot , F tot and S tot (assuming that condition (1.6) holds), the intersection of the steady state variety and the linear conservation equations of the corresponding network is described by system (4.3). It is easy to check that in order to get the scaling in (4.4), it is sufficient to rescale only the original constants k on j as follows: t j−n k on j , for j = 0, . . . , n − 1. Then, for these choices of constants, the system has at least 2[ n 2 ] + 1 positive steady states. Now, we will show how to obtain the more general rescaling in the statement. The existence of the positive constants B 1 , . . . , B n follows from the inequalities that define the cone C T of heights inducing the regular triangulation T and item (2) in Theorem 3.3. For instance, we can check that C T is defined by n inequalities: where , denotes the canonical inner product of R n+3 and m 1 = e 1 − 2e 3 + e 4 , m j = e j+1 − 2e j+2 + e j+3 , for j = 2, . . . , n − 1 and m n = e 2 + e n+1 − e n+2 − e n+3 , where e i denotes the i-th canonical vector of R n+3 . Fix ε ∈ (0, 1) n+3 . As (1.6) holds, item (2) in Theorem 3.3 says that there exist positive numbers B j for j = 1, . . . , n (depending on ε), such that the system has at least 2[ n 2 ] + 1 nondegenerate positive solutions, for any vector γ ∈ R n+3 satisfying γ m j < B j , for all j = 1, . . . , n. In particular, this holds if we take γ 1 = γ 2 = γ n+3 = 1 and If we call λ 0 = γ 3 and λ j = γ j+3 γ j+2 for j = 1, . . . , n − 1, the inequalities in (4.6) are equivalent to the conditions (1.7). Then, if λ j , j = 0, . . . , n−1, satisfy these inequalities, the rescaling of the given parameters k on j by λ j k on j for j = 0, . . . , n − 1, gives rise to a system with exactly 2[ n 2 ] + 1 positive steady states. The proof of Theorem 1.4 is similar to the previous one.
Proof of Theorem 1.4. Again, the number of positive steady states of our system is equal to the number of positive solutions of the system (2.4). Recall that the support of the system is In this case, the coefficient matrix C (2.6) is equal, after multiplying the columns by convenient positive numbers, to the matrix where M i = kcat i ν i Ftot + 1 and D i = 1, for each i ∈ J, and M i = 1 and D i = 0, for each i / ∈ J. We consider again the regular triangulation T in Lemma 4.1. Recall that J ⊂ J n , see (1.8), and therefore each j ∈ J has the same parity as n, in particular 0 ≤ j ≤ n − 2. For each j ∈ J, consider the simplices ∆ j = {(1, j), (1, j + 1), (0, 0)} and ∆ j+1 = {(1, j + 1), (1, j + 2), (0, 0)}. Note that if j = j then {∆ j , ∆ j+1 } and {∆ j , ∆ j +1 } are disjoint since j, j and n have the same parity.
The simplices are positively decorated by Csimple (and then by C) if and only if the submatrices are positively spanning, and this happens if and only if E tot M j − S tot < 0, where M j = kcat j ν j Ftot + 1, since j ∈ J. The simplex ∆ n = {(0, 1), (1, n), (0, 0)} is trivially positively decorated by Csimple. Then, by imposing the inequalities E tot M j − S tot < 0 for j ∈ J, which are equivalent to the ones in the statement (1.9), we can obtain 2|J| + 1 positively decorated simplices. Assume (1.9) holds. Given any height function h inducing the triangulation T , by item (1) in Theorem 3.3 there exists t 0 in R >0 , such that for all 0 < t < t 0 , the number of positive nondegenerate solutions of the deformed system with support A and coefficient matrix C t , with (C t ) ij = t h(α j ) c ij (with α j ∈ A, C = (c ij )) is at least 2|J| + 1. In particular if we choose h as in the proof of Lemma 4.1, there exists t 0 in R >0 , such that for all 0 < t < t 0 , the system has at least 2|J| + 1 positive solutions. If we change the variableē = t n e, we get the following system: Similarly as we did in the previous proof, if we scale the original parameters k on j , for j ∈ J, respectively, and if we keep fixed the values of the remaining rate constants and the values of the linear conservation constants E tot , F tot and S tot , the intersection of the steady state variety and the linear conservation relations is described by system (4.8). Then, for these choices of constants the system G J has at least 2|J| + 1 positive steady states. The general rescaling that appears in the statement can be obtained in a similar way as we did in the proof of Theorem 1.2.

Computer aided results
In this section we explore a computational approach to the multistationarity problem, more precisely we find new regions of mulstistationarity. The idea is to find good regular triangulations of the point configuration corresponding to the exponents of the polynomial system given by the steady state equations and conservation laws, and deduce from them some regions of multistationarity. We first give the idea and then apply it for the n-site phosphorylation system for n = 2, 3, 4, and 5, where we have successfully found several regions of multistationarity. This approach can be, in principle, applied to other systems if they satisfy certain hypotheses, see [1,Theorem 5.4], and are sufficiently small in order for the computations to be done in a reasonable amount of time.
The strategy is the following. Given a polynomial system with support A and matrix of coefficients C, one first computes all possible regular triangulations of A with the aid of a computer. The number of such triangulations can be very large depending on the size of A, thus the next step is to discard in each such triangulation the simplices that obviously will not be positively decorated by C. With the reduced number of triangulations one can now search through all of them for the ones giving the biggest number of simultaneously positively decorated simplices. Each set of k simultaneously positively decorated simplices gives a candidate for a region of multistationarity with k positive steady states. If one finds m of such sets, then it is possible to have up to m such regions. Have in mind, however, that among these regions can be repetitions.
Next we apply this to the n-site phosphorylation system with all intermediates and explain more concretely this procedure in this case.
In Corollary 1.5 we obtained regions of multistationarity with 2[ n 2 ] + 1 positive steady states each using only 1/4 of the intermediates, our objective now is to understand if it is possible to find more such regions with more intermediates. Consider the network G of the n-site phosphorylation system with all possible intermediates: In Section 4 of [1], the concentration at steady state of all species are given in terms of the species s 0 , e, f : This looks very similar to (2.3), where F tot is replaced by f , which is now a variable, and so we need to work in dimension 3. The main difference is that the networks G J considered in the previous sections have intermediates only in the E component and the network G we consider in this section has all intermediates.
Note from (1.2) that the support A of this system, which has 2n + 4 elements, ordering the variables as s 0 , e, f, is given by the columns of the following matrix Recall that if one multiplies a column of a matrix C by a positive number, then a simplex is positively decorated by C if and only if it is positively decorated by the modified matrix. So, in order to test whether a simplex with vertices in A is positively decorated by C is enough to test if it is positively decorated by the following matrix (2) With L 1 compute L 2 by discarding all simplices which do not have the last vertex (0, 0, 0). In fact we only need these simplices since a simplex not containing the last vertex cannot be positively decorated, because the corresponding coefficients of Csimple will be all positive. (3) Compute L 3 from L 2 by removing all simplices with the corresponding 3 × 4 submatrix of Csimple having a zero 3 × 3 minor. The reason for this is clear, such simplices will never be positively decorated by Csimple. (4) Compute L 4 from L 3 using the symmetry of Csimple. More precisely, change any index 4, 5, . . . , n + 3 on the simplex to 1 because on Csimple they yield the same column. Here we are using the easy-to-check fact that changing the order of indexes does not change the conditions for a simplex to be positively decorated. Step (1) can be done with the package TOPCOM inside SAGE [22], the other steps are quite simple to implement, for instance in MAPLE [19]. We show in the table below the number of elements in some of the lists and an approximation of the computation time for small values of n. The most computationally expensive part is to compute all regular triangulations, taking at least 90% of the time. These computations were done in a Linux virtual machine with 4MB of RAM and with 4 cores of 3.2GHz of processing. With a faster computer or more time one probably can do n = 6 or even n = 7 but probably not much more than this. For n = 5 just the file for the raw list L 1 of regular triangulations already has 10Mb.
An alternative path to Steps (6) and (7) is to set a number k and look for sets T ∈ L 5 and S ⊂ T with #S ≥ k such that there is viable N 0 , . . . , N n−1 such that all ∆ ∈ S are positively decorated by Csimple at the same time. We actually used this with k = 2[ n 2 ]+1. This other route depends upon a good guess one may previously have at how many positive steady states to expect. After Step (7) one has to determine if there are any repetitions among the candidates for regions of multistationarity in L 7 and also if there are any superfluous candidates of regions, that is conditions C 1 and C 2 such that C 1 implies C 2 . In our case we did it by hand since the #L 7 was quite small. Once Step (7) is done, one has a list of inequalities for each element S of L 7 . These come from the conditions imposing that the simplices in S are positively decorated by Csimple. We are going to use these conditions to describe the regions of mulstistationarity. Because of the uniformity of Csimple the only kind of conditions that appear are or the opposite inequalities, and these translate from the N i to the k cat i , cat i as follows Note that • conditions (III) i and (V) together imply (II) i ; • the opposite of condition (II) i together with condition (V) imply the opposite of (III) i ; • the opposite of condition (III) i together with condition (V) imply (IV) i ; • the opposite of condition (IV) i together with condition (V) imply (III) i ; • condition (III) i and the opposite of (III) j together imply (I) i,j . Using these properties it is easy to describe in a nice manner the regions of multistationarity and discard the repeated and superfluous ones. We sum up our findings on the following results which are proved in the same fashion as Theorems 1.2 and 1.4, once you have the regular triangulation obtained with the computer script. In the following propositions we describe the regions of multistationarity for n = 2, 3, 4 and 5.
Proposition 5.2. Let n = 2. Assume that S tot > E tot + F tot . Then there is a choice of reaction rate constants for which the distributive sequential 2-site phosphorylation system admits 3 positive steady states. More explicitly, given rate constants and total concentrations such that after rescaling of the k on 's and on 's the distributive sequential 2-site phosphorylation system has 3 positive steady states.
Proposition 5.3. Let n = 3. Assume that S tot > E tot + F tot . Then, there is a choice of rate constants for which the distributive sequential 3-site phosphorylation system admits at least 3 positive steady states. More explicitly, if the rate constants and total concentrations are in one of the regions described below , then after rescaling of the k on 's and on 's the distributive sequential 3-site phosphorylation system has at least 3 positive steady states.
Proposition 5.4. Let n = 3. If the rate constants and total concentrations are in one of the regions described below then after rescaling of the k on 's and on 's the distributive sequential 3-site phosphorylation system has at least 3 positive steady states.
Proposition 5.5. Let n = 4. Assume that S tot > E tot + F tot . Then, there is a choice of rate constants for which the distributive sequential 4-site phosphorylation system has at least 5 steady states. More explicitly, if the rate constants and total concentrations are in one of the regions described below , then after rescaling of the k on 's and on 's the distributive sequential 4-site phosphorylation system has at least 5 steady states.
Proposition 5.7. Let n = 5. If the rate constants and total concentrations are in one of the regions described below then after rescaling of the k on 's and on 's the distributive sequential 5-site phosphorylation system has at least 5 positive steady states.
Note that the conditions in this section describe different regions from the ones described by the inequalities in Theorem 1.2 and Theorem 1.4. For instance, in Propositions 5.2, 5.3, 5.5, 5.6 the inequalities between the reaction rate constants and total conservations constants do not involve the value of S tot . In Propositions 5.4 and 5.7, the conditions onse the linear conservation constants are also different (e.g. on Ftot Etot and Stot Etot −1 instead of the product F tot ( Stot Etot − 1)). The inequalities in Theorem 1.2 and Theorem 1.4 hold for reactions rate constants of a reduced system G J , but if we use Theorem 2.4 to extrapolate these conditions to the full n-site phosphorylation network, the regions are different as well.

Discussion
We developed in this paper both the theoretical setting based on [1,2] and the algorithmic approach that follows from it, to describe multistationarity regions in the space of all parameters for subnetworks of the n-site sequential phosphorylation cycle, where there are up to 2[ n 2 ] + 1 positive steady states with fixed linear conservation constants. We chose to assume that the subnetworks we consider have intermediate species only in the E component, but of course there is a symmetry in the network interchanging E with F , each S i with S n−i , the corresponding intermediates and rate constants, and completely similar results hold if we assume that there are only intermediates in the F component. These regions can be lifted to regions of multistationarity with at least these number of stoichiometrically compatible positive steady states of the full n-site sequential phosphorylation cycle. The main feature of our polyhedral approach is that it gives a systematic method applicable to other networks (for instance, to enzymatic cascades [14]), and it provides open conditions and not only choices of parameter values where high multistationarity occurs. Moreover, we show that it can be algorithmically implemented.
For instance, Theorem 1 in [24] states that for any positive values of S tot , E tot and F tot , there exists ε 0 > 0 such that if Etot Stot < ε 0 , then there exists a choice of rate constants such that the system admits 2[ n 2 ] + 1 positive steady states. Also, in the recent paper [10], the authors prove the existence of parameters for which there are 2[ n 2 ] + 1 positive steady states and nearly half of them are stable and the other half unstable, using a different degeneration of the full network. However, the maximum expected number of stoichiometrically compatible steady states for the n-site system is equal to 2n − 1 (this is an upper bound by [24]), which has been verified for n = 3 and 4 in [13]. Probably, it is not possible to find a region in parameter space with this number of positive steady states using degeneration techniques.
In [5] and [3] it is proved that if k cat 0 / cat 0 < k cat 1 / cat 1 then there exist constants E tot , F tot , S tot such that the distributive sequential 2-site phosphorylation presents multistationarity. In [9] the author gives new regions of multistationarity for n = 2, more precisely it is shown in Theorem 1 that given parameters K 1 , L 0 , L 1 , k cat 0 , k cat 1 , cat 0 , cat 1 , for any K 0 small enough there exist constants E tot , F tot , S tot such that the distributive sequential 2-site phosphorylation presents multistationarity. How small K 0 has to be taken is explicitly given in terms of the other parameters by a rather involved equation (a similar result is obtained interchanging K 0 and L 1 ). Our Proposition 5.2 gives a similar result, but with the advantage that we only need to rescale k on 0 , k on 1 , on 0 , on 1 .
Proposition 5.2 is in agreement with [4,Corollary 4.11] which establishes that in order for the distributive sequential 2-site phosphorylation to present multistationarity the total concentration of the substrate needs to be larger than either the concentration of the kinase or the phosphatase (S tot > F tot or S tot > E tot ). In the regions of multistationarity we found for n = 3, 4 and 5 this is the case as well. Propositions 5.2, 5.3, 5.4, 5.5, 5.6 and 5.7 are of the same flavor as Theorems 1.2 and 1.4, in the sense that all of them give sufficient conditions on the rate constants and linear conservation constants such that after rescaling of the k on 's and on 's, the distributive sequential n-site phosphorylation system is multistationary. The computational approach described in Section 5 can be used for other networks of moderate size.
In conclusion, to the best of our knowledge other results (for instance [3], [5] and [9]), only give regions of multistationarity as precise as ours for the distributive sequential 2-site phosphorylation system. The present work has several advantages: it gives new regions of multistationarity for the distributive sequential n-site phosphorylation system for any n; for small n it gives several distinct regions of multistationarity; it gives lower bounds in the number of positive steady states in each of such regions; and moreover, our approach is easily applicable to other networks, both analytically and computationally.