Core signalling motif displaying multistability through multi-state enzymes

Bistability, and more generally multistability, is a key system dynamics feature enabling decision-making and memory in cells. Deciphering the molecular determinants of multistability is thus crucial for a better understanding of cellular pathways and their (re)engineering in synthetic biology. Here, we show that a key motif found predominantly in eukaryotic signalling systems, namely a futile signalling cycle, can display bistability when featuring a two-state kinase. We provide necessary and sufficient mathematical conditions on the kinetic parameters of this motif that guarantee the existence of multiple steady states. These conditions foster the intuition that bistability arises as a consequence of competition between the two states of the kinase. Extending from this result, we find that increasing the number of kinase states linearly translates into an increase in the number of steady states in the system. These findings reveal, to our knowledge, a new mechanism for the generation of bistability and multistability in cellular signalling systems. Further the futile cycle featuring a two-state kinase is among the smallest bistable signalling motifs. We show that multi-state kinases and the described competition-based motif are part of several natural signalling systems and thereby could enable them to implement complex information processing through multistability. These results indicate that multi-state kinases in signalling systems are readily exploited by natural evolution and could equally be used by synthetic approaches for the generation of multistable information processing systems at the cellular level.

1 A model for an allosteric kinase 1

.1 Model description
We consider a reaction network consisting of an allosteric kinase for one substrate. We let K be the kinase that exists in two conformations: K r (relaxed state) and K t (tensed state). Each of the conformations acts as a kinase for a common substrate S. We let S p denote the phosphorylated form of the substrate. We assume that the intermediate kinase-substrate complexes, K r S and K t S, also undergo conformational change.

Summary of results
The results for the model with one allosteric kinase can be summarised in the following way. In subsection 1.4 we show that the steady states of the system can be given in terms of the concentration x 5 of the substrate S only. That is, for a given K tot knowing the value of x 5 at steady state allows us to calculate the value of the remaining concentrations at steady state. Further, we show that the system can have up to three positive steady states by choosing the reaction rate constants and total amounts appropriately.
If the inequality is not fulfilled, then there cannot be multistationarity for any choice of total amounts. Moreover, by inspecting the inequality, other necessary conditions for multistationarity might be deduced. For example, one of the following two constraints is necessary for multistationarity to occur: (a) κ 3 > κ 6 and η r κ 9 κ 10 > η t κ 8 κ 11 .
If a set of reaction rate constants fulfils the necessary and sufficient conditions for multistationarity, then the next question is to find total amounts for which it occurs. The linear inequalities in S tot and K tot given in (2) restrict the possible values considerably. However, it is also possible to give necessary and sufficient conditions involving all parameters, that is, the reaction rate constants and the total amounts. These conditions are given in equations (15) and (16) in subsection 1.6.2. In subsection 1.7 we discuss how to explicitly find parameter sets for which multistationarity arises, using the conditions discussed above, and illustrate it with an example.
In subsection 1.8 we explore the relationship between the species concentrations at steady state and the concentration of the substrate S (x 5 ) at steady state. To this end we use the explicit descriptions of the species concentrations at steady state as functions of x 5 , found in subsection 1.4. Furthermore we illustrate that the steady states of the system cannot be given in terms of the concentration x 6 of the modified substrate S p alone. The reason is that x 5 cannot be expressed in terms of x 6 when multistationarity occurs. This occurs because the substrates S and S p do not appear in a symmetric way in the reactions.
In subsection 1.9, we consider bifurcation plots in the multistationary setting. We study the effect of changing the total amounts of the substrate and the kinase on the number of steady states. We encounter here again the necessary and sufficient conditions from subsection 1.6.

Model reduction and multistationarity
In order to study the number of steady states of the networks we are interested in, we make use of a couple of results that allow us to simplify the network while keeping information on the number of steady states.
The first result concerns removal of intermediates. Intermediates are species in the network that do not appear interacting with any other species. For example the species K r S and K t S are intermediates.
Given a network, we can obtain a reduced network by "removing" some intermediates, one at a time. This is done in the following way. Say we want to remove an intermediate Y from the network. We remove all reactions of the original network that involve Y and add a reaction y → y whenever y → Y → y with y = y belongs to the original network. Here y and y are the reactant of a reaction y → Y and product of a reaction Y → y , respectively.
To illustrate this, we consider the removal of the intermediate K r S of the network in subsection 1.1. The reactions of the reduced network are obtained by considering all length 2 paths of the original network that go through K r S. We have four such paths: Clearly the process can be repeated now choosing another intermediate, for example K t S. In this way we can obtain reduced networks by removing several intermediates.
Our interest in reduction by removal of intermediates is explained in the following theorem: if the reduced network admits N positive steady states, then so does the original network. This theorem and its proof can be found in [4].
. Consider a chemical reaction network G and let G be the reduced network obtained by the removal of some intermediates as described above. If G has N non-degenerate positive steady states for some choice of reaction rate constants and total amounts, then there is a choice of reaction rate constants and total amounts such that G has at least N corresponding non-degenerate positive steady states. Moreover, if all the steady states are hyperbolic, then the correspondence preserves the asymptotic stability of the steady states.
The condition that the steady state is non-degenerate and hyperbolic is technical. A steady state is non-degenerate if the Jacobian map of the system when evaluated at the steady state and restricted to the stoichiometric subspace (defined by the conservation law equations equated to zero), is non-singular. A hyperbolic steady state is a steady state for which the eigenvalues of the Jacobian of the system when evaluated at the steady state and restricted to the stoichiometric subspace have non-zero real part. These conditions are fulfilled in our applications of the theorem.
The other result we are interested in concerns subnetworks. A subnetwork of a network is obtained by removing some of the reactions of the original network. For example, by making some reversible reactions irreversible, we obtain a subnetwork. Observe that removal of a reaction corresponds to setting its reaction rate constant to zero in the ODE system. For example, making the reversible reaction K r + S Under some conditions, knowing that a subnetwork admits N positive steady states is enough to conclude that this is also the case for the original network. A sufficient condition for this to hold is that the subnetwork and the network have the same species and conservation laws (or stated more precisely, the two networks have the same stoichiometric subspace). This theorem and its proof can be found in [6].
Theorem 2 ([6]). Let G be a chemical reaction network and G be a subnetwork of G such that G and G have the same stoichiometric subspace. If there exist reaction rate constants and total amounts such that G has N non-degenerate positive steady states, then there exist reaction rate constants and total amounts such that G has at least N non-degenerate positive steady states. Moreover, if all the steady states are hyperbolic, then the correspondence preserves the asymptotic stability of the steady states.
The theorem can be applied to subnetworks obtained by making some reversible reactions irreversible, because the conservation laws, and thus the stoichiometric subspace, are not changed.
We will use these two reduction strategies in sections 2 and 3. In both sections we find the number of steady states of a reduced network, obtained either by removal of intermediates, reactions or both. These two results allow us then to deduce that the original network of interest has at least the same number of steady states (for some choice of reaction rate constants and total amounts).
Note that reduced networks might not be biologically realistic. In this text we use reduction as a mathematical simplification tool to obtain results for the networks of interest.
The expressions for x 1 , x 2 , x 3 , x 4 , x 6 are positive provided x 5 is positive.
We note that this polynomial has at least one positive root since p(0) < 0 and p(+∞) > 0. Moreover, the reaction rate constants and the total amounts can be chosen such that p(x 5 ) has indeed three positive roots. Therefore, there exist reaction rate constants and total amounts such that the system has three positive steady states. This is first shown by giving a set of parameters such that the system has three positive steady states in subsection 1.8.
A formal mathematical proof of this result is given in subsection 2.4 where we study a generalisation of this model. The proof is based on showing that the result holds for a reduced model and then using Theorem 2 (Section 1.3).

Necessary conditions for multistationarity
In this section we use the Descartes' rule of signs on the number of positive roots of a polynomial. The rule says that the number of positive roots of a polynomial is smaller than or equal to the number of sign changes in the sequence of signs of the coefficients of the polynomial (ignoring the coefficients that are zero).
Proof. Following Descartes' rule of signs, a necessary condition for p(x 5 ) to have three positive roots is that the coefficients of the polynomial have alternating signs, that is, there must be three sign changes in the sequence of signs of the coefficients of the polynomial and none of these coefficients can be zero. Since the leading coefficient of p(x 5 ) is positive and the independent term is negative, a necessary condition is that the coefficient of degree 2 is negative (this gives the first condition in the statement) and the coefficient of degree 1 is positive (this gives the second condition in the statement).
An example of how such a sector might look like is given in Example 7 below.

Conditions involving only reaction rate constants
In order to find sufficient conditions for multistationarity involving only reaction rate constants, we apply the strategy introduced in [1]. The strategy is based on Brouwer Degree Theory. In [2], sufficient conditions for multistationarity, based on the reaction rate constants only, were found for a two-site phosphorylation cycle in which both the kinase and the phosphatase follow a sequential and distributive mechanism.
Proposition 4. For a given choice of reaction rate constants κ, there exist total amounts K tot and S tot such that the network admits multiple positive steady states if and only if where η r = κ 1 κ 2 + κ 3 and η t = κ 4 κ 5 + κ 6 are the inverses of the Michaelis-Menten constants for the different allosteric states of the kinase, i.e. for K r and K t respectively.
Equivalently, the network admits multiple positive steady states if and only if is positive.
Proof. We apply the strategy from [1]. The steps of the procedure are as follows: (i) Consider the function f κ (x) given by the two conservation laws and the expressions foṙ x 2 ,ẋ 3 ,ẋ 4 andẋ 6 , that is, the right-hand side of the expressions in (3) together with the second equation in (1): Compute the Jacobian matrix of f κ (x): and its determinant, det(J κ (x)).
(ii) Express x 2 , x 3 , x 4 and x 6 at steady state in terms of x 1 and x 5 . That is, consider the steady state equationsẋ 2 =ẋ 3 =ẋ 4 =ẋ 6 = 0 and solve them for x 2 , x 3 , x 4 , x 6 in terms of x 1 , x 5 . This gives the following expression: (iii) Substitute the values of x 2 , x 3 , x 4 , x 6 from (14) into the determinant of the Jacobian matrix in (13). The resulting expression is the quotient of a polynomial b κ (x 1 , x 5 ) in x 1 , x 5 and the polynomial χ κ (x 5 ). Note that all coefficients of the denominator are positive.
Based on Brouwer Degree Theory (see [1]), we conclude that multistationarity occurs if and only if the polynomial b κ (x 1 , x 5 ) is positive for some positive values of x 1 , x 5 . Specifically, for a given choice of reaction rate constants κ, there exist total amounts such that the system displays multistationarity if and only if b κ (x 1 , x 5 ) is positive for some positive values of x 1 , x 5 .
Viewed as a polynomial in Further all coefficients of b κ (x 1 , x 5 ) are polynomials in κ. All coefficients have negative sign, independently of the values of κ i , except for the coefficient of x 1 x 2 5 which is α(κ). Clearly, if this coefficient is also negative, then all coefficients of the polynomial are negative and thus b κ (x 1 , x 5 ) is negative for all positive values of x 1 , x 5 . As a consequence multistationarity cannot occur for this choice of reaction rate constants κ.
Assume now that α(κ) is positive. We show that in this case the polynomial b κ (x 1 , x 5 ) is positive for some values of x 1 , x 5 . Let x 5 = T and x 1 = T 2 . Then b κ (T 2 , T ) is a polynomial in T of degree 4 with leading positive coefficient α(κ). By letting T be arbitrarily large, b κ (T 2 , T ) becomes eventually positive. We conclude that b κ (x 1 , x 5 ) is positive for some values of x 1 , x 5 and hence that there exist total amounts such that the network with reaction rate constants κ has more than one positive steady state.
This shows that multistationarity occurs if and only if α(κ) is positive. After rearranging the terms of the coefficient α(κ), the inequality α(κ) > 0 is equivalent to the inequality given in the statement of the proposition, i.e. inequality (11).
By inspecting inequality (11), we can find some simple necessary conditions for multistationarity. For example: • Either κ 9 κ 10 = 0 or κ 8 κ 11 = 0 (or both are non-zero). That is, allosteric changes must occur both for the kinase and the kinase-substrate complexes.
• Since the left-hand side of the inequality must be positive for the inequality to hold, we deduce that one of the following conditions is necessary for multistationarity: (a) κ 3 > κ 6 and η r κ 9 κ 10 > η t κ 8 κ 11 .

Conditions involving reaction rate constants and total amounts
Here we provide necessary and sufficient conditions on all parameters (reaction rate constants and total amounts) of the system for multistationarity to occur. To obtain the conditions, we apply Sturm's Theorem: Theorem 5 (Sturm). Let p(x) be a real polynomial. Define recursively the Sturm sequence by Proposition 6. Given reaction rate constants κ and total amounts K tot and S tot , the network has three positive steady states if and only if the following inequalities are satisfied Proof. We are interested in the positive roots of the polynomial p(x) = p(x 5 ) in (9). In order to apply Sturm's theorem, we should take a = 0 and b so large that all positive roots are in (a, b] = (0, b]. If b is chosen arbitrarily large, then the signs of p 0 (b), . . . , p m (b) are determined by the leading coefficients of the polynomials p 0 , . . . , p m . Because b is an arbitrarily large number, we write b = +∞ and the sequence is written as p 0 (+∞), . . . , p m (+∞). Observe that a = 0 is not a root of p(x) and hence the hypothesis of Sturm's theorem applies.
Therefore, for σ(+∞) = 0 we require p 2 (+∞), p 3 (+∞) > 0. The polynomial p 3 (x) has degree zero, and hence p 3 (0) = p 3 (+∞). Therefore, we are left with 4 conditions on the parameters that fully characterise the region of the parameter space with three positive steady states, namely Using that a 0 > 0 and a 3 < 0, these conditions simplify to the following conditions:

Necessary and sufficient conditions in practice
In order to find explicit values of the parameters such that the system exhibits multistationarity, the procedure is the following: 1. First, use the necessary and sufficient condition given in (11) to find appropriate values for the reaction rate constants.
2. Second, substitute these values of the reaction rate constants into (15). This yields a system of 4 inequalities in K tot and S tot . The positive values of (K tot , S tot ) fulfilling the inequalities correspond to parameter sets for which there are three positive steady states. By the results above, there are always values of (K tot , S tot ) for which this is the case.
After the first step we might use the necessary conditions for multistationarity from subsection 1.5, that is, using the inequalities in (10) instead of the conditions in (15). This gives (simpler) regions of the parameter space of total amounts containing all pairs (K tot , S tot ) for which there is multistationarity. Remember though that not all pairs (K tot , S tot ) satisfying the inequalities in (10) yield multistationarity as the conditions are only necessary and not sufficient.
Example 7. Consider the set of parameters for which (11) is satisfied. The system of inequalities (10) is 110.1K tot + 3.06 < S tot < 144.16K tot + 0.65, and the pairs (K tot , S tot ) fulfilling the inequalities are highlighted in light blue in Figure 1(a,b).  The plot in Figure 1b illustrates that for very small values of K tot there is no value of S tot for which the inequalities are satisfied. The dot in Figure 1a corresponds to S tot = 591 and K tot = 5, for which there is multistationarity.
One set of total amounts fulfilling the above system of inequalities is S tot = 591 and K tot = 5 (the point plotted in the figure above). In the following two plots (Figure 2(a,b)) the yellow region is the region of common solutions to the first, second and fourth inequalities, and the blue region the solution to the third inequality. The intersection of the two regions is the solution set to the four inequalities. It is the small blue region inside the yellow region. Figure 2b    The region for which multistationarity exists is much smaller than the region given in Figure 1(a,b). Note that this example illustrates the region for S tot and K tot leading to multistationarity only for the specific set of parameters given in (17).
In practice, it is not straightforward to solve the Sturm's inequalities for S tot and K tot .

Describing the steady states
In this section we describe the intersection of the solution set to the steady state equations, with the linear space defined by x 1 + x 2 + x 3 + x 4 = K tot , using the expressions (4)- (8). To illustrate the results of this section we choose the set of parameters for which multistationarity occurs given in (17): κ 7 = 0.01, κ 8 = 0.8, κ 9 = 0.1, κ 10 = 0.01, κ 11 = 0.1.

The concentration [S p ] as a function of [S]
First, we discuss how the concentration of S p , x 6 , changes according to the concentration of S, x 5 , at steady state using (8).
We are interested in the steady states for a fixed value of S tot = x 3 + x 4 + x 5 + x 6 . Therefore, we also consider the rational function of x 6 obtained by substitution of (6) (x 3 = ϕ 3 (x 5 )) and (7) (x 4 = ϕ 4 (x 5 )) into the conservation law for S tot : This expresses x 6 by another rational function in x 5 whose numerator has degree three and denominator has degree two. We plot the two functions using the parameters in (18), K tot = 5 and S tot = 591. Figure 3b. Figure 3a shows the graph of the function ϕ 6 (x 5 ). Figure 3b shows the graph of the function ϕ 6 (x 5 ) (in blue) together with the function S tot − ϕ 3 (x 5 ) − ϕ 4 (x 5 ) − x 5 (in red). The intersection points of the two graphs in Figure 3b (marked with dots) are the pairs (x 5 , x 6 ) for the three steady states corresponding to the chosen values of the total amounts K tot and S tot .
We observe in Figure 3a that ϕ 6 (x 5 ) increases for small values of x 5 , until it reaches a maximum and then decreases towards a limit value for large x 5 . In the next proposition we establish that this shape of ϕ 6 (x 5 ) is actually necessary for multistationarity.
Recall that by Proposition 4, for a given choice of reaction rate constants κ, there exist total amounts K tot and S tot such that the network admits multiple steady states if and only if α(κ) > 0. We first list some properties of the function ϕ 6 (x 5 ).
Proposition 8. Let the reaction rate constants κ and the total amount K tot be given.
• If α(κ) ≤ 0 (i.e. the network is not multistationary), then ϕ 6 (x 5 ) is either an increasing function for all x 5 or there exists x 5 such that ϕ 6 (x 5 ) is increasing for x 5 < x 5 and is decreasing for Proof. It is straightforward to see that ϕ 6 (x 5 ) > 0 unless x 5 = 0 and to compute the limit value using (8). The derivative of ϕ 6 (x 5 ) with respect to x 5 has the form where λ 2 (κ), λ 3 (κ) are positive polynomials in the reaction rate constants, and with α(κ) as in (12) and λ 1 (κ) a positive polynomial in the reaction rate constants. If λ 1 (κ) is positive, then ϕ 6 (x 5 ) > 0 for all x 5 > 0 and hence ϕ 6 (x 5 ) is an increasing function.
If λ 1 (κ) is negative, then the numerator of ϕ 6 (x 5 ) is a second degree polynomial with negative leading coefficient K tot λ 1 (κ), and the rest of the coefficients are positive. This implies that the polynomial has a unique positive root, that is, there exists a unique x 5 > 0 for which the polynomial vanishes (and hence ϕ 6 ( x 5 ) = 0). The numerator is negative for x 5 > x 5 and positive for x 5 < x 5 . As a consequence, the derivative of ϕ 6 (x 5 ) is positive for x 5 < x 5 and negative for x 5 > x 5 . It follows that ϕ 6 (x 5 ) is increasing for x 5 smaller than x 5 and is decreasing for x 5 larger than x 5 .
The proposition states that x 6 = ϕ 6 (x 5 ) settles at a non-zero value when x 5 , and hence also S tot , are large (second statement in the proposition). When multistationarity occurs, Figure 3a illustrates that there are values of x 6 that correspond to two different values of x 5 (third statement in the proposition). It follows that we cannot express x 5 as a function of x 6 , and hence we cannot express the other concentrations at steady state in terms of x 6 alone.

The concentrations [K r ] and [K t ] as functions of [S]
We consider the rational functions x 1 = ϕ 1 (x 5 ) in (4) and x 2 = ϕ 2 (x 5 ) in (5). The plot of these functions using the parameters in (18) and K tot = 5 are shown in Figure 4(a,b). Proposition 9. Let the reaction rate constants κ and the total amount K tot be given.
For small values of x 5 , the sign of the derivative of ϕ 1 (x 5 ) and ϕ 2 (x 5 ) agrees with the sign of a 3 (κ) and b 3 (κ), respectively. If a 3 (κ) is negative, then ϕ 1 (x 5 ) decreases for all x 5 . If a 3 (κ) is positive, then ϕ 1 (x 5 ) vanishes at exactly one value x 5 . In this case ϕ 1 (x 5 ) increases for x 5 < x 5 and decreases for x 5 > x 5 . We argue analogously for the monotone behaviour of ϕ 2 (x 5 ) depending on the sign of b 3 (κ).

The concentrations [K r S] and [K t S] as functions of [S]
We consider the rational functions x 3 = ϕ 3 (x 5 ) in (6) and x 4 = ϕ 4 (x 5 ) in (7). The plot of these functions using the parameters in (18) and K tot = 5 are shown in Figure 5(a,b). Proposition 10. Let the reaction rate constants κ and the total amount K tot be given.
• We have lim • For both ϕ 3 (x 5 ) and ϕ 4 (x 5 ) it holds that: the function either increases for all x 5 , or there exists a value x 5 such that ϕ 3 (x 5 ) increases for x 5 < x 5 and decreases for x 5 > x 5 .
Proof. From (6) and (7), it follows that ϕ 3 (x 5 ) = 0 and ϕ 4 (x 5 ) = 0 only if x 5 = 0. Further, it follows that the numerator and the denominator of both ϕ 3 (x 5 ) and ϕ 4 (x 5 ) have degree two. By inspecting the leading coefficient of the numerator and the denominator of these functions, we find the limit values given in the statement. The derivatives of ϕ 3 (x 5 ) and ϕ 4 (x 5 ) are where c 2 (κ), c 3 (κ), d 2 (κ), d 3 (κ) are positive polynomials in the reaction rate constants, and c 1 (κ) = c 1 (κ)(κ 1 κ 10 (κ 11 − κ 9 ) (κ 5 + κ 6 ) + κ 4 κ 11 (κ 8 + κ 11 ) (κ 2 + κ 3 )), with c 1 (κ) and d 1 (κ) positive polynomials in the reaction rate constants. The only coefficients in the numerators of ϕ 3 (x 5 ) and ϕ 4 (x 5 ) that have undetermined sign are thus c 1 (κ) and d 1 (κ) and the derivatives can change sign at most once. It follows that for small values of x 5 , both derivatives take positive values and the functions are increasing. For large x 5 , the derivatives are positive or negative, depending on the signs of c 1 (κ) and d 1 (κ) (and thus on the reaction rate constants). Therefore, each of the functions ϕ 3 (x 5 ) and ϕ 4 (x 5 ) can either be increasing towards the limit (19) (as in Figure 5a) or be increasing towards a maximum value and then be decreasing towards the limit in (19) (as in Figure 5b).
These results show that for a steady state with a large concentration of S, the concentrations of K r S and K t S are close to a limit value, and the concentrations of K r and K t are close to zero. This confirms mathematically that saturation occurs: large amounts of substrate imply that the kinase is essentially only in bound form.

Bifurcation plots with K tot and S tot
In this subsection we investigate how the number of positive steady states depends on the total amounts of kinase and substrate.

Changing S tot
We fix first the total amount of kinase K tot and analyse ψ S (x 5 ). We start with an example. We let K tot = 5 and choose the reaction rate constants as in (18). We plot the value of S tot against x 5 at steady state. Figure 6b is the graph of ψ S (x 5 ), while Figure 6a is the bifurcation plot obtained by interchanging the axes of the plot in Figure 6b. Figure 6a. Figure 6b.
From Figure 6a, we conclude that multistationarity occurs for values of S tot in a certain interval (highlighted in blue in the figure): the values of S tot that correspond to three values of x 5 (which in turn give rise to three positive steady states). For low values of S tot the concentration of x 5 is low and for higher values of S tot the concentration of x 5 is high.
To understand Figure 6a, we study the function ψ S (x 5 ) in Figure 6b, since the bifurcation plot is simply obtained by interchanging the axes. We do this because we do not have an analytical expression of the type x 5 = Φ(S tot ).
Proposition 11. Fix the reaction rate constants κ and the total amount K tot .
• If there exists S tot such that the network has multiple positive steady states, then ψ S (x 5 ) has exactly one local maximum, one local minimum and the monotonicity pattern of ψ S (x 5 ) is increasing-decreasing-increasing (as shown in Figure 6b).
• If the network does not admit multiple positive steady states for any S tot , then ψ S (x 5 ) is an increasing function.
• If α(κ) > 0, then the network admits multiple positive steady states for some S tot , provided K tot is large enough.
Proof. The function ψ S (x 5 ) is a rational function whose numerator has degree three and the denominator has degree two. The coefficients of the numerator and denominator of ψ S (x 5 ) are positive polynomials in the reaction rate constants and K tot . Hence, ψ S (x 5 ) tends to infinity as x 5 goes to infinity. The derivative of ψ S (x 5 ) with respect to x 5 is the rational function where µ 1 (κ) and µ 2 (κ) are positive polynomials in the reaction rate constants, µ 4 (κ, K tot ) and µ 5 (κ, K tot ) are positive polynomials in the reaction rate constants and K tot (depending linearly on K tot ), q 2 (κ, x 5 ) is a degree 4 polynomial in x 5 whose coefficients are positive polynomials in the reaction rate constants and with µ 3 (κ) a positive polynomial in the reaction rate constants and α(κ) as in (12). Thus the only coefficient of the numerator of ψ S (x 5 ) that can be negative is µ 3 (κ, K tot ).
Since the independent and leading coefficients of the numerator of ψ S (x 5 ) are positive, ψ S (x 5 ) is positive for small and large enough values of x 5 and the function ψ S (x 5 ) is increasing.
Multistationarity occurs if and only if there are values of x 5 for which the corresponding values of S tot = ψ S (x 5 ) agree. This occurs if and only if ψ S (x 5 ) decreases in some interval (it cannot be an increasing function) and hence necessarily µ 3 (κ, K tot ) is negative. In this case the sequence of coefficients of the polynomial in the numerator of ψ S (x 5 ) has two changes of sign and the Descartes' rule of signs tells us that ψ S (x 5 ) = 0 has at most two solutions. It follows that the function ψ S (x 5 ) cannot change the increasing/decreasing behaviour more than twice. Combined with the discussion above, multistationarity implies that the monotonicity pattern of ψ S (x 5 ) is increasing-decreasing-increasing. We deduce that there is exactly one local maximum and one local minimum when the system displays multistationarity and that the network has three steady states (and not more). We conclude that the graph of ψ S (x 5 ) must be as in Figure 6b when multistationarity occurs, that is, it has an S-shape.
In Figure 7a, the values of S tot for which multistationarity occurs can be identified (highlighted in blue in the figure). In Figure 7b, we used a larger range of x 5 values in the plot. We find that x 6 increases with S tot towards a maximum value after which x 6 decreases. After the maximum value there is a transition phase where multistationarity occurs, and finally for larger S tot there is one steady state. In this case it is not straightforward to pursue a general analysis of the relationship between S tot and x 6 , as we did for S tot and x 5 , because we cannot express x 5 , and hence S tot , as a function of x 6 .

Changing K tot
We consider the value of S tot fixed and analyse K tot = ψ K (x 5 ). Figure 8 shows the bifurcation plot obtained by plotting the function ψ K (x 5 ) with S tot = 591 and the reaction rate constants as in (18), with the axes interchanged. The region with three positive steady states is highlighted in blue.

Figure 8.
To investigate the behaviour of ψ K (x 5 ), we plot the function ψ K (x 5 ) for two different domains of x 5 , see Figure 9(a,b). Figure 9a. Figure 9b.
Proposition 12. Fix the reaction rate constants κ and the total amount S tot .
• ψ K (x 5 ) is well defined and continuous for x 5 > 0, tends to plus infinity as x 5 tends to zero, is positive only when 0 < x 5 < S tot and it holds that ψ K (S tot ) = 0.
• If there exists K tot such that the network has multiple positive steady states, then ψ K (x 5 ) has exactly one local maximum, one local minimum and the monotonicity pattern of ψ K (x 5 ) is decreasing-increasing-decreasing (as shown in Figure 9a).
• If the network does not admit multiple positive steady states for any K tot , then ψ K (x 5 ) is a decreasing function.
• If α(κ) > 0, then the network admits multiple positive steady states for some K tot , provided S tot is large enough.
Proof. The function ψ K (x 5 ) is a rational function whose numerator has degree three and the denominator has degree two (see (21)). The denominator of ψ K (x 5 ) vanishes when x 5 = 0 and hence ψ K (x 5 ) is only defined and continuous for x 5 > 0. Since the independent term of the numerator is positive, we deduce that ψ K (x 5 ) (and thus K tot ) tends to plus infinity as x 5 tends to zero. We next analyse when ψ K (x 5 ) is positive. The denominator of ψ K (S tot ) is positive for all strictly positive values of x 5 . It is easy to see that ψ K (S tot ) = 0. The sequence of coefficients of the polynomial in the numerator of ψ K (x 5 ) has exactly one change of sign. Thus by Descartes' rule of signs, x 5 = S tot is the only positive root of the numerator. Since the independent term of the numerator is positive, the function ψ K (x 5 ) is positive for x 5 < S tot and negative for x 5 > S tot . The derivative of ψ K (x 5 ) with respect to x 5 is a rational function where γ 1 (κ), γ 2 (κ) are negative polynomials in the reaction rate constants, γ 4 (κ, S tot ), γ 5 (κ, S tot ) are negative polynomials in the reaction rate constants and S tot (they are linear in S tot ), q 3 (κ, x 5 ) is a degree 4 polynomial in x 5 whose coefficients are positive polynomials in the reaction rate constants, and where γ 3 (κ) is a polynomial in the reaction rate constants with positive and negative terms, and α(κ) is as in (12). The only coefficient of the numerator of ψ K (x 5 ) that can have both signs is γ 3 (κ, S tot ). Since the leading and independent terms γ 1 (κ) and γ 5 (κ, S tot ) are negative, we deduce that ψ K (x 5 ) is negative for small and large enough values of x 5 . Thus ψ K (x 5 ) decreases in this case.
If multistationarity occurs for some values of K tot , ψ K (x 5 ) must increase in some interval, where necessarily γ 3 (κ, S tot ) > 0. We argue as above for ψ S (x 5 ) to conclude that the function has exactly one local maximum and one local minimum. In this case the graph of ψ K (x 5 ) has the same S-shape as in the example graph in Figure 9a. If multistationarity does not occur, then ψ K (x 5 ) must be decreasing for all 0 < x 5 ≤ S tot .
The first statement of the proposition tells us that the function ψ K (x 5 ) needs to be considered for 0 < x 5 ≤ S tot . K tot tends to infinity as x 5 tends to zero and is zero when x 5 = S tot . This implies that the concentration of S at steady state, x 5 , equals S tot when K tot = 0 and tends to zero as K tot tends to infinity.
We plot also the value of x 6 against K tot . To do that, we use the expression of x 6 given in (8), where we substitute K tot with ψ K (x 5 ). Then we plot the pair of functions (K tot , x 6 ) for varying values of x 5 : Figure 10a. Figure 10b.
For low values of K tot , the concentration of x 6 increases, until it reaches a transitional phase, after which it tends towards an upper bound, the total amount of substrate S tot . This bound is found by using that x 5 = 0 when K tot tends to infinity. Specifically, it is found by letting x 5 = 0 in the expression of x 6 given by (8) after replacing K tot by ψ K (x 5 ).

Positive feedback loops
The following figure shows the DSR-graph (directed-species-reaction graph) for the network in Subsection 1.1. The DSR-graph is constructed as follows: The node set is {K r , K t , K r S, K t S, S, S p , R 1 , . . . , R 11 }.
There is an edge from a species X to a reaction R with label + if X is a reactant of R. There is an edge from a reaction R to a species X if the net production of X in R is different from zero. The label is + if the net production is positive and − otherwise.
We highlight in green and red the four positive feedback loops that contribute to multistationarity, as found following the method in [5]. In this section we consider a generalisation of the model studied in the previous section. Namely, we consider the network that consists of n allosteric kinases that act on the same substrate. The network studied in the previous section corresponds to the case n = 1. We show that the network with n allosteric kinases admits 2n + 1 steady states for some choice of reaction rate constants and total amounts. In order to prove this, we use the reduction techniques explained in subsection 1.3. Explicitly, we make some reversible reactions irreversible and make use of Theorem 2.
For n = 1 the reduction corresponds to setting κ 2 = κ 5 = 0 and κ 9 = κ 10 = 0 in the model in subsection 1.1. Biologically, we can interpret the first choice as assuming strong complex formation. The second choice is compatible with condition (11). Since we aim to study the effect of multiple kinases on multistationarity, we consider conditions that are compatible with multistationarity in the single kinase case. This reduced model is still multistationary as we will show below.

Model description
We study a reduction of the system consisting of n allosteric kinases competing for the same substrate.
Let K i , for i = 1, . . . , n, denote the n allosteric kinases. We use subindices r, t to denote the relaxed or tensed state (respectively) of each of them. The set of reactions given in the previous subsection are reproduced for the n allosteric kinases, after making some reversible reactions irreversible and renaming the reaction rate constants accordingly. That is, for i = 1, . . . , n, the reactions under consideration are as follows: In addition, there is a dephosphorylation reaction In the full model, the reactions with reaction rate constants κ i,1 , κ i,3 , κ i,5 and κ i,6 are reversible.
We denote the concentration of the species as follows: for i = 1, . . . , n. We proceed as in the previous section and model the dynamics of the concentrations over time under the law of mass action by the following system of ordinary differential equations:ẋ for i = 1, . . . , n. The system has n + 1 conservation laws: one for the substrate and one for each kinase. Namely, for i = 1, . . . , n, we have for some K i tot > 0, and for S tot > 0,

Summary of results
The results for the reduced model with n allosteric kinases competing for the same substrate can be summarised in the following way. In subsection 2.3 we show that the steady states of the system can be given in terms of the concentration x 5 of the substrate S only. That is, the concentrations of the other species at steady state can be found from the value of x 5 at steady state. In subsection 2.4, we show that there exist reaction rate constants and total amounts such that the system has exactly 2m + 1 positive steady states, for all m = 0, . . . , n. In particular this is true for m = n in which case there are 2n + 1 positive steady states. We show that 2n + 1 is the maximal possible number of non-negative steady states (each concentration is positive or zero). In subsection 2.5 we consider the stability of the steady states and show that if there are 2n + 1 positive steady states, then at least n of them are unstable.
These results extend to the full model (with the reactions with reaction rate constants κ i,1 , κ i,3 , κ i,5 and κ i,6 being reversible), due to Theorem 2.

Positive steady states
The positive steady states of the system are the solutions to the equationsẋ i,1 =ẋ i,2 =ẋ i,3 = x i,4 = 0, for i = 1, . . . , n, together withẋ 5 =ẋ 6 = 0, constrained by the conservation laws (22) and (23). We reason as in the previous section and disregard the steady state equationṡ x 5 = 0 andẋ i,1 = 0, for i = 1, . . . , n. Using the equationsẋ i,2 =ẋ i,3 =ẋ i,4 = 0 and (22), we obtain algebraic expressions for x i,1 , x i,2 , x i,3 , x i,4 at steady state, depending on the value of x 5 at steady state, analogous to the expressions (4)-(7): where q i (x) is the following polynomial in x 5 : These expressions are positive provided x 5 is positive. From the equationẋ 6 = 0 we obtain which, using expressions (26) and (27), is positive provided x 5 > 0. All concentrations are expressed as functions of x 5 . We replace x 6 and subsequently x i,3 , x i,4 in (23) by their expressions in (26)-(28) to obtain After multiplying this equation by n i=1 q i (x), the left-hand side becomes a polynomial p(x 5 ) of degree 2n + 1 in x 5 . As argued in the previous section, any positive root of the polynomial corresponds to a positive steady state. We note again that p(x 5 ) has at least one positive root since one can check that p(0) < 0 and p(+∞) > 0.

Existence of 2n + 1 positive steady states.
We have shown that the positive steady states of the reduced system with n allosteric kinases competing for the same substrate are determined by the positive solutions to a polynomial p(x 5 ) of degree 2n + 1. By the fundamental theorem of algebra, a polynomial of degree 2n + 1 has 2n + 1 complex roots counted with multiplicity. Therefore, such a polynomial can at most have 2n + 1 distinct positive real roots.
We show in this section that there exist choices of reaction rate constants κ i and total amounts K i tot , S tot such that the polynomial has exactly 2n + 1 distinct positive real roots. As a consequence, this proves that the reduced system with n allosteric kinases competing for the same substrate admits 2n + 1 positive steady states for some choice of reaction rate constants and total amounts and there cannot be more. As argued at the beginning of the section, Theorem 2 guarantees that this result holds for the general system, where the conformational changes of the kinases and the binding of the kinases and the substrate are considered reversible.
The proof of this statement consists of a series of simplifications and constructions analogous to those in [7].
We plug these expressions into the equation for α 2 and obtain the equation which is equivalent to the polynomial equation The right-hand side of the equation can be seen as a polynomial of degree 2 in κ 6 with negative leader term. If the independent term is positive, then the polynomial has one positive root, which determines the value of κ 6 . Hence, we want to show that there exist values of κ 2 , κ 4 such that −κ 2 α 5 (κ 2 α 2 α 3 − κ 4 α 1 α 4 − α 1 α 4 + α 2 α 3 ) > 0 or equivalently, that For a fixed value of κ 2 , this expression is a decreasing linear function on κ 4 . Therefore, we can find a positive value of κ 4 such that it is negative. We conclude that for any κ 2 > 0, we can find values of κ 1 , κ 3 , κ 4 , κ 5 , κ 6 > 0 and K tot > 0 satisfying (35)-(39).
As a consequence of Lemma 1, there exist values of the reaction rate constants and total amounts such that the network has 2n+1 positive steady states if we can find α i,1 , . . . , α i,5 > 0 such that has 2n + 1 positive solutions.
With this notation, we want to determine the number of positive roots of the polynomial obtained by clearing denominators in (40): The coefficient of degree 2n + 1 of p(x) is n i=1 α i,3 and the independent term of p(x) is −S tot n i=1 α i,5 . We set α i,4 = 0 and α i,1 = 0. Setting these two constants to zero, for i = 1, . . . , n, does not change the degree of the polynomial. By the continuity of the isolated roots of a polynomial as functions of the coefficients of the polynomial, if we can find α i,2 , α i,3 , α i,5 > 0 such that with α i,4 = α i,1 = 0, the polynomial p(x) has 2n + 1 distinct positive real roots, then for α i,4 , α i,1 small enough, the polynomial p(x) still has 2n + 1 distinct positive real roots.
We further let α i,3 = 1 for all i = 1, . . . , n, and S tot = 1. To ease the notation, we write a i = α i,2 and b i = α i,5 , such that the polynomial of interest becomes Lemma 2. There exist a i , b i > 0, for i = 1, . . . , n, such that p(x) has 2n + 1 positive roots.
Proof. The statement is a consequence of Lemmas 2, 3 and 4 in the Supplementary Information of [7].
We are ready to prove the main result on the number of positive steady states.
Theorem 13. For any n ≥ 1, there exists a choice of reaction rate constants and total amounts such that the reduced network with n allosteric kinases competing for the same substrate has 2n + 1 distinct positive steady states.
Remark 14 (Existence of 2m + 1 steady states). Consider the polynomial p(x) in (42) and let a k = 0 for a certain k. We obtain the polynomial Assuming b k > 0, the polynomialp(x) factorises into two factors: one factor of degree two with non-real roots and the other factor is of the form of (42) with degree 2(n−1)+1 = 2n−1. By Lemma 2, the latter factor admits 2n − 1 positive roots for some choice of a i , b i . Therefore, we conclude thatp(x) admits precisely 2n − 1 positive roots for some choice of parameters as well. By the continuity of the roots of a polynomial, this remains true for a k > 0 small enough. Proceeding as in the proof of Theorem 13, this implies that there exist reaction rate constants and total amounts such that the system has precisely 2n − 1 positive steady states.
We can repeat the argument by letting m of the parameters among a 1 , . . . , a n be equal to zero, and conclude that we can find reaction rate constants and total amounts such that the system has precisely 2m + 1 positive steady states for all m = 0, . . . , n.

n unstable steady states
In this subsection we show that, considering the 2n + 1 steady states ordered increasingly by their value of x 5 , then the steady states number 2, 4, . . . , 2n are unstable relative to the invariant subspaces described by the conservation laws (22) and (23).
Consider the Jacobian J f (x) of the function f in the ODE systemẋ = f (x). Since the system with n allosteric kinases competing for the same substrate has 4n + 2 variables and n + 1 conservation laws, the Jacobian of f in the ODE systemẋ = f (x) always has n + 1 zero eigenvalues. The eigenvectors of the remaining 3n + 1 eigenvalues (which could include zero) belong to the so-called stoichioimetric subspace, which has equations x i,1 + x i,2 + x i,3 + x i,4 = 0, i = 1, . . . , n, and dictate the dynamics around the steady state within this space. If the steady state is locally asymptotically stable relative to the invariant subspace described by the conservation laws, then the product of these 3n + 1 eigenvalues has sign (−1) 3n+1 . Therefore, if the sign of the product of these eigenvalues is (−1) 3n , then the steady state is necessarily unstable relative to the invariant subspace. We argue in the proof of the next theorem that this is the case for the steady states in even position 2, 4, . . . , 2n.
Theorem 15. The 2, 4, . . . , 2n-th steady states are unstable relative to the invariant subspace described by the conservation laws.
It follows from [8,Prop. 5.3] that the product of the 3n + 1 eigenvalues of the Jacobian with eigenvectors in the invariant subspace agrees with the determinant of the Jacobian of the function g : R 4n+2 → R 4n+2 , where (κ j,2 x j,3 + κ j,4 x j,4 ) − κ 7 x 6 .
Therefore, the steady state is unstable if the sign of the determinant of the Jacobian of g is (−1) 3n . We now apply the method described in [3] (see the Electronic Supplementary Material of the paper), to determine the sign of the determinant of the Jacobian of g from iterative eliminations. One can check that the expressions in (24)-(27) are obtained from iteratively eliminating x i,1 , . . . , x i,4 from the equations g 4(i−1)+1 (x) = · · · = g 4(i−1)+4 (x) = 0 respectively, which correspond to the conservation law with total amount K i tot together withẋ i,2 =ẋ i,3 = x i,4 = 0. For each i, we first eliminate x i,1 from g 4(i−1)+1 (x) = 0. We have that g 4(i−1)+1 (x) increases in x i,1 . We substitute next the value of x i,1 into g 4(i−1)+2 (x). The resulting expression decreases in x i,2 . We isolate x i,2 then from the equation g 4(i−1)+2 (x) = 0. We proceed in the same way to find x i,3 and x i,4 from g 4(i−1)+3 (x) = 0 and g 4(i−1)+4 (x) = 0. The functions g 4(i−1)+3 (x) and g 4(i−1)+4 (x) (after the appropriate substitutions are made) are decreasing in x i,3 and x i,4 respectively. Finally, x 6 is eliminated from g 4n+2 (x) = 0. Note that g 4n+2 (x) decreases in x 6 .
Let p(x 5 ) be the polynomial obtained after clearing denominators in (29). Then, by [3,Proposition 7], the sign of the determinant of the Jacobian of g at a steady state agrees with the sign of the derivative of p(x 5 ), p (x 5 ), times (−1) 3n+1 . Therefore, if p (x 5 ) is negative, then the corresponding steady state is unstable. Since p(0) is negative, the first real root of p(x 5 ) has positive derivative, and then the signs alternate. Therefore, the steady states corresponding to the 2, 4, . . . , 2n-th roots are unstable relative to the invariant subspace described by the conservation laws.

Allosteric kinases with several states
In this section we consider the case in which the allosteric kinase is specific to one substrate, but the kinase might have more than 2 states and transitions between all states occur. Let K i , i = 1, . . . , n denote the n states of the kinase. The general model is in this case: for all i, j = 1, . . . , n, i = j.
In order to analyse this system, we consider a reduced model. Then we use Theorem 1 and Theorem 2 to translate results on the number of positive steady states of the reduced model to the model of interest. Explicitly, we consider the cases where n = 3 and n = 4.
We proceed through three simplification steps to obtain the reduced network.