Enzyme sharing as a cause of multistationarity in signaling systems
EEnzyme sharing as a cause of multistationarity in signalingsystems
Elisenda Feliu and Carsten Wiuf ∗ Bioinformatics Research Centre, Aarhus University, C. F. Møllers All´e 8, DK-8000Aarhus, Denmark15/10/2018
Abstract
Multistationarity in biological systems is a mechanism of cellular decision making.In particular, signaling pathways regulated by protein phosphorylation display featuresthat facilitate a variety of responses to different biological inputs. The features thatlead to multistationarity are of particular interest to determine as well as the stabilityproperties of the steady states.In this paper we determine conditions for the emergence of multistationarity insmall motifs without feedback that repeatedly occur in signaling pathways. We derivean explicit mathematical relationship ϕ between the concentration of a chemical speciesat steady state and a conserved quantity of the system such as the total amount ofsubstrate available. We show that ϕ determines the number of steady states andprovides a necessary condition for a steady state to be stable, that is, to be biologicallyattainable. Further, we identify characteristics of the motifs that lead to multistationa-rity, and extend the view that multistationarity in signaling pathways arises frommultisite phosphorylation.Our approach relies on mass-action kinetics and the conclusions are drawn in fullgenerality without resorting to simulations or random generation of parameters. Theapproach is extensible to other systems. Keywords : steady state; kinase; stability; cross-talk; phosphorylation
Multistationarity (the existence of more than one steady state under particular biologicalconditions) in cellular systems can be seen as a mechanism for cellular decision mak- ∗ [email protected] a r X i v : . [ q - b i o . S C ] J u l Enzyme sharing and multistationarity ing. How it arises is therefore fundamental to the understanding of cell signaling, thatis, the communication of signals to regulate cellular activities and responses. Generally,cell signaling involves post-translational modifications of proteins, such as phosphorylation,acetylation, or methylation. These modifications change the state of a protein in a discretemanner, for example from an active to an inactive state.In eukaryotes, reverse phosphorylation is the most frequent form of protein modificationaffecting ∼
30% of all proteins in humans [1].
Kinases catalyze the transfer of phosphategroups to target proteins and phosphatases catalyze the reverse operation. After the com-pletion of the human genome project, genome analysis estimated the number of kinasesto ∼
500 [2], while the number of phosphatases is smaller by two thirds [1]. Two proteinphosphatases, PP-1 and PP-2A, account for the vast majority of all phosphatase activity[3] with more than 50 PP-1 targets being characterized [4].As a consequence, there is a substantial complexity in the interplay between enzymes(kinases and phosphatases) and substrates, exemplified by systems where protein substratesuse the same catalyzing enzymes (enzyme sharing), and systems where different enzymescatalyze the same reaction (enzyme competition). Competition and sharing are generalexamples of cross-talk between motifs.The aim of this work is to determine characteristics that lead to multistationarity.Following different modeling strategies it has already been shown that feedback in signalingnetworks as well as multisite phosphorylation can both account for multistationarity [5, 6,7]. We present a mathematical approach for analyzing the steady states of small systems.Our method leads to explicit conditions for when multistationarity occurs in terms of rateconstants and conserved total amounts of substrates and enzymes. Further the approachprovides means to study the stability of steady states.First we present the motifs that we analyze and then we develop the method to deter-mine multistationarity and to study stability. The paper concludes with some perspectivesand discussion.
We analyze the motifs shown in Fig. 5. The motifs are referred to as Motif (a)-(l) andprovide simple abstract representations of known cellular systems. Some examples moti-vating our choice of motifs are given in Table 1. A rich source of examples is found in thewell-studied MAPK cascades.To understand how multistationarity relates to enzyme usage we base our investigationon a motif that does not show multistationarity itself. Therefore we build the motifs froma one-site phosphorylation cycle which is monostable [16, 17, 18, 19] and shown in Motif
Enzyme sharing and multistationarity (a). A specific kinase (phosphatase) catalyzes phosphorylation (dephosphorylation) andall modifications can be reversed. In general, protein phosphoforms are denoted by S and P , Fig. 5. If one phosphoform is converted into another, an arrow is drawn and the enzyme( E or F ) catalyzing the reaction is indicated.Motifs (a)-(d) cover different possibilities for a one-site modification process. In Motif(b), the same enzyme catalyzes phosphorylation and dephosphorylation. Motifs (c) and(d) account for competition between kinases and/or phosphatases to catalyze the samemodification(s).In eukaryotes, phosphorylation of most proteins takes place in more than one site [20],potentially with different biological effects [21]. Combination of two one-site cycles intoa two-site sequential cycle yields three motifs: (e) all enzymes are different, (f) only onekinase but two phosphatases, and (g) one kinase and one phosphatase. By symmetry,Motifs Biological phenomena(b) A kinase acting also as phosphatase on the same substrate, e.g. HPrK/Pkinase-phosphatase in Gram positive bacteria [8].(c),(d) Several kinases and/or phosphatases acting on the same substrate, e.g. (i)Several kinases phosphorylate the alpha subunit eukaryotic initiation factor(eIF2 α ) at Ser51 [9]; (ii) The phosphatases MPK-1 and PTP-SL both modifyERK1 [10].(e) Multi-site phosphorylation by different kinases and phosphatases at each site,e.g. (i) Primed kinases, e.g. GSK-3 [11]; (ii) Akt1 is (de)activated throughthree-site sequential (de)phosphorylation by three different kinases (phos-phatases) [3].(f),(g) Multi-site phosphorylation with the same kinase and/or phosphatase respon-sible for all modifications, e.g. (i) Two-site phosphorylation of ERK catalyzedby MEK; (ii) Dephosphorylation of ERK2 catalyzed by DUSP6 [12].(h),(i) The same enzyme catalyzing the modification of two different substrates, e.g.the kinases ERK1, ERK2 and the kinase products of the p38 pathway catalyzephosphorylation of two substrates (the mitogen- and stress-activated proteinkinase (MSK) 1/2 and the MAP kinase signal-integrating kinase (MNK) 1/2)[13].(j),(k),(l) Cascades with several modification steps and substrates, e.g. (i) MAPKcascades; (ii) Protein kinase A (PKA) phosphorylates phosphorylase kinase,which in turn phosphorylates glycogen phosphorylase (with dephosphorylationcarried out by the same phosphatase, PP-1, in the two different layers) [14,Fig. 7.17],[15].Table 1: Cellular systems represented by Motifs (a)-(l) Enzyme sharing and multistationarity
Motif (f) represents as well a motif with one phosphatase but two kinases. We assumefor simplicity that both phosphorylation and dephosphorylation proceed in a sequentialand distributive manner [22]; that is, one site is (de)phosphorylated at a time in a specificorder.Motif (h) represents one-site modification of two substrates that share the same kinasebut use different phosphatases. This motif represents by symmetry also a system withshared phosphatase. If both the kinase and the phosphatase are shared, we obtain Motif(i). Finally, two one-site modification cycles can be combined in a cascade motif, where theactivated substrate of the first cycle acts as the kinase of the next. The interplay betweenenzymes is represented by three cascades: (j) dephosphorylation at each layer uses differentphosphatases, (k) the phosphatase is not layer specific, and (l) the kinase of the first layercatalyzes the modification in the second layer as well.
We assume that any modification S → S ∗ follows the classical Michaelis-Menten mecha-nism in which an intermediate complex Z is formed reversibly but dissociates into product One-site modification (a)(a) S S EF (b)(b) S S EE (c)(c) S S E , E F (d)(d) S S E , E F , F Two-site modification (e)(e) S S S E E F F (f)(f) S S S E EF F (g)(g) S S S E EF F
Modification of two substrates (h)(h) S S F P P F E (i)(i) S S P P EF Two-layer cascade (j)(j) S S E F P P F (k)(k) S S E F P P F (l)(l) S S EF P P F Figure 1:
Motifs composed of one or two one-site cycles. Motifs with purple label, and only these,admit multiple biologically meaningful steady states. S i and P i are substrates with i = 0 , , E, E , E denote kinases, and F, F , F phosphatases. In Motif (b) the kinaseand the phosphatase are the same enzyme. Enzyme sharing and multistationarity and enzyme G irreversibly: S + G a (cid:47) (cid:47) Z b (cid:111) (cid:111) c (cid:47) (cid:47) S ∗ + G The phosphate donor, generally ATP, is assumed to be in large constant concentrationand hence embedded into the rate constants. Imposing mass-action kinetics, the speciesconcentrations over time can be modeled by a system of polynomial differential equations.For example, in Motif (a) the equations are (here E also refers to the concentration of thekinase E , and similarly for the other species):˙ E = ( b E + c E ) X − a E ES ˙ X = − ( b E + c E ) X + a E ES ˙ F = ( b F + c F ) Y − a F F S ˙ Y = − ( b F + c F ) Y + a F F S ˙ S = b E X + c F Y − a E ES ˙ S = c E X + b F Y − a F F S , where X ( Y ) is the intermediate complex formed by the enzyme E ( F ) and the substrate S ( S ), and ˙ x denotes differentiation of x = x ( t ) with respect to time. For all motifs,there are conservation laws which define time-conserved quantities (total amounts), e.g.˙ E + ˙ X = 0 and so E = E + X is conserved. The total amounts are fixed by the initialconcentrations and determine the state space of the dynamical system. Motif (a) has threeconserved total amounts, namely F = F + Y and S = S + S + X + Y in addition to thatof E .The steady states of the system are solutions (potentially negative) to the polynomialequations obtained by setting all derivatives to zero with the constraints imposed by theconservation laws, once total amounts have been fixed. These laws imply that some steadystate equations are redundant, e.g. either ˙ E = 0 or ˙ X = 0 can be disregarded. We focuson the biologically meaningful steady states (BMSSs), that is, the steady states for whichall concentrations are non-negative (positive or zero). If at least two BMSSs exist for fixedtotal amounts, then the system is said to be multistationary.The specific form of the chemical reactions for Motifs (a)-(l) together with the cor-responding systems of differential equations are described in the Supplementary Material(SM), attached at the end of the main text. ϕ In this section we outline the procedure used to analyze the motifs. Details of the mathe-matical analysis are in the SM.The system of equations describing the steady states can be reduced substantially byelimination of variables [7, 23]. For the motifs considered here, elimination of variablesimplies that the steady states are characterized by a relation S = ϕ ( Y ) between theconcentration of one of the species, typically an intermediate complex Y , and the totalamount of a substrate S . The concentrations of the other species are given in terms of Y , Enzyme sharing and multistationarity usually as ratios of polynomials in Y . By imposing all concentrations to be non-negative, Y is restricted to a set Γ of possible values. Further, for any S ≥
0, there is at least one
BMSS, that is, S = ϕ ( Y ) for some Y in Γ. The function ϕ is continuous and differentiablein Γ and depends on the rate constants and the total amounts, except for S .The number of BMSSs can be found from the analysis of ϕ . If ϕ is strictly increasingor decreasing in Γ, ϕ is one-to-one and hence, for a given total amount S , there is acorresponding unique Y at steady state. Consequently, multistationarity cannot occur(Fig. 6a).Figs. 6b and 6c show situations where multistationarity occurs. If ϕ has increasing anddecreasing parts, or if Γ is not connected, then Y (cid:54) = Y with ϕ ( Y ) = ϕ ( Y ) = S mightexist. Hence, there are at least two BMSSs with the same S .These two figures represent substantially different switch responses. In Fig. 6c, there isonly one BMSS for low S . An increase of S to S max causes the system to switch to a ‘high’steady state (high Y ) under the assumption that the green steady states in the figure arestable. If S is decreased again to S min , then the system switches back to a ‘low’ steadystate. In Fig. 6b, there is one BMSS for low S . An increase of S keeps the system in thefirst branch of ϕ and thus it will behave as a monostationary system.Interestingly, the derivative ϕ (cid:48) ( Y ) of ϕ ( Y ) provides means to determine whether somesteady states are unstable. Unstable steady states are unattainable under biological con-ditions. Specifically, we find that either the regions in which ϕ is increasing or those inwhich it is decreasing must correspond to unstable steady states, see Section 5.In summary, the function ϕ determines whether multiple BMSSs exist and encodesinformation about the stability of steady states. In the next section we analyze ϕ for Motifs(e) and (f). We show how enzyme sharing in a two-site cycle (f) leads to multistationarity, s y YS (a) S min s y y y YS (b) S max S min s y y y YS (c) Figure 2:
Possible shapes of ϕ in Γ (colored regions: magenta=unstable BMSSs; green=(possible)stable BMSSs). (a) ϕ is increasing and for any s , there is one BMSS ( y ) such that ϕ ( y ) = s .(b) Γ consists of two disconnected regions. For s < S min there is one BMSS; for s = S min thereare precisely two; and for s > S min there are three; ϕ is also defined in the white region but someconcentrations become negative. (c) ϕ is in part decreasing, in part increasing. For S min < s < S max there are three BMSSs; for s = S min or s = S max , there are two; and for s < S min or s > S max ,there is one. Enzyme sharing and multistationarity as opposed to a two-site cycle with different enzymes (e). A detailed analysis of all motifsis given in the SM.
Motifs (a)-(e), (h), and (j) have exactly one BMSS for any choice of rate constants andtotal amounts. In all cases, the function ϕ is increasing in Γ. The procedure is very similarin all cases and is thus only illustrated for Motif (e). We take some effort in explaining thedetails as the procedure might have general applicability.Motif (e) consists of three phosphoforms of the substrate, S , S , S , with subscriptindicating the number of phosphorylated sites. The chemical reactions of the system are: S + E a ,E (cid:47) (cid:47) X b ,E (cid:111) (cid:111) c ,E (cid:47) (cid:47) S + E S + E a ,E (cid:47) (cid:47) X b ,E (cid:111) (cid:111) c ,E (cid:47) (cid:47) S + E S + F a ,F (cid:47) (cid:47) Y b ,F (cid:111) (cid:111) c ,F (cid:47) (cid:47) S + F S + F a ,F (cid:47) (cid:47) Y b ,F (cid:111) (cid:111) c ,F (cid:47) (cid:47) S + F We denote the inverse of the Michaelis-Menten constants of E i by κ i,E = a i,E / ( b i,E + c i,E )and of F i by κ i,F = a i,F / ( b i,F + c i,F ). The ratio of the catalytic constants of phosphataseand kinase is denoted by µ i = c i,F /c i,E .The system has five conserved total amounts, which are assumed to be positive: Fourfor enzymes, E i = E i + X i and F i = F i + Y i ( i = 1 , S = S + S + S + X + X + Y + Y . The steady state equations can be rewritten as X = κ ,E E S Y = κ ,F F S X = µ Y (1) X = κ ,E E S Y = κ ,F F S X = µ Y . The last column gives X i in terms of Y i . The total amounts E i , F i give E i , F i in terms of Y i as well: E i = E i − µ i Y i , F i = F i − Y i . Further, if E i , F i > E i = 0 or F i = 0cannot be solutions to Eq. (26). It follows that the concentrations E i , F i are positive if andonly if Y i is in Γ i = [0 , ξ i ) with ξ i = min( F i , E i /µ i ).We further isolate S , S from the first row in Eq. (26) and S from the second andobtain S = µ Y κ ,E ( E − µ Y ) , S i = Y i κ i,F ( F i − Y i ) (2)for i = 1 ,
2. Then, S , S (resp. S ) are non-negative increasing continuous functions of Y in Γ (resp. Y in Γ ). The remaining equation, X = κ ,E E S , gives Y in terms of Y : Y = f ( Y ) = κ ,E E Y µ ( κ ,F ( F − Y ) + κ ,E Y ) . (3) Enzyme sharing and multistationarity
The function f is non-negative increasing and continuous in Γ . Further, for Y to be inΓ , it is required that Y is in Γ = [0 , ξ ) ⊆ Γ with ξ = min( ξ , f − ( ξ )).Finally, using Eq. (29), we find that X and S are increasing functions of Y in Γ.Therefore, using the formulas above, all concentrations at steady state are non-negative ifand only if Y is in Γ. We conclude that the BMSSs of the system satisfy S = S + S + S + X + X + Y + Y = ϕ ( Y )for Y in Γ. Since ϕ is a sum of increasing continuous functions in Y , then so is ϕ .Additionally, ϕ (0) = 0 and ϕ ( Y ) tends to infinity as Y tends to ξ . Thus, ϕ has the formin Fig. 6a with a unique Y for any given S , that is, there is one BMSS. We consider a two-site phosphorylation system with one kinase but different phosphatasesfor each phosphoform, as shown in Motif (f). Multistationarity has been observed numer-ically in this system [6]. The system reduces to Motif (e) by setting E = E and we usethe notation introduced previously. The conservation laws are the same with the exceptionthat there is only one kinase law, E = E + X + X . Define ξ i = min( F i , E/µ i ) andΓ i = [0 , ξ i ).The system of equations to be solved is similar to Eq. (26) with E = E i . Thus, we startby writing X i , E, F i as functions of Y , Y . Since E, F i must be positive at any BMSS, werequire 0 ≤ Y i < F i and µ Y + µ Y < E . For these values we obtain S = µ Y κ ,E ( E − µ Y − µ Y ) , S i = Y i κ i,F ( F i − Y i )for i = 1 ,
2, which are non-negative increasing continuous functions of Y i . Using X = κ ,E ES , we obtain Y as a non-negative continuous function of Y in Γ : Y = f ( Y ) = κ ,E ( E − µ Y ) Y µ ( κ ,F ( F − Y ) + κ ,E Y ) . This function resembles that in Eq. (29) except for the quadratic term in the numerator,which is a consequence of the conservation law for E involving both Y and Y . Further, f might not be increasing for all Y .Let Γ = { Y ∈ Γ , such that f ( Y ) ∈ Γ } . Using the formulas derived above, all con-centrations at steady state are non-negative if and only if Y is in Γ. Hence, for anyBMSS, S = S + S + S + X + X + Y + Y = ϕ ( Y )with Y in Γ. The function ϕ is continuous with ϕ (0) = 0 but Γ might not be a connectedinterval. Enzyme sharing and multistationarity (cid:45) (cid:45) F ≤ NN < F ≤ f ( α ) f ( α ) < F < MM ≤ F Figure 3:
The function ϕ for Motif (f) for different values of F and fixed E > µ F and Λ > N = ( E − µ F ) /µ . The values M and α depend on E , F and the rate constants (see text).The set Γ is disconnected only when N < F ≤ f ( α ). For large F , ϕ approaches the black line.The vertical bars mark the boundary of Γ. Define Λ = (1 + κ ,E /κ ,F ) µ F − E . If Λ ≤
0, then f is an increasing function inΓ and we conclude that there is exactly one BMSS. If Λ >
0, then f has a unique localmaximum for some α in Γ and all cases in Fig. 6 can occur. By varying the value of F while keeping the other constants fixed, we obtain (Fig. 3):(i) F ≤ ( E − µ F ) /µ (orange): Γ = [0 , α ) with f ( α ) = F . The function f , andthus ϕ , are increasing and there is one BMSS (Fig. 6a).(ii) ( E − µ F ) /µ < F ≤ f ( α ) (green): Γ = [0 , α ) ∪ ( α , ξ ) with α ≤ α ≤ α and f ( α ) = f ( α ) = F . Hence, f is increasing in [0 , α ), decreasing in ( α , ξ ) andmultistationarity occurs (Fig. 6b).When f ( α ) < F there is an M such that:(iii) f ( α ) < F < M (purple): Γ = [0 , ξ ). The function ϕ has a decreasing part andmultistationarity occurs (Fig. 6c).(iv) M ≤ F (blue): Γ = [0 , ξ ). The function ϕ is increasing and there is one BMSS(Fig. 6a). Motifs (f), (g), (i), (k) and (l) exhibit multistationarity for some choices of total amountsand rate constants (Fig. 4). The regions for which multistationarity occurs are detailed inthe SM. In Motifs (i), (k) and (l) it appears only as in Fig. 6c, while in Motifs (f) and (g)both the forms in Fig. 6b and Fig. 6c occur.
It is remarkable that in Motifs (f), (k) and (l) multistationarity occurs for any set ofrate constants and depends only on the initial conditions (that is, the total amounts).Thus, multistationarity can occur in these systems independently of the specific kinetics.In contrast, multistationarity in Motif (i) depends on the rate constants and hence not allkinetics exhibit multistationarity. The same appears to be the case for Motif (g) [24, 25].The common characteristic of these motifs is that a single enzyme is responsible forcatalyzing two different substrate modifications, which at the same time are linked (Fig. 4).Indeed, in Motifs (f) and (g) the substrates are linked through S , which is a modified aswell as an unmodified substrate for the shared enzyme E . For the Motifs (k) and (l), thelink is given by S which is a modified substrate and a kinase, and the common enzymesare F and E , respectively. In Motif (i) the kinase E is common and the phosphatase F provides the link (or vice versa ). In contrast, in Motif (h) an enzyme is responsible fortwo different modifications, but there is no link between the two substrates. Consequently,multistationarity cannot be observed.Multistationarity can arise from two opposing dynamics acting on the same substrate Two-site modification . Shared enzyme: E ; Link: S (f)(f) S S S E E F F F (cid:29) E , F < M (g)(g) S S S E E
F Fµ F > E > µ F Two substrates . Shared enzyme: E ; Link: F (i)(i) S S P P EF µ F > E > µ Fµ F > E > µ F ( µ − µ )( µ δ η − µ δ η ) < Two-layer cascade . Shared enzyme: E ; Link: S (k)(k) S S E F P P FS large, F (cid:29) E (l)(l) S S EF P P F E (cid:29) µ F , F large Figure 4:
Conditions for multistationarity are given. The shared enzyme is marked with a coloredsquare; the link is marked with a colored circle; predominant modifications are marked in bold.The symbol (cid:29) is short for ‘much larger’. (Fig. 4). For example, in Motif (f), if F is much bigger than E and F < M , then thereare multiple BMSSs. Thus, since the amount of phosphatase in the first cycle is muchlarger than the amount of kinase, S is pushed towards the unmodified form S , while inthe second cycle S is driven towards the fully modified form S (because F < M ).In Motif (i), provided the conditions on the parameters are fulfilled (Fig. 4), multista-tionarity occurs if either µ F > E > µ F or µ F < E < µ F . It implies that in one cyclethe phosphatase ‘wins’, while in the other the kinase does. BMSSs are defined as steady states for which all concentrations are non-negative. However,a steady state is biologically attainable only if it is (asymptotically) stable, that is, nearbytrajectories are attracted to it. We show here for our motifs that if ϕ (cid:48) ( Y ) < ϕ ( Y ) = S , then it is unstable. For a system of ordinary differential equations in R m , a steady state z is asymptotically sta-ble if all eigenvalues of the Jacobian evaluated at z have negative real parts [26, Thm. 1.1.1].Since the Jacobian is a real matrix, the complex eigenvalues come in pairs of conjugatesand their product is a positive number. If m is odd and all eigenvalues have negative realparts, their product, and hence the determinant of the Jacobian, must be negative. If m is even and z stable, then the product of the eigenvalues must be positive. Thus, the signof the determinant of the Jacobian provides a necessary condition for a steady state to bestable and a sufficient condition for it to be unstable.For x = ( x , . . . , x n ) let x ( j ) = ( x , . . . , (cid:98) x j , . . . , x n ) ( x with x j removed). We makethe following observation (see SM for a proof): Let f = ( f , . . . , f n ) : Ω ⊆ R n → R n be adifferentiable function defined on an open set Ω and z such that f ( z ) = 0. Assume that x j can be eliminated from the equation f i = 0 in a neighborhood around z ; that is, there existsa differentiable function ψ : Ω ( j ) ⊆ R n − → R , z ( j ) in Ω ( j ) = { x ( j ) | x ∈ Ω } , such that x j = ψ ( x ( j ) ) if f i ( x ) = 0. Define ¯ f : Ω ( j ) → R n − by ¯ f k ( x ) = f k ( x , . . . , x j − , ψ ( x ) , x j , . . . , x n − )for all k (cid:54) = i and let ¯ J denote the associated Jacobian. Then, the determinant of theJacobian of f at z satisfies( − i + j ∂f i ∂x j ( z ) det( ¯ J ( z ( j ) )) = det( J ( z )) . (4) The relation between the sign of the determinant of the Jacobian and stability, togetherwith Eq. (73), lead to a criterion to detect unstable steady states. For each motif, let x = ( x , . . . , x n ) be the species concentrations, ˙ x i = h i ( x ) the differential equations and A = g ( x ) , . . . , A c = g c ( x ) the equations for the total amounts. We choose the order ofthe species such that x i , i = 1 , . . . , c , can be isolated from A i = g i ( x ) and the steady stateequation ˙ x i = 0 becomes redundant. For fixed total amounts, A , . . . , A c , the steady statesare the solutions to the system f ( x ) = 0 of n equations in n variables with f i ( x ) = g i ( x ) − A i for i = 1 , . . . , c and f i ( x ) = h i ( x ) for i = c + 1 , . . . , n .Let J ( z ) denote the Jacobian of f at z . In the SM we prove: If z is a steady state,that is, f ( z ) = 0, and either (i) n − c is even and det( J ( z )) < n − c is odd anddet( J ( z )) >
0, then z is unstable. The proof relies on the observation made about theeigenvalues and Eq. (73).The function ϕ of our motifs is derived through successive elimination of variablesprecisely from the system of equations f ( x ) = 0. Using Eq. (73), the sign of det( J ( z )) ata steady state z can be traced back from the sign of the derivative of ϕ (the Jacobian of asystem with one equation) by considering the equation number ( i ), the equation variable( j ), and the sign of ∂f i /∂x j after each elimination.To exemplify the procedure, consider Motif (f), where n = 10 and c = 4. The systemis (see SM for details): f ( x ) = E + X + X − E f ( x ) = X − κ ,E ES f ( x ) = F + Y − F f ( x ) = X − µ Y f ( x ) = F + Y − F f ( x ) = X − µ Y f ( x ) = S + S + S + X + X + Y + Y − S f ( x ) = Y − κ ,F F S f ( x ) = X − κ ,E ES f ( x ) = Y − κ ,F F S with species x = ( E, F , F , S , X , X , S , S , Y , Y ) . The function ϕ is f in terms of Y after successive eliminations. Let (cid:15) k = ± − ) or not (+) after the k -th elimination. Then,sign( ϕ (cid:48) ( Y )) = ( (cid:81) k (cid:15) k ) sign(det( J ( z ))), where z is the steady state with z = Y .The order and sign of the eliminations are shown in Table 2. We find that (cid:81) k (cid:15) k = 1,implying that the sign of ϕ (cid:48) ( Y ) agrees with the sign of the determinant of the Jacobianof f evaluated at the corresponding steady state. Since n − c = 6 is even, we concludethat the values of Y for which ϕ is decreasing, that is, ϕ (cid:48) ( Y ) <
0, correspond to unstablesteady states. Further, it follows that unstable points come between other steady statesthat presumably are stable.
The Routh-Hurwicz criterion [27] gives sufficient and necessary conditions for the Jacobianto have all eigenvalues with negative real parts. Thus, the (asymptotic) stability of a steady k Elimination Behavior a (cid:15) k b k Elimination Behavior (cid:15) k f , E ) (1 , , +) + 6 ( f , S ) (2 , , − ) +2 ( f , F ) (1 , , +) + 7 ( f , S ) (2 , , − ) +3 ( f , F ) (1 , , +) + 8 ( f , S ) (2 , , − ) +4 ( f , X ) (4 , , +) + 9 ( f , Y ) (2 , , − ) +5 ( f , X ) (4 , , +) + a ( i, j, σ ) indicates that i, j are the indices of the equation and variable iteratively beingeliminated and σ corresponds to whether ¯ f i ( f i after substitution of the previous eliminations) isincreasing ( σ = +) or decreasing ( σ = − ) as a function of x j . b obtained as obtained as σ ( − i + j . Table 2: Elimination of variables for Motif (f). After each elimination the system ¯ f isrewritten to correctly determine the sign of ∂ ¯ f i /∂x j before the next eliminationstate can be determined by this criterion. For the Motifs (a)-(e) and (h) the criterion isfulfilled and the unique BMSS is asymptotically stable. We have not been able to determinethis for Motif (j). We have investigated small motifs without feedback that account for cross-talk, enzymecompetition, sharing and specificity in post-translational modification systems and deter-mined some features that lead to multistationarity in signaling pathways.Bistability, and generally multistability, in biological systems is seen as a mechanismof cellular decision making. Compared to systems with a single steady state, the presenceof multiple stable steady states provides a possible switch between different responses andincreased robustness with respect to environmental noise. Our study has been driven bythe observation that biological systems deviate from a one-to-one correspondence betweenenzymes and the modifications they catalyze. This phenomenon, known as cross-talk andenzyme sharing, can cause multistationarity and hence be essential for regulating signalingsystems.Our work extends the view of multistationarity as arising from multisite phosphoryla-tion [7] to the view that multistationarity is driven by a single enzyme that catalyzes linkedsubstrates. Two opposing dynamics acting on the same substrate is a recurrent charac-teristic of multistationarity. These observations await a precise mathematical formulationand an investigation of its generality.
Our approach is conceptually simple and reduces to the study of analytical propertiesof a function ϕ that relates a conserved total amount and the concentration of a speciesat steady state. The graph ( ϕ ( Y ) , Y ) can be seen as a bifurcation diagram with oneparameter, S . When monostationarity occurs, analysis of ϕ is quite straightforward, whilea more in-depth analysis is required when multistationarity occurs. An advantage of thisapproach is that unstable steady states are readily detected from the form of ϕ .The existence of ϕ is not guaranteed in general. For larger systems, a detailed studyindependently of the rate constants cannot be pursued. However, we have shown that forcascades of arbitrary length (extentions of Motif (j)), the function ϕ exists and propertiesof the full cascade can be derived from properties of the building block, the one-site cycle[23]. We are currently working on extending this approach to other systems.There are some mathematical characterizations of monostationarity in chemical reac-tion networks that relate to our work. The theory of monotone systems [28] characterizessystems in which there is only one BMSS and at the same time gives conditions for globalstability of the BMSS. However, a condition for the theory to be applicable here is that nospecies take part in more than two reactions [29]. This condition is only fulfilled for Motif(a) and hence the theory cannot be applied here.The only general theory of applicability is that of injective systems [30, 31, 32]. Themotifs that do not allow multiple steady states are in fact injective in the sense of [30]when modeled as a continuous flow stirred tank reactor . This fact implies that at most one steady state exists [33]. The advantage of this theory is that monostationarity is derived formore general kinetics than mass-action [34]. However, when restricted to mass-action, ourapproach is as simple as checking for injectivity and provides in addition simple rationalfunctions that enable further comprehensive studies of variation in species concentrations,such as stimulus-response curves and signal amplification [23, 35]. Acknowledgements.
EF is supported by a postdoctoral grant from the “Ministerio deEducaci´on” of Spain and the project MTM2009-14163-C02-01 from the “Ministerio deCiencia e Innovaci´on”. CW is supported by the Lundbeck Foundation, Denmark. NeilBristow, Freddy Bugge Christiansen and Michael Knudsen are thanked for commenting onthe manuscript.
References
1. Cohen, P. 1989 The structure and regulation of protein phosphatases.
Annu. Rev.Biochem. , 453–508.2. Manning, G., Whyte, D. B., Martinez, R., Hunter, T., Sudarsanam, S. 2002 Theprotein kinase complement of the human genome. Science , 1912–1934.
3. Xiao, L., Gong, L. L., Yuan, D., Deng, M., Zeng, X. M., Chen, L. L., Zhang, L.,Yan, Q., Liu, J. P., Hu, X. H., et al.
CellDeath Differ. , 1448–1462.4. Cohen, P. T. 2002 Protein phosphatase 1–targeted in many directions. J. Cell Sci. , 241–56.5. Kapuy, O., Barik, D., Domingo Sananes, M. R., Tyson, J. J., Nov´ak, B. 2009 Bistabilityby multiple phosphorylation of regulatory proteins.
Prog. Biophys. Mol. Biol. (1-3),47–56.6. Markevich, N. I., Hoek, J. B., Kholodenko, B. N. 2004 Signaling switches and bistabilityarising from multisite phosphorylation in protein kinase cascades.
J. Cell Biol. ,353–359.7. Thomson, M., Gunawardena, J. 2009 Unlimited multistability in multisite phosphory-lation systems.
Nature , 274–277.8. Chaptal, V., Vincent, F., Gueguen-Chaignon, V., Monedero, V., Poncet, S., Deutscher,J., Nessler, S., Morera, S. 2007 Structural analysis of the bacterial HPr ki-nase/phosphorylase V267F mutant gives insights into the allosteric regulation mecha-nism of this bifunctional enzyme.
J. Biol. Chem. , 34952–34957.9. de Haro, C., Mendez, R., Santoyo, J. 1996 The eIF-2alpha kinases and the control ofprotein synthesis.
FASEB J. , 1378–1387.10. Shaul, Y. D., Seger, R. 2007 The mek/erk cascade: from signaling specificity to diversefunctions. Biochim. Biophys. Acta (8), 1213–26.11. Cohen, P. 2000 The regulation of protein function by multisite phosphorylation–a 25year update.
Trends Biochem. Sci. (12), 596–601.12. Ferrell, J. E., Bhatt, R. R. 1997 Mechanistic studies of the dual phosphorylation ofmitogen-activated protein kinase. J. Biol. Chem. (30), 19008–16.13. Keshet, Y., Seger, R. 2010 The map kinase signaling cascades: a system of hundredsof components regulates a diverse array of physiological functions.
Methods Mol. Biol. , 3–38.14. Fell, D. 1997
Understanding the control of metabolism Frontiers in Metabolism , vol-ume 2. Portland Press 1st edition.15. Cohen, P. 1992 Signal integration at the level of protein kinases, protein phosphatasesand their substrates.
Trends Biochem. Sci. , 408–413.
16. Goldbeter, A., Koshland, D. E. 1981 An amplified sensitivity arising from covalentmodification in biological systems.
Proc. Natl. Acad. Sci. USA , 6840–6844.17. Goldbeter, A., Koshland, D. E. 1984 Ultrasensitivity in biochemical systems controlledby covalent modification. Interplay between zero-order and multistep effects. J. Biol.Chem. , 14441–14447.18. Bluthgen, N., Bruggeman, F. J., Legewie, S., Herzel, H., Westerhoff, H. V., Kholo-denko, B. N. 2006 Effects of sequestration on signal transduction cascades.
FEBS J. , 895–906.19. Salazar, C., H¨ofer, T. 2006 Kinetic models of phosphorylation cycles: a systematicapproach using the rapid-equilibrium approximation for protein-protein interactions.
Biosystems , 195–206.20. Olsen, J. V., Blagoev, B., Gnad, F., Macek, B., Kumar, C., Mortensen, P., Mann, M.2006 Global, in vivo, and site-specific phosphorylation dynamics in signaling networks. Cell , 635–648.21. Wu, R. C., Qin, J., Yi, P., Wong, J., Tsai, S. Y., Tsai, M. J., O’Malley, B. W. 2004Selective phosphorylations of the SRC-3/AIB1 coactivator integrate genomic reponsesto multiple cellular signaling pathways.
Mol. Cell , 937–949.22. Salazar, C., H¨ofer, T. 2009 Multisite protein phosphorylation–from molecular mecha-nisms to kinetic models. FEBS J. , 3177–3198.23. Feliu, E., Knudsen, M., Andersen, L. N., Wiuf, C. 2011 An algebraic approach tosignaling cascades with n layers. Bull. Math. Biol. doi:10.1007/s11538-011-9658-0.24. Conradi, C., Flockerzi, D., Raisch, J. 2008 Multistationarity in the activation of aMAPK: parametrizing the relevant region in parameter space.
Math. Biosci. ,105–131.25. Wang, L., Sontag, E. D. 2008 On the number of steady states in a multiple futile cycle.
J. Math. Biol. , 29–52.26. Wiggins, S. 2003 Introduction to applied nonlinear dynamical systems and chaos Textsin Applied Mathematics , volume 2. New York: Springer-Verlag 2nd edition.27. Hurwitz, A. 1996 ¨Uber die Bedingungen, unter welchen eine Gleichung nur Wurzelnmit negativen reellen Theilen besitzt.
In Stability theory (Ascona, 1995) Internat. Ser.Numer. Math. , volume 121 pp. 239–249. Basel: Birkh¨auser. Reprinted from Math.Ann.
28. Smith, H. L. 1995
Monotone dynamical systems Mathematical Surveys and Mono-graphs , volume 41. Providence, RI: American Mathematical Society. An introductionto the theory of competitive and cooperative systems.29. Angeli, D., De Leenheer, P., Sontag, E. 2010 Graph-theoretic characterizations ofmonotonicity of chemical networks in reaction coordinates.
J. Math. Biol. , 581–616.30. Craciun, G., Feinberg, M. 2005 Multiple equilibria in complex chemical reaction net-works. I. The injectivity property. SIAM J. Appl. Math. (5), 1526–1546 (electronic).31. Craciun, G., Feinberg, M. 2006 Multiple equilibria in complex chemical reaction net-works. II. The species-reaction graph. SIAM J. Appl. Math. (4), 1321–1338 (elec-tronic).32. Craciun, G., Tang, Y., Feinberg, M. 2006 Understanding bistability in complex enzyme-driven reaction networks. Proc. Natl. Acad. Sci. USA , 8697–8702.33. Craciun, G., Feinberg, M. 2006 Multiple equilibria in complex chemical reaction net-works: extensions to entrapped species models.
Syst Biol (Stevenage) , 179–186.34. Banaji, M., Donnell, P., Baigent, S. 2007 P matrix properties, injectivity, and stabilityin chemical reaction systems. SIAM J. Appl. Math. (6), 1523–1547.35. Feliu, E., Knudsen, M., Wiuf, C. 2011 Signaling cascades: consequences of varyingsubstrate and phosphatase levels . Advances in Systems Biology. Springer Verlag. Inpress.
Supplementary Material
A Introduction
In this Supplementary Material (SM) we provide details of the analysis of multistationaritygiven in the main text, as well as proofs of the results mentioned there. We go throughall motifs and the reader will easily see that many arguments are repeated as the motifsshare common structures. We have tried to keep the analysis of each motif independentof the analyses of the other motifs. However, the motifs progress from simple one-sitemodification cycles to more complex motifs and the line of thought transpires most easilyfrom the simple motifs. It is therefore advisable to study the simple motifs before the morecomplicated motifs. To keep the SM self-contained, Figures 1 and 2 of the main text arereproduced here.Motif (g) (a futile cycle with two parts) has been studied extensively in the literature.Motif (j) (a linear cascade with two layers) was studied for arbitrary length in [23], wherewe showed that the system admits only one steady state. It is briefly covered here forcompleteness.
One-site modification (a)(a) S S EF (b)(b) S S EE (c)(c) S S E , E F (d)(d) S S E , E F , F Two-site modification (e)(e) S S S E E F F (f)(f) S S S E EF F (g)(g) S S S E EF F
Modification of two substrates (h)(h) S S F P P F E (i)(i) S S P P EF Two-layer cascade (j)(j) S S E F P P F (k)(k) S S E F P P F (l)(l) S S EF P P F Figure 5:
Motifs covered in this work. Figure 1 of the main text. s y YS (a) S min s y y y YS (b) S max S min s y y y YS (c) Figure 6:
Possible shapes of ϕ in Γ (colored regions: magenta=unstable BMSSs; green=(possible)stable BMSSs). (a) ϕ is increasing and for any s , there is one BMSS ( y ) such that ϕ ( y ) = s .(b) Γ consists of two disconnected regions. For s < S min there is one BMSS; for s = S min thereare precisely two; and for s > S min there are three; ϕ is also defined in the white region but someconcentrations become negative. (c) ϕ is in part decreasing, in part increasing. For S min < s < S max there are three BMSSs; for s = S min or s = S max , there are two; and for s < S min or s > S max ,there is one. Figure 2 in the main text. A.1 Notation and preliminaries
A continuous and differentiable function with continuous derivative is said to be C . Wedenote by R + the set of positive real numbers and by R + the set of non-negative realnumbers. A rational function is a quotient of two polynomials. An increasing (decreasing)function f fulfills f ( x ) < f ( y ) ( f ( x ) > f ( y )) for x < y , i.e. we take increasing (decreasing)to mean strictly increasing (decreasing). The notation x ∈ A means that x belongs to theset A .We make frequent use of the Implicit Function Theorem in two dimensions to relatetwo variables to each other and to find derivatives of implicit functions. Details about theImplicit Function Theorem can be found in text books on functional analysis.For x = x ( t ) a real function of t , we denote by ˙ x the derivative of x with respect to t , dx/dt . A.2 Constants and variables
We consider the rate constants to be fixed and positive, i.e. in R + . The constants are a ∗∗ , b ∗∗ , c ∗∗ and those derived from these: κ ∗ , ∗ , the inverse of the Michaelis-Menten constantsand µ ∗∗ , the ratios of the catalytic constants of phosphatase and kinase. Here, to ease thepresentation, η ∗ denotes κ ∗ ,E for a kinase E and δ ∗ denotes κ ∗ ,F for a phosphatase F .The total amounts are likewise considered fixed and positive. Species concentrationsare considered variables of the system. For example, in X − ηES , η is a fixed unknownconstant and X, E and S are variables that depend on each other, e.g. X can be considereda function of E and S , X = X ( E, S ). That is, the differential equations describing thesystem is a set of polynomials with coefficients in the ring R [ a ∗∗ , b ∗∗ , c ∗∗ , S ∗ , P ∗ , E ∗ , F ∗ ] or R [ η ∗ , δ ∗ , µ ∗ , S ∗ , P ∗ , E ∗ , F ∗ ] (where some constants might not be present in all systems,e.g. P is not in Motif (a)) and variables being a finite list of species concentrations, E ∗ , F ∗ , S ∗ , P ∗ , X ∗ , Y ∗ (where some variables might not be present in all systems, e.g. P is not in Motif (a)). A.3 BMSSs and total amounts
We only consider the system at steady state, that is all differential equations are put to zeroand solved for the species concentrations under the constraints imposed by the conservationlaws. Only solutions to the system of equations for which all species concentrations arenon-negative, i.e. positive or zero, are of interest. These solutions are called BiologicallyMeaningful Steady States, or BMSSs. Also, we assume that the total amounts are positive (non-zero), i.e. enzymes and substrates are always present in the system in some form(e.g. in some phosphoform or intermediate complex).One can prove for each specific motif that if all total amounts are positive, then allspecies concentrations at a BMSS are positive as well (i.e. if a BMSS exists then all speciesconcentrations must be positive). It follows that if 1) the species concentrations are non-negative (i.e. positive or zero), 2) the total amounts are positive, and 3) the conservationlaws and the steady state equations are fulfilled, then the species concentrations constitutein fact a BMSS and hence are positive. This observation is frequently used in the following.
A.4 Method
For all motifs we follow the same procedure. We take the set of differential equationsdescribing the systems together with the equations for the total amounts (the conservationlaws) and solve for one variable, say an intermediate complex Y . Specifically, we chooseone equation for a total amount, e.g. S = S + S + X + Y , and use the other equations(differential equations and equations for total amounts) to provide expressions for thevariables at steady state in terms of Y , i.e. we find expressions such that S = S ( Y ), S = S ( Y ) and X = X ( Y ) are functions of Y . These functions only depend on Y , therate constants and the total amounts, excluding S . The functions are substituted into theexpression for S to obtain a function ϕ ( Y ) that relates Y to S , i.e. S = ϕ ( Y ). The analyticform of ϕ determines how many BMSSs the system has for a given set of rate constantsand total amounts. For example, if ϕ is increasing there can only be one Y correspondingto a given S . See Figure 6 for illustration.We allow Y to be zero in which case S = ϕ ( Y ) also is zero (it follows from the con-struction of ϕ for each motif). This is done for simplicity as we otherwise would have toconsider the limit of ϕ ( Y ) for Y → A.5 Technicalities
Some manipulations of the equations are left to the reader if they only involve standardelimination, rewriting and differentiation techniques. In general, it is a good idea to dothe calculations yourself (either by hand or in Mathematica
T M , Maple or similar) as manydetails are left out due to space and readability constraints. We have done all derivationsin Mathematica
T M and checked them manually. Proofs of propositions are collected in anappendix to keep the presentation simpler. Equation numbers are to the immediate rightof the equations and only equations that are used later are numbered.
B Steady states of the motifs
B.1 One-site phosphorylation cycles
Motif (a).
We recently used this motif as the building block of linear cascades andshowed that it admits only one BMSS [23]. We follow the approach taken in [23] andsummarize it below. Essentially this is the approach taken for all motifs.The chemical reactions of the motif are S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F The corresponding system of differential equations is:˙ S = b E X + c F Y − a E ES ˙ E = ( b E + c E ) X − a E ES ˙ F = ( b F + c F ) Y − a F F S ˙ S = c E X + b F Y − a F F S (5)˙ X = − ( b E + c E ) X + a E ES (6)˙ Y = − ( b F + c F ) Y + a F F S . (7)It follows that ˙ E + ˙ X = 0, ˙ F + ˙ Y = 0 and ˙ S + ˙ S + ˙ X + ˙ Y = 0. Hence, the conservationlaws are given by E = E + X, F = F + Y, S = S + S + X + Y. Due to the conservation laws some equations are redundant, for example ˙ F = 0 and ˙ Y = 0are the same equation. If fixed total amounts are given, we have to solve a system of 6equations in 6 variables consisting of the equations for the total amounts and, for instance,equations (5)-(7). From the latter equations, we obtain the equivalent system X − ηES = 0 , Y − δF S = 0 , X − µY = 0 (8)with constants η = a E / ( b E + c E ), δ = a F / ( b F + c F ) and µ = c F /c E . It follows that at steady state E = E − µY and F = F − Y . Note that if both E, F arepositive (i.e. non-zero), then we cannot have E = 0 or F = 0 at steady state (for exampleif E = 0, then X = 0 from (6), hence E = 0). Let ξ = min( F , E/µ ) and Γ = [0 , ξ ). The variables
E, F are both positive and Y is non-negative only if Y is in Γ. Further it followsfrom (8) that S = µYη ( E − µY ) , and S = Yδ ( F − Y ) . (Note that division by zero does not occur as E = E − µY and F = F − Y are both greaterthan zero.) These two functions are non-negative, increasing and C for Y ∈ Γ. When Y tends to ξ , one of them tends to + ∞ , while the other remains bounded.Let ∆ = µF − E , so that Γ = [0 , F ) if ∆ ≤
0, and Γ = [0 , E/µ ) if ∆ >
0. Using theconservation law for S we obtain: Result 1 (Motif (a)) . Let a one-site modification cycle be given with positive total amounts
S, E, F . Then the system has a unique BMSS. Specifically, the BMSS satisfies S = ϕ ( Y ) for Y in Γ = [0 , ξ ) , ξ = min( F , E/µ ) , where S = ϕ ( Y ) = µYη ( E − µY ) + Yδ ( F − Y ) + (1 + µ ) Y is an increasing rational C function which tends to infinity as Y tends to ξ and fulfills ϕ (0) = 0 . Remark 1.
Since ϕ is a rational function, the equation S = ϕ ( Y ) can be written inpolynomial form by elimination of denominators. In the present case p ( Y ) = µδY ( F − Y ) + ηY ( E − µY ) + ηδ ( (1 + µ ) Y − S )( F − Y )( E − µY ) , which is a third degree polynomial in Y . Note that p (0) < p ( ξ ) > p ( ζ ) < < ξ < ζ = max( F , E/µ ), and p ( Y ) tends to + ∞ as Y tends to + ∞ ; hence p ( Y ) hasthree positive roots. However, only the first root is in Γ, and it corresponds to the onlyBMSS of the system. In some systems, several positive roots are BMSSs. This remark isapplicable in all cases below whenever ϕ is a rational function. Remark 2.
We argued above that E and F are non-zero if E and F are positive. Thisensures that we do not divide by zero when constructing ϕ . For the remaining variableswe only need to ensure non-negativity, cf. Section A.3. For the other motifs, we make useof similar reasoning. Motif (b).
Consider the system in Motif (b) where the two catalyzing enzymes are thesame. The chemical reactions of the system are given by S + E a (cid:47) (cid:47) X b (cid:111) (cid:111) c (cid:47) (cid:47) S + E S + E a (cid:47) (cid:47) Y b (cid:111) (cid:111) c (cid:47) (cid:47) S + E The corresponding system of differential equations is: ˙ S = b X + c Y − a ES ˙ E = ( b + c ) X + ( b + c ) Y − a ES − a ES ˙ S = c X + b Y − a ES (9)˙ X = − ( b + c ) X + a ES (10)˙ Y = − ( b + c ) Y + a ES . (11)We find that ˙ E + ˙ X + ˙ Y = 0 and ˙ S + ˙ S + ˙ X + ˙ Y = 0. Hence, the conservation lawsare given by E = E + X + Y, S = S + S + X + Y. If the total amounts are given, we need to solve a system of 5 equations in 5 variablesconsisting of those for the total amounts and, for instance, equations (9)-(11). From thelatter equations we obtain an equivalent system given by X − ηES = 0 , Y − δES = 0 , X − µY = 0with constants η = a / ( b + c ), δ = a / ( b + c ) and µ = c /c . Note that if
E >
0, then E = 0 cannot be a solution of the steady state equations and thus we require E > E = E − ( µ + 1) Y . For E > Y ≥ Y must be inΓ = [0 , E/ (1 + µ )). Further, we find S = µYη ( E − ( µ + 1) Y ) , S = Yδ ( E − ( µ + 1) Y ) . These functions are continuous and increasing for Y in Γ. In addition, all concentrations S , S , X, E, F are non-negative as functions of Y if and only if Y ∈ Γ. Using the conser-vation law for S we obtain: Result 2 (Motif (b)) . Let a one-site modification cycle be given with one enzyme actingas kinase as well as phosphatase. Further, assume that the total amounts
S, E are positive.Then, the system has a unique BMSS.Specifically, the BMSS satisfies S = ϕ ( Y ) for Y in Γ = [0 , ξ ) , ξ = E/ (1 + µ ) , where S = ϕ ( Y ) = µYη ( E − ( µ + 1) Y ) + Yδ ( E − ( µ + 1) Y ) + (1 + µ ) Y is an increasing, rational C function which tends to infinity as Y tends to ξ and fulfills ϕ (0) = 0 . Motif (c).
Consider now a one-site modification cycle with two competing kinases. Thechemical reactions of the system are given by ( k = 1 , S + E k a Ek (cid:47) (cid:47) X kb Ek (cid:111) (cid:111) c Ek (cid:47) (cid:47) S + E k S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F The corresponding system of differential equations is: ˙ E k = ( b Ek + c Ek ) X k − a Ek E k S ˙ F = ( b F + c F ) Y − a F F S ˙ S = b E X + b E X + c F Y − a E E S − a E E S ˙ X k = − ( b Ek + c Ek ) X k + a Ek E k S (12)˙ Y = − ( b F + c F ) Y + a F F S (13)˙ S = c E X + c E X + b F Y − a F F S (14) with k = 1 ,
2. It follows that ˙ X k + ˙ E k = 0, ˙ Y + ˙ F = 0 and ˙ S + ˙ S + ˙ X + ˙ X + ˙ Y = 0,leading to the conserved total amounts ( k = 1 , E k = E k + X k , F = F + Y, S = S + S + X + X + Y. Therefore, if total amounts are given, the steady states of the system are solutions to asystem of 8 equations in 8 variables consisting of the equations for the total amounts (4equations) and, for instance, equations (12) for k = 1 ,
2, (13) and (14). From the latterequations we obtain an equivalent system given by X k − η k E k S = 0 , Y − δF S = 0 , c E X + c E X − c F Y = 0 (15)for k = 1 ,
2, with constants η k = a Ek / ( b Ek + c Ek ) and δ = a F / ( b F + c F ). Note that if E k , F >
0, then neither E k = 0 nor F = 0 are solutions of the steady state equations.Thus, E k , F > E k = E k − X k and F = F − Y . Hence, for E, F > Y ≥ ≤ X k < E k , 0 ≤ Y < F at BMSS. We find S = Yδ ( F − Y ) , (16)which is increasing in Y . Likewise, we obtain the following relation (with non-zero deno-minators) X η ( E − X ) = S = X η ( E − X ) , and it follows that X = φ X ( X ) = η E X η ( E − X ) + η X . This function is increasing and C for X in Γ = [0 , E ). Additionally, 0 ≤ X = φ X ( X ) Y < F .Since S is a non-negative, increasing C function of Y (for Y < F , (16)), it is also anon-negative, increasing function of X ∈ Γ. In conclusion, all steady state concentrationsare non-negative if and only if X ∈ Γ. Additionally, either S or S tend to infinity as X approaches ξ . Using the conservation law for S we obtain: Result 3 (Motif (c)) . Let a one-site modification cycle be given with two different competingkinases and one phosphatase. Further, assume that positive total amounts S, E , E , F aregiven. Then the system has a unique BMSS.Specifically, the BMSS satisfies S = ϕ ( X ) for X in Γ = [0 , ξ ) , where S = ϕ ( X ) = X η ( E − X ) + φ Y ( X ) δ ( F − φ Y ( X )) + φ X ( X ) + X + φ Y ( X ) is an increasing, rational C function which tends to infinity as X tends to ξ and fulfills ϕ (0) = 0 . Motif (d). Consider a one-site modification cycle with two competing kinases and twocompeting phosphatases. The chemical reactions of the system are given by S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F The corresponding system of differential equations is: ˙ X k = − ( b Ek + c Ek ) X k + a Ek E k S (17)˙ Y k = − ( b Fk + c Fk ) Y k + a Fk F k S (18)˙ S = b E X + c F Y − a E E S + b E X + c F Y − a E E S ˙ S = c E X + b F Y − a F F S + c E X + b F Y − a F F S (19) ˙ E k = ( b Ek + c Ek ) X k − a Ek E k S ˙ F k = ( b Fk + c Fk ) Y k − a Fk F k S with k = 1 , 2. It follows that ˙ X k + ˙ E k = 0, ˙ Y k + ˙ F k = 0 and ˙ S + ˙ S + ˙ X + ˙ X + ˙ Y + ˙ Y = 0,leading to the following conserved total amounts ( k = 1 , E k = E k + X k , F k = F k + Y k , S = S + S + X + X + Y + Y . Therefore, if total amounts are given, we have to solve a system of 10 equations in 10variables consisting of the equations for the total amounts (5 equations) and, for instance,equations (17), (18) for k = 1 , X k − η k E k S = 0 , Y k − δ k F k S = 0 , c E X + c E X − c F Y − c F Y F = 0 (20)for k = 1 , η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ).If E k , F k > 0, then E k , F k (cid:54) = 0 at steady state. As for the previous motif, we have E k = E k − X k and F k = F k − Y k ; hence 0 ≤ X k < E k and 0 ≤ Y k < F k is required forany BMSS. The concentration S can be expressed in two different ways as increasing C functions: As a function of X and as a function of X . When these two expressions areequated we obtain the relation (similar to the relation obtained for the previous motif) X = φ X ( X ) = η E X η ( E − X ) + η X . It is a non-negative, increasing C function for X ∈ [0 , E ) such that E = E − φ X ( X )also is a positive function of X ∈ [0 , E ). Note that φ X ( E ) = E is well-defined.Similarly, S can be expressed as increasing C functions of Y and of Y , respectively,which provide the relation Y = φ Y ( Y ) = δ F Y δ ( F − Y ) + δ Y . It is a non-negative, increasing C function for Y ∈ [0 , F ) and F = F − φ Y ( Y ) is apositive function of Y ∈ [0 , F ).Finally, the last relation in (20) gives c E φ X ( X ) + c E X = c F φ Y ( Y ) + c F Y . The left side is an increasing C function in X , the right side an increasing C functionin Y which tends to infinity as Y tends to infinity. Hence, there exists an increasing C function g ( X ) = Y defined on X ∈ [0 , E ) relating X to Y .In summary, the concentrations E , E , X , S are non-negative functions of X if andonly if X ∈ [0 , E ). The concentrations F , F , Y , S are non-negative if and only if Y ∈ [0 , F ). Hence, to ensure that all concentrations are non-negative, we require Y = g ( X ) < F . Since g is increasing, it is bounded above by g ( E ) (well-defined). Therefore,let Γ = [0 , E ) if g ( E ) ≤ F and Γ = [0 , g − ( F )) otherwise. Using the conservation lawfor S we obtain: (Motif (d)) . Let a one-site modification cycle with two competing kinases and twocompeting phosphatases be given. Further, assume that the total amounts S, E , E , F , F are positive. Then, the system has a unique BMSS.Specifically, the BMSS satisfies S = ϕ ( X ) for X in Γ = [0 , ξ ) , where S = ϕ ( X ) = X η ( E − X ) + g ( X ) δ ( F g ( X )) + φ X ( X ) + X + φ Y ( g ( X )) + g ( X ) is an increasing C function which tends to infinity as X tends to ξ and fulfills ϕ (0) = 0 . The function g is not rational, hence neither ϕ is rational. B.2 Two-site phosphorylation cycles Motif (e). First we consider a two-site phosphorylation system in which modifications arecarried out by different kinases and phosphatases for each phosphoform. For simplicity, weassume that both phosphorylation and dephosphorylation occur sequentially. The chemicalreactions of the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F The differential equations describing the system are: ˙ X k = a Ek E k S k − − ( b Ek + c Ek ) X k ˙ Y k = a Fk F k S k − ( b Fk + c Fk ) Y k , ˙ S = b E X + c F Y − a E E S ˙ E k = ( b Ek + c Ek ) X k − a Ek E k S k − (21)˙ F k = ( b Fk + c Fk ) Y k − a Fk F k S k (22)˙ S = c E X + b F Y − a F F S (23)˙ S = c E X + b E X + b F Y + c F Y − ( a F F + a E E ) S (24) with k = 1 , 2. We have ˙ X k + ˙ E k = 0, ˙ Y k + ˙ F k = 0 and ˙ S + ˙ S + ˙ S + ˙ X + ˙ X + ˙ Y + ˙ Y = 0,which lead to the following conserved total amounts ( k = 1 , E k = E k + X k , F k = F k + Y k , S = S + S + S + X + X + Y + Y . (25)Therefore, if total amounts are given, we have to solve a system of 11 equations in 11variables consisting of those in (25) and, for instance, equations (21)-(24). From the latterequations, we obtain X k = η k E k S k − , Y k = δ k F k S k , k = 1 , with constants η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). Equation (22) and (23) for k = 2 give c E X − c F Y = 0 . This relation together with equation (24), (21) for k = 2, and(22) for k = 1 give c E X − c F Y = 0 . Therefore, we have that X k = µ k Y k , µ k = c Fk /c Ek . Let ∆ k = µ k F k − E k , ξ k = min( F k , E k /µ k ) and Γ k = [0 , ξ k ). Note that if E k , F k > E k , F k (cid:54) = 0. Since E k = E k − µ k Y k and F k = F k − Y k at steady state,any BMSS must satisfy Y k ∈ Γ k to ensure non-negativity of X k , Y k and positivity of E k , F k as functions of Y k for each k .For fixed values of S , X , Y , the steady state values of the remaining variables satisfythe steady state equations of Motif (a) (a one-site phosphorylation cycle) with species S , S , X , Y , E , F and total amounts E , F and S − S − X − Y . Therefore, usingResult 1, the BMSSs of the system satisfy S = ϕ ( Y ) + ( S + X + Y )with Y ∈ Γ . Here ϕ denotes the function ϕ in Result 1.Using the second equality in (26) for k = 2 together with the conservation law for F ,we obtain ϕ ( Y ) = S + X + Y = Y δ ( F − Y ) + µ Y + Y , (27)which is a non-negative, increasing C function of Y ∈ Γ . Consequently, the BMSSs ofthe system satisfy the relation S = ϕ ( Y ) + ϕ ( Y ) , (28)for Y ∈ Γ and Y ∈ Γ . The right hand side is an increasing function in both variables.From the equation X = η E S together with S expressed as a function of Y (similarto that of S in equation (27)) we obtain X = η ( E − X ) S = η ( E − X ) Y δ ( F − Y ) . Rewriting this equation yields X = η E Y δ ( F − Y ) + η Y , and hence Y = f ( Y ) = η E Y µ δ ( F − Y ) + µ η Y . (29)The function f is an increasing function, which is non-negative and C for Y ∈ Γ . Addi-tionally, f is defined for Y = F with f ( F ) = E /µ . By substitution of f ( Y ) into (28), the BMSSs of the system satisfy S = ϕ ( Y ) = ϕ ( Y ) + ϕ ( f ( Y ))for Y ∈ Γ and f ( Y ) ∈ Γ . Since f is increasing, continuous and f (0) = 0, this condition isequivalent to Y ∈ Γ = [0 , ξ ) with ξ = min( ξ , f − ( ξ )) (note that f − ( ξ ) is well-defined).The function ϕ is a rational function, which is increasing and C in Γ and either ϕ or ϕ ◦ f tends to infinity as Y tends to ξ . Additionally, the BMSS concentrations of allother species derived from the formulas above are non-negative functions of Y if and onlyif Y ∈ Γ.Using the conservation law for S we obtain: Result 5 (Motif (e)) . Let a two-site phosphorylation cycle be given with different kinasesand phosphatases. Further assume that the total amounts E k , F k , S , k = 1 , are positive.Then, the system has a unique BMSS.Specifically, the BMSS satisfies S = ϕ ( Y ) for Y in Γ = [0 , ξ ) , where S = ϕ ( Y ) = ϕ ( Y ) + ϕ ( f ( Y )) is an increasing rational C function, which tends to infinity as Y tends to ξ and fulfills ϕ (0) = 0 . Motif (f ). Next, we consider a two-site phosphorylation system where phosphorylation iscatalyzed by the same kinase at both sites but dephosphorylation is catalyzed by differentphosphatases. Again, we assume sequential (de)phosphorylation. The chemical reactionsof the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + ES + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F The differential equations describing the system are the following: ˙ S = b E X + c F Y − a E ES ˙ Y k = − ( b Fk + c Fk ) Y k + a Fk F k S k ˙ E = ( b E + c E ) X + ( b E + c E ) X − ( a E S + a E S ) E ˙ S = c E X + b E X + b F Y + c F Y − ( a F F + a E E ) S (30)˙ S = c E X + b F Y − a F F S (31)˙ X k = − ( b Ek + c Ek ) X k + a Ek ES k − (32)˙ F k = ( b Fk + c Fk ) Y k − a Fk F k S k (33) with k = 1 , 2. The conservation laws are given by E = E + X + X , F k = F k + Y k , S = S + S + S + X + X + Y + Y , k = 1 , 2. Therefore, if total amounts are given, the system to be solved consists of 10equations in 10 variables which are the equations for the total amounts and, for instance,equations (30)-(33) for k = 1 , 2. From (32) and (33) we obtain X k = η k ES k − , Y k = δ k F k S k , k = 1 , η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). Proceeding as in the previoussystem, we find that X k = µ k Y k , µ k = c Fk /c Ek . Let ∆ k = µ k F k − E , ξ k = min( F k , E/µ k ) and Γ k = [0 , ξ k ). If E, F k > 0, then at steadystate E, F k (cid:54) = 0. From E = E − X − X and F k = F k − Y k we see that positivity of E, F k and non-negativity of Y k requires that (at least) Y k ∈ Γ k . If Y k ∈ Γ k , then X k is also anon-negative, increasing function of Y k .The situation resembles the situation of the previous system where catalysis is mediatedby two different kinases. In both systems we have X = µ Y and S = Y δ ( F − Y ) and sothe concentration of S is a function of Y . It follows that ϕ ( Y ) = S + X + Y = Y δ ( F − Y ) + µ Y + Y is a non-negative, increasing continuous function of Y ∈ Γ .For a fixed value of Y , the steady state values of the remaining variables (except X and S ) satisfy the steady state equations of a one-site phosphorylation cycle with species S , S , X , Y , E, F and total amounts E − X = E − µ Y , F and S − ϕ ( Y ). UsingResult 1 with ϕ ( · , Y ) denoting ϕ for a fixed Y , the BMSSs satisfy the relation S = ϕ ( Y , Y ) + ϕ ( Y )for any Y ∈ Γ and 0 ≤ Y < min( F , ( E − µ Y ) /µ ). Under these conditions, allconcentrations of the other chemical species are non-negative functions of Y , Y . Notethat the total amount of enzyme E − µ Y is part of the function ϕ and hence ϕ dependson Y . Indeed, we have S = µ Y η ( E − µ Y − µ Y ) . The function ϕ is increasing in Y and in Y .The equation X = η ES combined with S expressed as a function of Y provide(after isolation of X ) the following relation at steady state Y = f ( Y ) = X /µ = η ( E − µ Y ) Y µ δ ( F − Y ) + µ η Y . Y Y (A) Y Y (B) f (cid:72) F (cid:76) Y Y (C) Figure 7: Different behaviors of the function f ( Y ) = Y . The region Γ is marked ingreen. (A) corresponds to Proposition 1 (i) with Y ∈ Γ = [0 , F ). (B) and (C) differ inwhether (B) Y ∈ Γ = [0 , E/µ ) or (C) Y ∈ Γ = [0 , F ); in both cases f ( E/µ ) = 0. If F > f ( α ) (the top point of the curve) then Proposition 1 (ii)(c) applies. If F ≤ f ( α )then Proposition 1 (ii)(b) applies for (B). In (C), if F ≤ f ( F ) (the dashed line) thenProposition 1 (ii)(a) applies, whereas Proposition 1 (ii)(b) applies if f ( F ) < F ≤ f ( α ).This function resembles that in (29) except from the quadratic term in the numeratorwhich gives f a very different analytic form from that in (29).The function f is C and takes non-negative values for Y < ξ . Therefore, a BMSSsatisfies S = ϕ ( Y ) = ϕ ( Y , f ( Y )) + ϕ ( f ( Y )) , (34)for Y ∈ Γ , such that f ( Y ) ∈ Γ . Since both ϕ , ϕ are increasing functions, the behaviorof ϕ needs to be understood from the behavior of f .If we let Γ = { Y ∈ Γ | f ( Y ) ∈ Γ } , then the concentrations of all the chemical species are non-negative when expressed interm of Y , if and only if Y ∈ Γ. If Y = 0, then f (0) = 0 and it follows that ϕ (0) = 0;in particular, it follows that 0 ∈ Γ. Note that ϕ is C in Γ; however Γ might not be aconnected interval as we will see below.The function f is C and non-negative for Y ∈ Γ . Depending on whether it is monotoneor not, different forms of ϕ are expected. These behaviors can be found by computing thederivative of f (see Figure 7). LetΛ = (1 + η /δ ) µ F − E. Observe that if Λ ≤ ≤ Proposition 1. The following statements hold:(i) If Λ ≤ , then f is an increasing C function of Y ∈ Γ = [0 , F ) . Further, we have Γ = [0 , min( F , f − ( ξ )) (with f − ( ξ ) set to + ∞ if not defined). (ii) If Λ > , then there is α ∈ Γ = [0 , ξ ) such that f (cid:48) ( α ) = 0 , f is increasing in [0 , α ) and decreasing in ( α, ξ ) . In this case, we obtain:(a) If E − µ F − µ F ≥ , then there is α ≤ α such that f is increasing and C in Γ = [0 , α ) and f ( α ) = F .(b) If E − µ F − µ F < and F ≤ f ( α ) , then Γ = [0 , α ) ∪ ( α , ξ ) with α ≤ α ≤ α and f ( α ) = f ( α ) = F . Hence, f is increasing in [0 , α ) and decreasing in ( α , ξ ) .(c) If F > f ( α ) , then Γ = [0 , ξ ) . All possibilities are covered in the proposition. Indeed, the condition in (ii)(c) implies E − µ F − µ F < 0: If not, according to (ii)(a) there would be α ≤ α such that F = f ( α ) ≤ f ( α ) < F , which is a contradiction.In the cases (i) and (ii)(a), f is an increasing C function in a connected intervalΓ = [0 , ξ ). Hence, by composition of functions, ϕ in (34) is also an increasing C functionof Y ∈ Γ with ϕ (0) = 0. Additionally, in both cases ϕ tends to infinity as Y tends to ξ since either Y = f ( Y ) tends to F or Y to F . We conclude that there is exactly oneBMSS in each of these cases.The cases (ii)(b) and (ii)(c) require further analysis. In case (ii)(b), letΓ (cid:48) = [0 , α ) , Γ (cid:48)(cid:48) = ( α , ξ ) . The values α , α correspond to the values of Y for which f ( Y ) = F and hence they arethe zeros of the denominator of S : S = Y δ ( F − Y ) = Y δ ( F − f ( Y )) . In Γ (cid:48) , f is increasing and therefore ϕ is also increasing. It tends to infinity as Y tends to α . Hence, for any value of S , there is a BMSS corresponding to a Y ∈ Γ (cid:48) .In Γ (cid:48)(cid:48) , f is a decreasing function and hence it is uncertain when/whether ϕ is increasing.However, we find that (A) when Y tends to α from the right, the function ϕ ◦ f tends to+ ∞ , while ϕ = ϕ ( · , f ( · )) is bounded; (B) when Y tends to ξ from the left, the function ϕ ◦ f is bounded, while ϕ tends to + ∞ (in fact, either S or S expressed in terms of Y does). It follows that in the interval Γ (cid:48)(cid:48) = ( α , ξ ), the function ϕ starts decreasing frominfinity and ends increasing towards infinity. By continuity, there exists a minimum S min of ϕ in this interval. We see that for S > S min , at least two values of Y satisfy ϕ ( Y ) = S ;hence at least two BMSSs exist in this interval. All together, we conclude that at leastthree BMSSs occur in this case, one in Γ (cid:48) and two in Γ (cid:48)(cid:48) . Note that when S = S min thenthere are at least two (not at least three) as the two in Γ (cid:48)(cid:48) coincide.Finally, let us consider case (ii)(c), that is the case Λ > F > f ( α ) and consequentlyΓ = [0 , ξ ). The function ϕ is increasing at Y = 0 and tends to infinity as Y approaches For fixed values of F and E with Λ > 0, multistationarity occurs for certain values of F and S . The figure shows the function f for ∆ < ϕ asfunction of F .If F ≥ M (blue region) there is only one BMSS.If M > F > f ( α ) (purple region) there is mul-tistationarity, (ii)(c) in Proposition 1. If f ( α ) ≥ F ≥ − ∆ /µ (green region), then multistation-arity also occurs, (ii)b. Finally for − ∆ /µ > F (orange region) there is only one BMSS, (ii)(a). If∆ ≥ 0, then the orange region cannot occur. Y Y Figure 8: Different shapes of ϕ depending on E, F , F . ξ . Multistationarity can only occur in the form of Figure 6(c) and there will be differentvalues of Y corresponding to the same value of S , only if the function ϕ is not alwaysincreasing; that is if the derivative of ϕ has more than two zeros. Equivalently, if thereexists Y for which ϕ (cid:48) ( Y ) < 0, then we are guaranteed multistationarity. This result isstated in the following proposition: Proposition 2. Assume that Λ > . Then there exists a value M > f ( α ) such that, for F ≥ M , the derivative of ϕ is positive (except potentially for a finite number of pointswhere it is zero). Hence, ϕ is increasing. Further for M > F > f ( α ) , there exist valuesof Y for which ϕ (cid:48) ( Y ) < . We conclude that for all F ≥ M , there is exactly one BMSS for any value of S , whilefor all M > F > f ( α ) there is multistationarity for certain values of S . In particular,multistationarity occurs for S satisfying S min ≤ S ≤ S max , where S min is the smallest localminimum (excluding 0) and S max the largest local maximum of ϕ (excluding infinity), cf.Figure 6(c).Based on numerical examples we have made the following observation: When multi-stationarity occurs, the function ϕ has one local minima and one local maxima, β ≤ β ,resulting in at most three BMSSs. When F increases, β increases, while β decreases.For some F , β = β and the decreasing part of ϕ is lost and ϕ becomes increasing. Infact, the value of F for which β = β is M .All together, this implies that for fixed values of E, F , the value of F determineswhether values of S for which multistationarity occurs exist. Figure 8 shows how thenumber of steady states changes with F .We conclude that: (Motif (f)) . Let a two-site phosphorylation cycle be given with one kinase cata-lyzing phosphorylation at both sites and two different phosphatases catalyzing dephospho-rylation. Further, assume that positive total amounts E, F , F , S are given.Then, the BMSSs satisfy S = ϕ ( Y ) for Y ∈ Γ = [0 , ξ ) , where ϕ is a C function with ϕ (0) = 0 .Let Λ = (1 + η /δ ) µ F − E . Then we have:(i) If Λ ≤ or E − µ F − µ F ≥ then there is a unique BMSS.(ii) If Λ > then there exists values of F and S such that the system has more than oneBMSS. Further, there is an upper bound to F for which multistationarity can occur.This bound is independent of S . Motif (g). Next, we consider a two-site phosphorylation cycle, where phosphorylation iscatalyzed by the same kinase at both sites and dephosphorylation is catalyzed by the samephosphatase at both sites. Again, we assume sequential phosphorylation. The chemicalreactions of the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + ES + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F The differential equations describing the system are the following: ˙ S = b E X + c F Y − a E ES ˙ E = ( b E + c E ) X + ( b E + c E ) X − ( a E S + a E S ) E ˙ F = ( b F + c F ) Y + ( b F + c F ) Y − ( a F S + a F S ) F ˙ S = c E X + b E X + b F Y + c F Y − ( a F F + a E E ) S (35) ˙ S = c E X + b F Y − a F F S (36)˙ X k = − ( b Ek + c Ek ) X k + a Ek ES k − (37)˙ Y k = − ( b Fk + c Fk ) Y k + a Fk F S k (38) with k = 1 , 2. As in the previous system, the conservation laws are given by E = E + X + X , F = F + Y + Y , S = S + S + S + X + X + Y + Y . Therefore, if total amounts are given, the system to be solved consists of 9 equations in9 variables which are the equations for the total amounts and, for instance, equations(35)-(38). From (37) and (38) we obtain X k = η k ES k − , Y k = δ k F S k , k = 1 , η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). From (36) and (38) with k = 2,and from (35), (36), (38) with k = 1 , k = 2, we obtain X k = µ k Y k , µ k = c Fk /c Ek . If E, F > E, F (cid:54) = 0 at steady state. We isolate E, F from the correspondingconservation laws and write them as functions of Y , Y . Since the denominators are non-zero, we find S = µ Y η ( E − µ Y − µ Y ) , S = Y δ ( F − Y − Y ) . From Y = δ F S = δ ( F − Y − Y ) S we obtain Y = δ S ( F − Y )1+ δ S . Now from the equality X = η ES = η ( E − µ Y − µ Y ) S , we find that Y = η S ( E + δ ( E − µ F ) S ) µ + µ ( δ + η ) S + δ η ( µ − µ ) S . We can therefore write all concentrations as functions of S . Let p ( S ) = E + δ ( E − µ F ) S p ( S ) = µ F + η ( µ F − E ) S p ( S ) = µ + µ ( δ + η ) S + δ η ( µ − µ ) S . Then we have Y = δ S p ( S ) p ( S ) , Y = η S p ( S ) p ( S ) , S = µ δ S p ( S ) µ η p ( S ) , S = η S p ( S ) δ p ( S ) . This leads to the following relation for the BMSSs S = ϕ ( S ) = S + S + S + X + X + Y + Y with ϕ (0) = 0 and where each term is considered a function of S , as defined above. Forthis motif, a ϕ -function in terms of an intermediate complex cannot be obtained: Therelationship between S and Y or Y is in general not invertible.We next seek to determine the values of S that lead to non-negative concentrations ofthe species; that is to determine the region Γ, where S is defined. From (39) we see thatif Y k , S , S , S are non-negative, then so are E, F, X k at steady state. The roots of p , p are ξ = Eδ ( µ F − E ) , ξ = µ Fη ( E − µ F ) , respectively. Since p , p are in the denominators of S , S , these polynomials are notallowed to vanish. Let ∆ = µ F − E and ∆ = µ F − E . Then(i) If ∆ < 0, then p ( S ) > S < ξ . If ∆ ≥ 0, then p ( S ) > S ≥ (ii) If ∆ ≤ 0, then p ( S ) > S > 0. Otherwise, if ∆ > 0, then p ( S ) > S < ξ .The polynomial p ( S ) has degree 2 and is positive whenever S is positive and µ ≥ µ .If µ < µ , then there is exactly one positive root ξ . Therefore, we have the followingsituations:A. If µ ≥ µ (∆ ≥ ∆ ), we require p ( S ) , p ( S ) > 0. Three different scenarios occur:1) ∆ > > 0, i.e. E < µ F ≤ µ F , Γ = [0 , ξ ).2) ∆ ≤ ≥ 0, i.e. µ F ≤ E ≤ µ F , Γ = [0 , + ∞ ).3) ∆ < < 0, i.e. µ F ≤ µ F < E , Γ = [0 , ξ ).B. If µ > µ (∆ > ∆ ), we have1) ∆ > ≥ 0, i.e. E ≤ µ F < µ F , Γ = [0 , (cid:101) ξ ) for (cid:101) ξ = min( ξ, ξ ).2) ∆ ≥ < 0, i.e. µ F < E ≤ µ F , Γ = [0 , (cid:101) ξ ) ∪ ( ξ, + ∞ ) for (cid:101) ξ = min( ξ , ξ , ξ )and ξ = max( ξ , ξ , ξ ).3) ∆ < < 0, i.e. µ F < µ F < E , Γ = [0 , (cid:101) ξ ) for (cid:101) ξ = min( ξ, ξ ).Observe that in all cases, the function ϕ is C and non-negative in Γ. In the cases A.1,A.3, B.1 and B.3, when S approaches the upper limit of Γ, then ϕ tends to infinity, sinceat least one of Y , Y , S , S does. In case A.2 the function tends to infinity as S tends toinfinity. In all these cases, multistationarity would appear in the form of Figure 6(c), thatis, ϕ should decrease for some values of S .In case B.2, Γ is not connected so we let Γ = Γ (cid:48) ∪ Γ (cid:48)(cid:48) . When S tends to (cid:101) ξ , then ϕ tends to infinity. It follows that for every S , there is one BMSS located in Γ (cid:48) . Additionally,when S approaches ξ from the right, ϕ also tends to infinity, implying that ϕ comes downfrom infinity in Γ (cid:48)(cid:48) . When S tends to infinity, ϕ tends to infinity. This implies that incase B.2 the function resembles that in Figure 6(b), potentially with more increasing anddecreasing parts. In any case, when S is large, multistationarity occurs and there are atleast three steady states.Let us consider the derivatives with respect to S of the following summands in ϕ : ∂ ( Y + Y ) ∂S = µ ( η E + δ µ F ) + 2 δ η µ ( µ − µ ) F S + η δ ( µ − µ )∆ S p ( S ) ∂ ( µ Y + µ Y ) ∂S = µ ( η E + δ µ F ) + 2 δ η µ ( µ − µ ) ES − η δ µ ( µ − µ )∆ S p ( S ) ∂S ∂S = δ µ µ E F + 2 η E ∆ S − µ η δ ∆ ∆ S ( η µ p ( S )) ∂S ∂S = η µ E F − η δ µ F ∆ S − η δ ∆ ∆ S ( δ p ( S )) In case A.2 these derivatives are all positive since ∆ ≤ 0, ∆ ≥ µ − µ ≥ ϕ is decreasing for some value of S .Thus, further analysis of the derivatives or the function ϕ itself is required. For example, incase A.3, the following choices of numerical values provide multiple steady states: F = 3, E = 10, δ = 10, δ = 100, µ = 1, µ = 3, η = 0 . η = 100.We conclude that: Result 7 (Motif (g)) . Let a two-site phosphorylation cycle be given with phosphorylationat both sites being catalyzed by the same kinase and dephosphorylation at both sites beingcatalyzed by the same phosphatase. Further, assume that the total amounts E, F , S arepositive. Then, any BMSSs satisfy S = ϕ ( S ) for S ∈ Γ (with Γ defined above) where ϕ is a C function with ϕ (0) = 0 .Additionally,(i) If µ F ≤ E ≤ µ F , then for any total amount S there is a unique BMSS.(ii) If µ F < E < µ F , then there exists a value S min such that for any S > S min thesystem has at least three BMSSs. B.3 Modification of two different substrates Motif (h). In this system, two cycles are connected through a joint catalyzing kinase.The chemical reactions of the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E P + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) P + ES + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F P + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) P + F The differential equations describing the system are the following: ˙ S = b E X + c F Y − a E ES ˙ S = b F Y + c E X − a F F S (40)˙ X = − ( b E + c E ) X + a E ES (41)˙ Y = − ( b F + c F ) Y + a F F S (42)˙ F = ( b F + c F ) Y − a F F S ˙ E = ( b E + c E ) X + ( b E + c E ) X − ( a E S + a E P ) E. ˙ P = b E X + c F Y − a E EP ˙ P = b F Y + c E X − a F F P (43)˙ X = − ( b E + c E ) X + a E EP (44)˙ Y = − ( b F + c F ) Y + a F F P (45)˙ F = ( b F + c F ) Y − a F F P The conservation laws are given by E = E + X + X , F k = F k + Y k , S = S + S + X + Y , P = P + P + X + Y , with k = 1 , 2. Therefore, if total amounts are given, the system to be solved consists of 11equations in 11 variables which are the equations for the total amounts (5 equations) and,for instance, equations (40)-(45). From (41), (42), (44) and (45) we obtain X = η ES , X = η EP , Y = δ F S , Y = δ F P , with constants η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). From these equations, (40) and(43) we obtain X k = µ k Y k , µ k = c Fk /c Ek . If E, F k > 0, then at steady state E, F k (cid:54) = 0. Since F k = F − Y k we require 0 ≤ Y k < F k for any BMSS. Using the conservation laws for P and F we have P = µ Y η E + Y δ ( F − Y ) + µ Y + Y , and it follows that E is an increasing C function of Y E = g ( Y ) = µ δ Y ( F − Y ) η ( δ P ( F − Y ) − Y − δ ( µ + 1) Y ( F − Y ))The numerator is non-negative for 0 ≤ Y < F . The denominator is a degree two polyno-mial in Y with positive independent and leading coefficients. For Y = F , the denominatoris negative and thus has two positive real roots. It follows that for E to be non-negativewe require Y ∈ [0 , ξ ) where ξ < F is the first positive root of the denominator. We havethat g (0) = 0 and that g ( Y ) goes to infinity as Y tends to ξ . To sum up, for 0 ≤ Y < ξ ,the steady state values of F , P , P , X , E are non-negative as well.Using the conservation law for E , it follows that E = E + X + X = g ( Y ) + µ Y + µ Y , (46)and Y is a decreasing C function of Y , h ( Y ) = Y , defined on [0 , E/µ ] such that h ( E/µ ) = 0 and h (0) < ξ . Consequently, g ( h ( Y )) = E is a decreasing C functiondefined on [0 , E/µ ] and since E > 0, we require Y ∈ [0 , E/µ ).Let Γ = [0 , ξ ) with ξ = min( F , E/µ ) and consider the conservation law for S . Weobtain S = ϕ ( Y ) = µ Y η g ( h ( Y )) + Y δ ( F − Y ) + µ Y + Y . Then, ϕ is an increasing C function defined on Γ with image R + and all steady stateconcentrations are non-negative if and only if Y ∈ Γ. Therefore, we have shown: (Motif (h)) . Let two one-site modification cycles with joint kinase and distinctphosphatases be given. Further, assume that the total amounts S, P , E, F k , k = 1 , , arepositive. Then, the system has a unique BMSS.Specifically, the BMSS satisfies S = ϕ ( Y ) for Y ∈ Γ = [0 , ξ ) , where S = ϕ ( Y ) = µ Y η g ( h ( Y )) + Y δ ( F − Y ) + µ Y + Y is an increasing C function which tends to infinity as Y tends to ξ and fulfills ϕ (0) = 0 . Remark 3. In equation (46), if we isolate Y instead of Y , we would obtain a relationshipbetween Y and Y , Y = h ( Y ) and the steady state relation in terms of Y , S = ψ ( Y ).In this case, Y ∈ ( ξ, ξ (cid:48) ) where ξ = h ( ξ ) and ξ (cid:48) is the pre-image of E of the increasingfunction µ Y + g ( Y ). If ξ = E/µ , then ξ = 0. Motif (i). In this system, two cycles are connected through a joint catalyzing kinase anda joint catalyzing phosphatase. The chemical reactions of the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E P + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) P + ES + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F P + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) P + F The differential equations describing the system are the following: ˙ S = b E X + c F Y − a E ES ˙ S = b F Y + c E X − a F F S (47)˙ X = − ( b E + c E ) X + a E ES (48)˙ Y = − ( b F + c F ) Y + a F F S (49)˙ F = ( b F + c F ) Y + ( b F + c F ) Y − ( a F S + a F P ) F ˙ E = ( b E + c E ) X + ( b E + c E ) X − ( a E S + a E P ) E. ˙ P = b E X + c F Y − a E EP ˙ P = b F Y + c E X − a F F P (50)˙ X = − ( b E + c E ) X + a E EP (51)˙ Y = − ( b F + c F ) Y + a F F P (52) The conservation laws are given by E = E + X + X , F = F + Y + Y , S = S + S + X + Y , P = P + P + X + Y . If total amounts are given, the system to be solved consists of 10 equations in 10 variableswhich are chosen to be the equations for the total amounts (4 equations) and equations(47)-(52). As usual we derive X = η ES , X = η EP , Y = δ F S , Y = δ F P , X k = µ k Y k , k = 1 , 2, with constants η k = a Ek / ( b Ek + c Ek ), δ k = a Fk / ( b Fk + c Fk ) and µ k = c Fk /c Ek .Let ξ k = min( F , E/µ k ), k = 1 , 2. For fixed Y ∈ [0 , ξ ), the equations above providethe steady state equations corresponding to a one-site cycle with species E, F, S , S , X , Y and total amounts S, E − µ Y , F − Y . Analogously, for fixed Y ∈ [0 , ξ ), we obtain theequations for a one-site cycle with species E, F, P , P , X , Y and total amounts P , E − µ Y , F − Y . Therefore, using Result 1 we obtain that the steady states are solutions tothe system S = ϕ ( Y , Y ) = µ Y η ( E − µ Y − µ Y ) + Y δ ( F − Y − Y ) + µ Y + Y P = ϕ ( Y , Y ) = µ Y η ( E − µ Y − µ Y ) + Y δ ( F − Y − Y ) + µ Y + Y for any Y , Y and for the solution to be a BMSS we require that Y + Y < F and µ Y + µ Y < E . All species concentrations derived with Y , Y fulfilling these conditionsare non-negative.Note that ϕ , ϕ are increasing functions of both Y and Y . Therefore, using ϕ for afixed S , there is a decreasing C function f ( Y ) = Y , defined for Y ∈ [0 , ξ ). Further, since f ( Y ) is the steady state value of Y in the first cycle with total amounts S, E − µ Y , F − Y ,the derived concentrations E, F, S , S , X are also non-negative. Thus, so are P , P , X for Y ∈ [0 , ξ ).Let Γ = [0 , ξ ) so that all concentrations are non-negative if and only if Y ∈ Γ. Weconclude that the steady states of the system are described by a C function P = ϕ ( Y ) = ϕ ( f ( Y ) , Y )defined for Y ∈ Γ.The behavior of the function ϕ determines the presence/absence of multiple steadystates. Note that when Y tends to ξ , then Y = f ( Y ) tends to zero and ϕ tends toinfinity, implying that ϕ increases as we approach the upper limit of Γ. Additionally, at Y = 0, f ( Y ) is finite, ϕ (0) = 0 and hence the function is increasing at 0 as well (itis positive for Y > ϕ (cid:48) ( Y ) < Y ∈ Γ.To proceed, let ψ = ( ϕ , ϕ ) : R → R and let J ψ be the Jacobian matrix of ψ , thatis, J ψ = (cid:32) ∂ϕ ∂Y ∂ϕ ∂Y ∂ϕ ∂Y ∂ϕ ∂Y (cid:33) . Then, a simple observation (proved in Proposition 7) shows that ϕ (cid:48) ( Y ) < J ψ ( f ( Y ) , Y )) < 0. Therefore, if det( J ψ ( Y , Y )) < E, F , Y , Y , thenfor total amounts S = ϕ ( Y , Y ) and P = ϕ ( Y , Y ), the system exhibits multistationarity. Let σ = ( µ − µ )( µ δ η − µ δ η ) and Γ = { ( Y , Y ) ∈ R | Y + Y Result 9 (Motif (i)) . Let two one-site modification cycles with joint kinase and jointphosphatase be given. Further, assume that the total amounts S, P , E, F are positive. Then,the BMSSs satisfy P = ϕ ( Y ) for Y ∈ Γ , where ϕ is a C function which tends to infinityas Y tends to ξ and fulfills ϕ (0) = 0 .Let σ = ( µ − µ )( µ δ η − µ δ η ) . Then: • The function ϕ is always increasing if either ( i ) σ ≥ or ( ii ) σ < and either (a) µ − µ > , together with µ F ≥ E or E ≥ µ F , or (b) µ − µ > , together with µ F ≥ E or E ≥ µ F . • If σ < and either µ F > E > µ F or µ F > E > µ F , then, there exist values Y ∈ Γ for which ϕ (cid:48) ( Y ) < . Hence multistationarity occurs. In this case the totalamounts S, P are required to be large. B.4 Cascade motifs Motif (j). We consider here the combination of two one-site modification cycles in acascade motif with a specific phosphatase acting in each layer. The chemical reactions ofthe system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E P + S a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) P + S S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F P + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) P + F The differential equations describing the system are the following: ˙ S = b E X + c F Y − a E ES ˙ P = b E X + c F Y − a E S P ˙ E = ( b E + c E ) X − a E ES ˙ F = ( b F + c F ) Y − a F F S ˙ F = ( b F + c F ) Y − a F F P ˙ X = a E ES − ( b E + c E ) X (53)˙ X = a E S P − ( b E + c E ) X (54)˙ Y = a F F S − ( b F + c F ) Y (55)˙ Y = a F F P − ( b F + c F ) Y (56)˙ P = b F Y + c E X − a F F P (57)˙ S = c E X + b F Y + ( b E + c E ) X − ( a F F + a E P ) S . (58) The conservation laws are given by E = E + X , F k = F k + Y k , S = S + S + X + X + Y , P = P + P + X + Y for k = 1 , 2. If total amounts are given, the system to be solved consists of 11 equationsin 11 variables which are the equations for the total amounts and, for instance, equations(53)-(58). As usual we obtain X = η ES , X = η S P , Y = δ F S , Y = δ F P with constants η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). Equations (56),(57) ( k = 2)and (54),(55),(58), ( k = 1) provide the relations X k = µ k Y k , µ k = c Fk /c Ek for k = 1 , E, F , F , X , X are solved in terms of Y , Y from the total amounts E, F , F and the two relations above. If E, F k > 0, then E, F k (cid:54) = 0 at steady state. Hence,any BMSS satisfies 0 ≤ Y k < F k , 0 ≤ Y < E/µ .Let ξ = min( F , E/µ ). Then we find S = µ Y η ( E − µ Y ) and S = Y δ ( F − Y ) , which are non-negative, increasing, continuous functions of Y ∈ [0 , ξ ). It follows that ϕ ( Y ) = S + S + X + Y is a non-negative, increasing, continuous function of Y ∈ [0 , ξ ). Also, it tends to infinityas Y tends to ξ and thus the image of ϕ over [0 , ξ ) is R + . From the conservation lawfor S we find Y = f ( Y ) = 1 µ ( S − ϕ ( Y )) , which is a decreasing function in Y ∈ [0 , ξ ). The concentration Y is non-negative provided ϕ ( Y ) ≤ S . Let ξ = ϕ − ( S ) ( ϕ is increasing, hence invertible). Then Y , S , S , X areall functions of Y and non-negative provided Y ∈ [0 , ξ ]. We have that f ( ξ ) = 0 and f (0) = S/µ . It follows from the inverse function theorem that f can be inverted so thatthere exists a continuous decreasing function Y = g ( Y ) for Y ∈ [0 , S/µ ] with g (0) = ξ and g ( S/µ ) = 0. Then S , S , X , Y are all functions of Y and non-negative provided Y ∈ [0 , S/µ ]. Note that if S > S = 0 is not a solution of the steady stateequations and hence any BMSS satisfies S > 0. Since S = g ( Y ) / ( δ ( F − g ( Y ))) werequire Y < S/µ .Let ξ = min( F , S/µ ) and Γ = [0 , ξ ). Next we find the relations P = µ Y η S , and , P = Y δ ( F − Y ) , which are non-negative continuous functions of Y ∈ Γ. Since S is increasing in Y andthus decreasing in Y , we have that both P , P are increasing in Y .We have seen that all concentrations at steady state are non-negative if and only if Y ∈ Γ. Further, when Y tends to ξ , either P or P tend to infinity. Using the conservationlaw for P we obtain: Result 10 (Motif (j)) . Let a cascade of one-site modification cycles be given with positivetotal amounts S, P , E, F , F . Then the system has a unique BMSS. Specifically, the BMSSsatisfies P = ϕ ( Y ) for Y in Γ = [0 , ξ ) , where P = ϕ ( Y ) = δ µ Y ( F − g ( Y )) g ( Y ) + Y δ ( F − Y ) + (1 + µ ) Y is an increasing C function, which tends to infinity as Y tends to ξ and fulfills ϕ (0) = 0 . Motif (k). We consider here the combination of two one-site modification cycles in acascade motif where the phosphatase is not layer specific; that is the same phosphataseacts in both layers. The chemical reactions of the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E P + S a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) P + S S + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F P + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) P + F The differential equations describing the system are the following: ˙ S = b E X + c F Y − a E ES ˙ P = b E X + c F Y − a E S P ˙ E = ( b E + c E ) X − a E ES ˙ F = ( b F + c F ) Y + ( b F + c F ) Y − ( a F S + a F P ) F ˙ S = c E X + b F Y + ( b E + c E ) X − ( a F F + a E P ) S (59) ˙ X = a E ES − ( b E + c E ) X (60)˙ X = a E S P − ( b E + c E ) X (61)˙ Y = a F F S − ( b F + c F ) Y (62)˙ Y = a F F P − ( b F + c F ) Y (63)˙ P = b F Y + c E X − a F F P . (64) The conservation laws are given by E = E + X , F = F + Y + Y , S = S + S + X + X + Y , P = P + P + X + Y . If total amounts are given, the system to be solved consists of 10 equations in 10 variableswhich are the equations for the total amounts and, for instance, equations (59)-(64). Fromthe latter equations, we obtain X = η ES , X = η S P , Y = δ F S , Y = δ F P with constants η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). Equations (63),(64) ( k = 2)and (59),(61) and (62) ( k = 1) provide X k = µ k Y k , µ k = c Fk /c Ek for k = 1 , E, F > 0, then E, F (cid:54) = 0 at steady state. Therefore, any BMSS satisfies0 ≤ Y k < F and Y < E/µ . If S = 0, then it follows from the steady state equationsthat S = X = X = Y = 0 and thus S = 0. Hence, if S > S > Y = X /µ < S/µ for any BMSS.Let ξ = min( S/µ , F ), ξ = min( E/µ , F ) and Γ k = [0 , ξ k ), k = 1 , 2. For fixed Y ∈ Γ , let ξ ( Y ) = min( E/µ , F − Y ) and note that ξ ( Y ) ≤ ξ . For Y fixed, thesteady states satisfy the steady state equations of a one-site phosphorylation cycle withspecies S , S , X , Y , E, F and positive total amounts E , F − Y and S − µ Y (these areindependent of P ). Using Result 1 with ϕ ( Y , Y ) denoting ϕ ( Y ) for the fixed Y , theBMSSs satisfy the relation S = Φ( Y , Y ) = ϕ ( Y , Y ) + µ Y for 0 ≤ Y < ξ ( Y ) ≤ ξ . Under these conditions, the concentrations S , S , E, X , Y , X , Y , F are non-negative. Note that S = µ Y η ( E − µ Y ) , S = Y δ ( F − Y − Y ) , (with non-zero denominators) and ϕ ( Y , Y ) = S + S + X + Y , such that only S depends on Y . Since the BMSS is unique in the one-site cycle, it follows that for every Y ∈ Γ there exists Y ∈ Γ ( Y is non-zero provided S > 0) satisfying S = Φ( Y , Y ).Thus, there is a function f so that for any BMSS, Y = f ( Y ) . (65)The function f is C and decreasing for Y ∈ Γ . Indeed, this follows from the implicitfunction theorem and the fact that Φ is C in Γ , the derivative with respect to Y ispositive and the derivative with respect to Y is positive too. This can be checked by directcomputation.By construction, we have 0 < f ( Y ) < ξ ( Y ) for every Y . Using Remark 1, the value Y = f ( Y ) is the first positive root of the polynomial in Y (for fixed Y ) obtained fromΦ( Y , Y ) by elimination of denominators. Also, if Y = 0 then Y = f (0) is the BMSS valueof a one-site phosphorylation system with total amounts S, E, F and hence it is positive.Consider now the conservation law for P : P = µ Y η S = µ δ ( F − Y − Y ) Y η Y , P = Y δ F = Y δ ( F − Y − Y ) , and hence P = ϕ ( Y , Y ) is a function of Y and Y , which is C with respect to eachvariable. Further, each summand P , P , X and Y of P is non-negative if Y ∈ Γ and soall concentrations at steady state are non-negative if and only if Y ∈ Γ . Accordingly, welet Γ = Γ .The BMSSs are characterized by the relation P = ϕ ( Y ) = ϕ ( f ( Y ) , Y ) . The function ϕ is C for Y ∈ Γ. Further, ϕ (0) = 0 and as Y tends ξ , f ( Y ) tends to zeroand hence ϕ ( Y ) tends to + ∞ . Consequently, there is at least one BMSS.Determination of the number of BMSSs follows from the behavior of the function ϕ inΓ. Since Γ is a connected interval, multistationarity can only occur if ϕ (cid:48) ( Y ) < Y ∈ Γ (Figure 6(c)). In this case multiple BMSSs occur for P min ≤ P ≤ P max forsome P min < P max . Analysis of the behavior of this function is not straightforward butstill some general conclusions can be derived.Clearly Y + µ Y is an increasing function of Y . The derivatives of P , P with respectto Y are ∂P ∂Y = µ δ η f ( f ( F − f − Y ) − f (cid:48) Y ( F − Y )) and ∂P ∂Y = F − f + f (cid:48) Y δ ( F − f − Y ) . Either ∂P ∂Y or ∂P ∂Y must be negative for ϕ (cid:48) ( Y ) < = µ F − E . Let σ = min( S, P ) and σ µ = min(1 + µ , µ / . If ∆ > and σ µ ( F − E/µ ) > σ then ∂P ∂Y ( Y ) , ∂P ∂Y ( Y ) > for all Y ∈ Γ . Consequently, multistationarity cannot occur. Finally, let us analyze the situation where S is large. Proposition 5. The following statements hold:(i) Fix Y > F − E/µ . Then, as S tends to + ∞ , ϕ (cid:48) ( Y ) becomes positive, ϕ (cid:48) ( Y ) > .(ii) Assume that ∆ > . Let ∆ = 27 µ µ δ η E ( µ F − E ) − δ ( δ µ µ F − ( µ δ + η (1 + µ )) E ) . If ∆ < , then ϕ decreases for some Y ≤ F − E/µ as S becomes large. Implicitly Proposition 5 assumes that P is large too: In (i), Y is fixed and S becomeslarge restricting the possible values of Y and P . In (ii), the same is in play. Thus,contradictory conclusions cannot be reached from the two propositions. For example, if F (cid:29) E , then ∆ < S (and thus P ) large,multistationarity exists. If F (cid:29) E and P , S fixed, Proposition 4 ensures monostationarity.Also, it follows from Proposition 5 (i) that if F − E/µ < 0, multistationarity for S largecannot occur.Multistationarity can also occur if ∆ ≤ 0. For instance, the total amounts F = 4 , E =10 , S = 50 , P = 195 and rate constants µ = 15 , µ = 2 , η = η = δ = 1 , δ = 0 . Y = 2 . , . , . 26, respectively. However, contrary tothe situation with the reversed inequality, ∆ > 0, multistationarity cannot occur for S large (Proposition 5(i)).We have the following result: Result 11 (Motif (k)) . Let a cascade be given with one phosphatase acting in both layers.Further, assume that the total amounts E, F , S, P are positive.(i) If F (cid:29) E and S large, then there exist values P min < P max such that for all P min ≤ P ≤ P max the system admits more than one BMSS.(ii) If F − E/µ > σ/σ µ with σ = min( S, P ) and σ µ = min(1 + µ , µ / , then multista-tionarity cannot occur.(iii) If F − E/µ < , then multistationarity cannot occur for large S . We conclude that multistationarity occurs in this motif. Note that statement (iii)implies that if multistationarity occurs for some values E/µ > F and some S, P , thenincreasing S eventually makes the system monostationary. Observe also that according to(i) and (ii) together, P min is required to be so large that the condition F − E/µ > P min /σ µ is not fulfilled. This motif is a combination of Motif (c) and Motif (j). The chemical reactionsof the system are: S + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) S + E P + S a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) P + S P + E a E (cid:47) (cid:47) X b E (cid:111) (cid:111) c E (cid:47) (cid:47) P + ES + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) S + F P + F a F (cid:47) (cid:47) Y b F (cid:111) (cid:111) c F (cid:47) (cid:47) P + F The differential equations describing the system are the following: ˙ P = b E X + b E X + c F Y − a E S P − a E EP ˙ E = ( b E + c E ) X + ( b E + c E ) X − ( a E S + a E P ) E ˙ F = ( b F + c F ) Y − a F F S ˙ F = ( b F + c F ) Y − a F F P ˙ P = b F Y + c E X + c E X − a F F P (66)˙ S = c E X + b F Y + ( b E + c E ) X − ( a F F + a E P ) S (67) ˙ S = b E X + c F Y − a E ES ˙ X = a E ES − ( b E + c E ) X (68)˙ X = a E S P − ( b E + c E ) X (69)˙ X = a E EP − ( b E + c E ) X (70)˙ Y = a F F S − ( b F + c F ) Y (71)˙ Y = a F F P − ( b F + c F ) Y (72) The conservation laws are given by E = E + X + X , F k = F k + Y k , S = S + S + X + X + Y , P = P + P + X + Y + X . If total amounts are given, the system to be solved consists of 12 equations in 12 variableswhich are the equations for the total amounts and, for instance, equations (66)-(72). Weobtain X = η ES , X = η S P , X = η EP , Y = δ F S , Y = δ F P with constants η k = a Ek / ( b Ek + c Ek ) and δ k = a Fk / ( b Fk + c Fk ). Equations (67),(69),(71), and(66),(72) provide X = µ Y , µ − X + µ − X − Y = 0with µ = c F /c E , µ = c F /c E and µ = c F /c E .We write X as a function of Y and isolate E, F , F using the conservation laws forthe total amounts E, F , F . If E, F k > E, F k (cid:54) = 0 at steady state and it followsthat E − X − X > ≤ Y k < F k and Y < E/µ . It follows as well that X < µ F . If in addition, S > P > S , S , X , Y , X , P , P , Y , X > ξ = min( F , E/µ ), Γ = (0 , ξ ) and Γ = (0 , µ F ). It is convenient for this motifto exclude 0 in the intervals. A necessary condition for positivity (i.e. a BMSS) is thus Y ∈ Γ , X ∈ Γ . We proceed to write S as an increasing C functions of Y ∈ Γ and P as a C function increasing in X ∈ Γ and decreasing in Y ∈ Γ : S = Y δ ( F − Y ) , P = δ X ( F − Y ) η Y . Using X = η EP , we express X as a positive C function, increasing in X ∈ Γ anddecreasing in Y ∈ Γ : X = Φ X ( Y , X ) = δ η X ( F − Y )( E − µ Y ) η Y + δ η X ( F − Y ) . We proceed to write Y as a positive C function, increasing in X ∈ Γ and decreasing in Y ∈ Γ : Y = Φ Y ( Y , X ) = µ − X + µ − Φ X ( Y , X ) . For Y < F , we require ( Y , X ) ∈ Γ (cid:48) withΓ (cid:48) = { ( Y , X ) | Y ∈ Γ , X ∈ Γ , Φ Y ( Y , X ) < F } . From P = Y / ( δ F ) we have P = Φ Y ( Y , X ) δ ( F − Φ Y ( Y , X )) , which is increasing in X , decreasing in Y and positive and continuous provided ( Y , X ) ∈ Γ (cid:48) . Finally, S = µ Y η ( E − µ Y − Φ X ( Y , X ))which is positive and C in Γ (cid:48) .To sum up: All concentrations at steady state are positive if only if ( Y , X ) ∈ Γ (cid:48) . Tofind a final relation between X , Y , we consider the total amount P . For Y , X in Γ (cid:48) , P = ϕ P ( Y , X ) , where ϕ P a positive C function. Note also that for every fixed value of Y ∈ Γ , ϕ P ( Y , · )increases in X , is well-defined at X = 0 where it is zero and tends to infinity as Φ Y ( Y , X )tends to F . Therefore, for P > Y ∈ Γ , there exists 0 < X satisfying P = ϕ P ( Y , X ) and Φ Y ( Y , X ) < F . Thus, there exists a function f defined on Γ forwhich X = f ( Y )at steady state and ( Y , f ( Y )) ∈ Γ (cid:48) . Since ϕ P is increasing in X and decreasing in Y ,the function f is C and increasing in Γ := Γ . All concentrations at steady state are positive if and only if Y ∈ Γ. The concentrations S , X , Y are increasing functions of Y and independent of X . X is increasing in Y (after substitution of f ). S might be increasing or decreasing in Y since Φ X ( Y , X ) isnot increasing in Y . In any case, if we insert these values into S , we obtain that the BMSSsare given by a relation S = ϕ ( Y ) , where ϕ is a positive C function defined on Γ. When Y tends to ξ , the function tendsto infinity, since f ( Y ) is bounded by µ F . When Y tends to zero, then X tends to zeroas well as can be seen from the equations Y = δ F S and X = η S P and the fact that P is bounded by P . Therefore, there is at least one BMSS. Multistationarity can onlyoccur in the form of Figure 6(c), implying the existence of lower and upper bounds on S , S min < S max , for the existence of multistationarity.If S is increasing in Y , then so is ϕ and there is exactly one BMSS. The derivative of S with respect to Y is ∂S ∂Y = µ ( η E + δ η X ( µ F − E )) η η ( E − µ Y ) + ∂S ∂X f (cid:48) . Since ( ∂S /∂X ) f (cid:48) > 0, we have that if (a) µ F − E ≥ µ F − E < η E + δ η F ( µ F − E ) ≥ 0, then ∂S ∂Y > Proposition 6. Fix the total amount F . Then there exist total amounts P , S, F , E satisfying µ F − E < for which ϕ (cid:48) ( Y ) < for some Y . It follows from the discussion above the proposition that η E + δ η F ( µ F − E ) < F and E are chosen so large that this condition isfulfilled. The proof of this proposition is provided in Appendix A. Result 12 (Motif (l)) . Consider a cascade with different phosphatases acting in the twolayers and where the kinase of the first layer also acts in the second layer. Assume thatthe total amounts E, F , F , S, P are positive.(i) If (a) µ F − E ≥ or (b) µ F − E < and η E + δ η F ( µ F − E ) ≥ , then thesystem has exactly one BMSS.(ii) For any F and E, F large and satisfying µ F − E < and η E + δ η F ( µ F − E ) < , there exist values P , S for which the system displays multistationarity. Fur-ther, in this case, there exists S min < S max for which multiple steady states occur for S min ≤ S ≤ S max . C Stability analysis We provide here the mathematical details of the stability analysis. We prove equation [3]in the main text and show that for the motifs exhibiting multistationarity, the equilibriumpoints for which ϕ (cid:48) < C.1 The determinant of the Jacobian We prove here equation [3] of the main text. For any 1 ≤ j ≤ n , we define the projection π ( j ) of R n to R n − that removes the j -th coordinate by π ( j ) : R n → R n − ( x , . . . , x n ) (cid:55)→ ( x , . . . , (cid:98) x j , . . . , x n ) . For simplicity, let x ( j ) = π ( j ) ( x ) for any x ∈ R n and let Ω ( j ) = π ( j ) (Ω) ⊆ R n − for any openset Ω ⊂ R n . Note that the latter is also an open set (projection maps are open maps).For a differentiable function f = ( f , . . . , f n ) : Ω ⊆ R n → R n , we denote by J f ( x ) theJacobian of f at x ∈ Ω. The ( i, j ) entry of J f ( x ) is ∂f i /∂x j ( x ).In the next proposition, Ω is an open neighborhood of z , suitably chosen. Proposition 7. Let f = ( f , . . . , f n ) : Ω ⊆ R n → R n be a differentiable function definedon an open set Ω such that there is z ∈ Ω with f ( z ) = 0 . Assume that x j can be eliminatedfrom the equation f i = 0 ; that is, there exists a differentiable function ψ : Ω ( j ) ⊆ R n − → R such that x j = ψ ( x ( j ) ) on { x ∈ Ω | f i ( x ) = 0 } . Define ¯ f : Ω ( j ) → R n − by ¯ f k ( x ) = f k ( x , . . . , x j − , ψ ( x ) , x j , . . . , x n − ) for all k (cid:54) = i . Then, the determinant of the Jacobian of f at z satisfies ( − i + j ∂f i ∂x j ( z ) det( J ¯ f ( z ( j ) )) = det( J f ( z )) . (73) Proof. It is sufficient to show that the proposition holds for i = j = 1. For other values of i and j the result is obtained by reorganizing rows and columns and keeping track of thesign of the determinant. For i = j = 1, the proposition states that ∂f ∂x ( z ) det( J ¯ f ( z (1) )) = det( J f ( z )) . where z (1) = ( z , . . . , z n ) and z = ( z , . . . , z n ).If the Jacobian is written in column vector notation, we have that J f ( z ) = ( D f ( z ) , . . . , D n f ( z )) where D k f = (cid:16) ∂f ∂x k , . . . , ∂f n ∂x k (cid:17) T . That is, D k f is the vector of the partial derivatives of f with respect to x k . By assumption, z = ( ψ ( z (1) ) , z , . . . , z n ) ∈ R n . Observe that if ( ∂f /∂x )( z ) = 0 thenfrom ∂f ∂x ( z ) ∂ψ∂x k ( z (1) ) + ∂f ∂x k ( z ) = 0, we have that ( ∂f /∂x k )( z ) = 0 for all k . Consequently,det( J f ( z )) = 0 (the matrix has one row of zeroes) and hence the proposition holds.Hence, we can assume that ( ∂f /∂x )( z ) (cid:54) = 0. By implicit differentiation, ∂ψ∂x k ( z (1) ) = − ∂f /∂x k ∂f /∂x ( z ) . (74)Additionally, by the chain rule we have ( k (cid:54) = 1) ∂ ¯ f l ∂x k ( z (1) ) = ∂f l ∂x k ( z ) + ∂f l ∂x ( z ) ∂ψ∂x k ( z (1) ) . Let f (1) = ( f , . . . , f n ) be the projection of f onto the last n − J ¯ f ( z (1) ) = ( D ¯ f , . . . , D n ¯ f )( z (1) ) = ( D f (1) , . . . , D n f (1) )( z )+ (cid:16) D f (1) ( z ) ∂ψ∂x ( z (1) ) , . . . , D f (1) ( z ) ∂ψ∂x n ( z (1) ) (cid:17) , where matrices are written as column vectors.Note that ∂ψ∂x k ( z (1) ) is a scalar and D f (1) ( z ) ∂ψ∂x k ( z (1) ) is the vector with l -th component ∂f l ∂x ( z ) ∂ψ∂x k ( z (1) ). Further, using the multilinear expansion of a determinant, we obtaindet( J ¯ f ( z (1) )) = det(( D f (1) , . . . , D n f (1) )( z ))+ n (cid:88) k =2 det( D f (1) ( z ) , . . . , D f (1) ( z ) ∂ψ∂x k ( z (1) ) , . . . , D n f (1) ( z ))+ (cid:88) { k ,...,k l }⊂{ ,...,n } l> det( a , . . . , a n ) , where a s = D s f (1) ( z ) for s (cid:54) = k , . . . , k l and a s = D f (1) ( z ) ∂ψ∂x s ( z (1) ) otherwise. In each ofthese summands, there are at least two terms of the second form, say D f (1) ( z ) ∂ψ∂x k ( z (1) )and D f (1) ( z ) ∂ψ∂x m ( z (1) ) for some k (cid:54) = m . If the two columns are non-zero, they are linearlydependent and the determinant of each of the matrices ( a , . . . , a n ) is zero. Further, notethat det( D f (1) ( z ) , . . . , D f (1) ( z ) ∂ψ∂x k ( z (1) ) , . . . , D n f (1) ( z ))= ∂ψ∂x k ( z (1) ) det(( D f (1) , . . . , D f (1) , . . . , D n f (1) )( z ))= ( − k − ∂ψ∂x k ( z (1) ) det(( D f (1) , . . . , (cid:92) D k f (1) , . . . , D n f (1) )( z )) Finally, using (74), we obtain ∂f ∂x ( z ) det( J ¯ f ( z (1) )) = n (cid:88) k =1 ( − k − ∂f ∂x k ( z ) det( D f (1) , . . . , (cid:92) D k f (1) , . . . , D n f (1) )( z )= det( J f ( z ))by considering the development of the determinant of J f ( z ) along the first row. Application. Let a dynamical system in Ω ⊆ R n be given˙ x = f ( x )with x = ( x , . . . , x n ) and f = ( f , . . . , f n ). Let z = ( z , . . . , z n ) be an equilibrium point,i.e., f ( z , . . . , z n ) = 0. Let det( J f ( z )) be the determinant of the Jacobian of f at z . Wemake the following observation: • If n is odd and det( J f ( z )) > 0, then z is unstable. • If n is even and det( J f ( z )) < 0, then z is unstable.Indeed, if z is (asymptotically) stable and det( J f ( z )) (cid:54) = 0, all eigenvalues have negativereal parts. Since det( J f ( z )) is a real number (it is the determinant of a real matrix),the complex eigenvalues come in pairs of conjugates a + bi, a − bi and the product of theeigenvalues is a positive number. If n is odd and all eigenvalues have negative real parts,their product must be negative and hence the determinant of the Jacobian, which equalsthe product of the eigenvalues, must be negative. If n is even and z (asymptotically) stable,then the product of the eigenvalues must be positive.An equilibrium point satisfies f ( z , . . . , z n ) = 0. If elimination of variables can beapplied then we can use Proposition 7 to track the sign of the determinant and potentiallydetect instability. Chemical reaction networks. Consider any of the motifs, or in general, any chemicalreaction network with conservation laws. The conserved total amounts imply that thedynamics of the associated dynamical system takes place in a fixed subspace of R n . Ingeneral, we have a dynamical system ˙ x = f ( x )and a series of (independent) conservation laws g i ( x ) − c i = 0 , i = 1 . . . , k (75) where g i a linear function of x satisfying g i ( ˙ x ) = 0 and c i ∈ R (these correspond to thetotal amounts). By independence we mean that the rank of this linear system is k . Theconservation laws do not depend on the rate constants and g i does not depend on the totalamounts.The existence of conservation laws implies that the determinant of the Jacobian of f at any point x is zero, since the matrix has linear relations among the rows. Therefore,stability of equilibrium points cannot be analyzed directly, but need to be considered insidethe stoichiometry class they belong to.Since (75) is a linear system of rank k , Gauss elimination allows the elimination of k variables. For simplicity, we can rename the variables such that the eliminated variablesare x , . . . , x k and those that remain are x k +1 , ..., x n . Apply the same renaming to thefunctions f i , such that if x l is now variable x j , function f l is labeled f j . By elimination,there exist (polynomial) functions x i = ψ i ( x k +1 , . . . , x n , c ) , i = 1 , . . . , k such that f i ( ψ , . . . , ψ k , x k +1 , . . . , x n ) = 0. Here c = ( c , . . . , c k ) is the vector of initial totalamounts.For a fixed stoichiometry class c = ( c , . . . , c k ) and ¯ x = ( x k +1 , . . . , x n ), let¯ f j (¯ x, c ) = f j ( ψ (¯ x, c ) , . . . , ψ k (¯ x, c ) , x k +1 , . . . , x n ) , j = k + 1 , . . . , n. To investigate stability of an equilibrium point ( z , . . . , z n ) belonging to the stoichiometryclass c with z i − ψ i ( z k +1 , . . . , z n , c ) = 0, we consider the the eigenvalues of the Jacobian ofthe function ( ¯ f k +1 (¯ z, c ) , . . . , ¯ f n (¯ z, c )) , evaluated at ¯ z = ( z k +1 , . . . , z n ). This function corresponds to the system called reduced inthe main text.By Proposition 7 the sign of the determinant of this Jacobian at ¯ z is exactly the signof the determinant of the Jacobian of the system S f = (cid:40) x i − ψ i (¯ x, c ) = 0 i = 1 , . . . , kf i ( x ) = 0 i = k + 1 , . . . , n evaluated at z . Indeed, the process leading from this system to the reduced system consistsof successive eliminations with i = j = 1 and the derivative of the eliminated functioncorresponding to the eliminated variable is 1 (and thus positive).Note that the reduced system has n − k variables. Let J ( S f ) denote the Jacobian of S f and let z be an equilibrium point with total amounts c . We conclude that: Result 13. With the notation introduced above: • If n − k is odd and det( J ( S f )( z )) > , then z is unstable. • If n − k is even and det( J ( S f )( z )) < , then z is unstable. C.2 Instability in the multistationary motifs Here we show how to apply Result 13 to each of the multistationary motifs: (f), (g), (i),(k), (l). Using Proposition 7, we find for all motifs but Motif (l) thatsign( ϕ (cid:48) ( z i )) = sign(det( J ( S f )( z ))) , where z i is the variable of ϕ (typically an intermediate complex Y ) and n − k is even. ForMotif (l), n − k is odd and there is a change of sign. Hence, using Result 13, the steadystate z is unstable whenever ϕ is decreasing. In particular, let (cid:15) = ± (cid:15) = +1or (cid:15) = − 1) such that (cid:15) · sign( ϕ (cid:48) ( z i )) = sign(det( J ( S f )( z ))) . Remark 4. The sign of the determinant of a matrix remains unchanged if a linear com-bination of rows are added to another row (in fact the determinant does not change). Ifwe multiply a row by a negative number, the determinant changes sign. For example, ifthe equation − ( b E + c E ) X + a E ES = 0 is transformed into X − ηES = 0, then the signof the determinant changes. If two such equations are transformed in this way, then thedeterminant remains with unchanged sign. In the sequel, the sign remains unchanged if thenumber of transformations of this type is even, or equivalently, if the number of constants δ ∗ , η ∗ is even.In the sequel, elimination is tracked using a table as in the main text and (cid:15) = (cid:81) n − l =1 (cid:15) l .The column ‘Behavior’ in the table shows ( i, j, s ) where i , respectively j , are the indices ofthe equation, respectively the variable, that iteratively are being eliminated and s indicateswhether ¯ f i ( f i after substitution of the previous eliminated variables) is increasing ( s is +)or decreasing ( s is − ) as function of x j .Motif (f) is covered in the main text and it is thus skipped here. Motif (g). Consider the variables in the following order( x , x , x , x , x , x , x , x , x ) = ( E, F, S , X , X , S , S , Y , Y ) . The variables of the reduced system are X , X , S , S , Y , Y and the variables eliminatedusing the conservation laws are E, F, S . The sign of the determinant of the Jacobian ofthe system S f is the same as the sign of the determinant of the Jacobian of the system f ( x ) = E + X + X − E f ( x ) = X − η ES f ( x ) = F + Y + Y − F f ( x ) = X − µ Y f ( x ) = S + S + S + X + X + Y + Y − S f ( x ) = X − µ Y f ( x ) = X − η ES f ( x ) = Y − δ F S f ( x ) = Y − δ F S . With the notation introduced above, ¯ x = ( X , X , S , S , Y , Y ) and ψ (¯ x ) = − X − X + E, ψ (¯ x ) = − Y − Y + F , ψ (¯ x ) = − S − S − X − X − Y − Y + S. Here n = 9 and k = 3, such that n − k = 6 is even. The function ϕ is equation f withall variables but x = S eliminated. The eliminations we performed in Section B.2 (Motif(g)) are summarized in the following table: l Elimination Behavior (cid:15) l l Elimination Behavior (cid:15) l f , E ) (1 , , +) + 6 ( f , S ) (2 , , − ) +2 ( f , F ) (1 , , +) + 7 ( f , S ) (3 , , − ) +4 ( f , X ) (4 , , +) + 8 ( f , Y ) (3 , , +) +5 ( f , X ) (4 , , +) + 9 ( f , Y ) (2 , , +) +Since (cid:15) = 1, we conclude that sign( ϕ (cid:48) ( z )) = sign(det( J ( S f )( z ))) for any equilibriumpoint z with z = S . Motif (i). Consider the variables in the following order( x , x , x , x , x , x , x , x , x , x ) = ( E, F, S , P , X , X , S , P , Y , Y ) . The variables of the reduced system are X , X , S , P , Y , Y and the variables eliminatedusing the conservation laws are E, F, S , P . The sign of the determinant of the Jacobianof the system S f is the same as the sign of the determinant of the Jacobian of the system f ( x ) = E + X + X − E f ( x ) = X − η ES f ( x ) = X − µ Y f ( x ) = F + Y + Y − F f ( x ) = X − η EP f ( x ) = Y − δ F S f ( x ) = S + S + X + Y − S f ( x ) = X − µ Y f ( x ) = Y − δ F P f ( x ) = P + P + X + Y − P Let ¯ x = ( X , X , S , P , Y , Y ). Then ψ (¯ x ) = − X − X + E ψ (¯ x ) = − S − X − Y + Sψ (¯ x ) = − Y − Y + F ψ (¯ x ) = − P − X − Y + P . Here n = 10 and k = 4, such that n − k = 6 is even. The function ϕ is equation f with allvariables but x = Y eliminated. The eliminations we performed in Section B.3 (Motif(i)) are summarized in the following table: l Elimination Behavior (cid:15) l l Elimination Behavior (cid:15) l f , E ) (1 , , +) + 6 ( f , P ) (5 , , − ) − f , F ) (1 , , +) + 7 ( f , S ) (3 , , − ) − f , X ) (5 , , +) + 8 ( f , P ) (3 , , − ) − f , X ) (5 , , +) + 9 ( f , Y ) (1 , , +) +5 ( f , S ) (5 , , − ) − Note that the last elimination f corresponds to ϕ ( Y , Y ) − S = 0, which is increasing in Y . Since (cid:15) = (cid:81) l =1 (cid:15) l = 1, we conclude that sign( ϕ (cid:48) ( z )) = sign(det( J ( S f )( z ))) for anyequilibrium point z and z = Y . Motif (k). Consider the variables in the following order( x , x , x , x , x , x , x , x , x , x ) = ( E, F, S , P , X , X , S , P , Y , Y ) . The variables of the reduced system are X , X , S , P , Y , Y and the variables eliminatedusing the conservation laws are E, F, S , P . The sign of the determinant of the Jacobianof the system S f is the same as the sign of the determinant of the Jacobian of the system f ( x ) = E + X − E f ( x ) = X − η ES f ( x ) = X − µ Y f ( x ) = F + Y + Y − F f ( x ) = X − η S P f ( x ) = Y − δ F S f ( x ) = S + S + X + Y + X − S f ( x ) = X − µ Y f ( x ) = Y − δ F P .f ( x ) = P + P + X + Y − P If ¯ x = ( X , X , S , P , Y , Y ) then ψ (¯ x ) = − X + E ψ (¯ x ) = − S − X − Y − X + Sψ (¯ x ) = − Y − Y + F ψ (¯ x ) = − P − X − Y + P . Here n = 10 and k = 4, such that n − k = 6 is even. The function ϕ is equation f with allvariables but x = Y eliminated. The eliminations we performed in Section B.4 (Motif(k)) are summarized in the following table: l Elimination Behavior (cid:15) l l Elimination Behavior (cid:15) l f , E ) (1 , , +) + 6 ( f , P ) (5 , , − ) − f , F ) (1 , , +) + 7 ( f , S ) (3 , , − ) − f , X ) (5 , , +) + 8 ( f , P ) (3 , , − ) − f , X ) (5 , , +) + 9 ( f , Y ) (1 , , +) +5 ( f , S ) (5 , , − ) − Note that the last elimination f corresponds to Φ( Y , Y ) − S = 0. Since (cid:15) = (cid:81) l =1 (cid:15) l = 1, weconclude that sign( ϕ (cid:48) ( z )) = sign(det( J ( S f )( z ))) for any equilibrium point z and z = Y . Motif (l). Consider the following order of the variables( x , x , x , x , x , x , x , x , x , x , x , x ) = ( E, F , F , S , P , X , X , X , S , P , Y , Y ) . The variables of the reduced system are X , X , X , S , P , Y , Y and the variables elimi-nated using the conservation laws are E, F , F , S , P . The sign of the determinant of theJacobian of the system S f is opposite the sign of the determinant of the Jacobian of thesystem f ( x ) = E + X + X − E f ( x ) = X − η S P f ( x ) = F + Y − F f ( x ) = X − η EP f ( x ) = F + Y − F f ( x ) = X − µ Y f ( x ) = S + S + X + Y + X − S f ( x ) = X − ( µ /µ ) X − µ Y f ( x ) = P + P + X + Y + X − P f ( x ) = Y − δ F S f ( x ) = X − η ES f ( x ) = Y − δ F P . Indeed, the reason for the change in sign comes from Remark 4, since there is an odd numberof equations that are rearranged by multiplication of negative numbers (corresponding to η , η , η , δ , δ ).Let ¯ x = ( X , X , X , S , P , Y , Y ). Then ψ (¯ x ) = − X − X + E ψ (¯ x ) = − S − X − Y − X + Sψ (¯ x ) = − Y + F ψ (¯ x ) = − P − X − Y − X + Pψ (¯ x ) = − Y + F . Here n = 12 and k = 5, such that n − k = 7 is odd. In this case, we should see that (cid:15) = 1.The function ϕ is equation f with all variables but x = Y eliminated. The eliminationswe performed in Section B.4 (Motif (l)) are summarized in the following table: l Elimination Behavior (cid:15) l l Elimination Behavior (cid:15) l f , E ) (1 , , +) + 7 ( f , X ) (4 , , +) − f , F ) (1 , , +) + 8 ( f , Y ) (4 , , − ) +3 ( f , F ) (1 , , +) + 9 ( f , P ) (4 , , − ) +4 ( f , X ) (6 , , +) − 10 ( f , S ) (3 , , − ) − f , S ) (7 , , − ) − 11 ( f , X ) (2 , , +) − f , P ) (4 , , − ) − The last elimination f corresponds to ϕ P ( Y , X ) − S = 0, which is increasing in X .Since (cid:15) = (cid:81) l =1 (cid:15) l = 1, we conclude that sign( ϕ (cid:48) ( z )) = − sign(det( J ( S f )( z ))) for anyequilibrium point z and z = Y . Together with the fact that n − k is odd, the steadystates for which ϕ (cid:48) < C.3 The monostationary motifs For Motifs (a)-(d), (e), (h) and (j), the above procedure is non-conclusive. We only showthis for Motif (j), since for the other motifs we can prove that the steady state is asymp-totically stable using the Routh-Hurwitz criterion. Motif (j). Consider the variables in the following order( x , x , x , x , x , x , x , x , x , x , x ) = ( E, F , F , S , P , X , X , S , P , Y , Y ) . The variables of the reduced system are X , X , S , P , Y , Y and the variables eliminatedusing the conservation laws are E, F , F , S , P . The sign of the determinant of the Jaco-bian of the system S f is the same as the sign of the determinant of the Jacobian of thesystem f ( x ) = E + X − E f ( x ) = X − η ES f ( x ) = X − µ Y f ( x ) = F + Y − F f ( x ) = X − η S P f ( x ) = Y − δ F S f ( x ) = F + Y − F f ( x ) = X − µ Y f ( x ) = Y − δ F P .f ( x ) = S + S + X + Y + X − Sf ( x ) = P + P + X + Y − P Here n = 11 and k = 5, such that n − k = 6 is even. The function ϕ is equation f with allvariables but x = Y eliminated. The eliminations we performed in Section B.4 (Motif(j)) are summarized in the following table: l Elimination Behavior (cid:15) l l Elimination Behavior (cid:15) l f , E ) (1 , , +) + 6 ( f , S ) (3 , , − ) − f , F ) (1 , , +) + 7 ( f , S ) (4 , , − ) − f , F ) (1 , , +) + 8 ( f , P ) (3 , , − ) − f , X ) (5 , , +) + 9 ( f , P ) (3 , , − ) − f , X ) (5 , , +) + 10 ( f , Y ) (1 , , +) +Since (cid:15) = (cid:81) l =1 (cid:15) l = 1, we conclude that sign( ϕ (cid:48) ( z )) = sign(det( J ( S f )( z ))) for any equi-librium point z and z = Y . Since n − k is even, stable steady states have positivedeterminant, and since ϕ is always increasing, nothing can be concluded. Routh-Hurwitz. We have computationally checked that Motifs (a)-(d), (e) and (h) haveasymptotically stable BMSSs. The Routh-Hurwitz criterion establishes when all roots of apolynomial have real negative parts from a condition on the coefficients of the polynomial.Specifically, consider a real polynomial p ( z ) = α z n + α z n − + · · · + α n − z + α n , where we can assume that α > 0. One constructs the following matrix: A = α α α . . . . . . α α α α . . . α α α . . . α α α . . . ... ... ... ... ... α n . That is, the entry ( i, j ) of the matrix is α i +2( j − i ) with the convention that if i + 2( j − i ) > n or i + 2( j − i ) < 0, then α i +2( j − i ) = 0. The criterion states that if all leading principalminors of the matrix have positive sign, then all roots of the polynomial p ( z ) have negativereal parts.In our case, the polynomial to be analyzed is the characteristic polynomial of theJacobian of the reduced system at a steady state. For example, by eliminating the variables E and S , the reduced system of Motif (b) is˙ X = − ( b + c ) X + a ( E − X − Y )( S − S − X − Y )˙ Y = − ( b + c ) Y + a ( E − X − Y ) S ˙ S = c X + b Y − a ( E − X − Y ) S . The Jacobian can easily be computed using a program that handles symbolic compu-tations (e.g. Mathematica T M ). For Motif (b) the Jacobian J is − b − c − a ( E + S − S − X − Y ) − a ( E + S − S − X − Y ) − a ( E − X − Y ) − a S − b − c − a S a ( E − X − Y ) c + a S b + a S − a ( E − X − Y ) Note that for any BMSS, the terms E + S − S − X − Y and E − X − Y are positive. Theleading principal minors of the Routh-Hurwitz matrix are polynomials in the entries of J .The task of computationally determining the sign of these minors is greatly simplified ifwe substitute back E = E + X + Y , and S = S + S + X + Y into J and write: J = − b − c − a ( E + S ) − a ( E + S ) − a E − a S − b − c − a S a Ec + a S b + a S − a E . The entries of the matrix J are now in the set of variables V = { E, X, Y, S , S , a , a , b , b , c , c } and all variables in V take positive values at any BMSS. Thus, each of the entries of J has afixed sign for any BMSS. The coefficients of the characteristic polynomial are polynomialsin V and thus, the leading principal minors of the Routh-Hurwitz matrix are also polyno-mials in V . If the coefficients of these polynomials are all positive, then we are guaranteedthat they take positive values for any choice of positive rates and any set of positive concen-trations. This computation can easily be implemented and checked with Mathematica T M .If we used the matrix J involving the terms E + S − S − X − Y and E − X − Y , thenthe corresponding polynomials take values in { E, S, X, Y, S , a , a , b , b , c , c } but theircoefficients are no longer positive. Positivity checking requires a convenient grouping ofterms.This procedure shows stability for all motifs but Motif (j) in which case some negativeterms in the expansion of some of the minors occur. These are computationally difficult tohandle and we have not been able to show that the minors all are positive. However, byrandom generation of values, it seems that the steady state is stable, though it remains tobe proven. D Proofs Proposition 1. Figure 7 illustrates the different cases of the proposition. Recall that Γ = { Y ∈ Γ | f ( Y ) ∈ Γ } and ∆ = µ F − E . The derivative of f is f (cid:48) ( y ) = η ( δ EF − δ µ F y + ( δ − η ) µ y ) µ ( δ ( F − y ) + η y ) The denominator is always positive. The numerator is a degree two polynomial in y which mighthave positive real roots and is positive for y = 0. We consider y ≤ ξ = min( F , E/µ ).The numerator evaluated at y = F is η ( δ E − ( δ + η ) µ F ) F . It is positive if and only ifΛ = (1 + η /δ ) µ F − E ≤ . If this is the case, then ξ = F and the numerator is bounded from below by η µ ( δ ( F − y ) + η ( F − y )) > 0. Hence both the numerator and f (cid:48) ( y ) are positive for all y < ξ .Therefore, for Λ ≤ f is increasing for all y ∈ Γ = [0 , F ) (see Figure 7(A)). Since f (0) = 0, the condition f ( Y ) ∈ Γ is equivalent to Y < f − ( ξ ) and hence Γ = [0 , min( F , f − ( ξ )).This proves (i).To prove (ii), note that f (cid:48) ( F ) < > 0. Further, f (cid:48) ( E/µ ) < E ( δ − η ) − δ µ F < 0. This is clearly the case when ∆ ≥ 0. If ∆ < 0, then using Λ > 0, it follows that E ( δ − η ) − δ µ F < η ( − E + µ F ) < . In either case, f (cid:48) ( E/µ ) < 0. It follows that if Λ > 0, then, whatever the sign of ∆ , we have f (cid:48) ( ξ ) < f (cid:48) has exactly one positive real root α in Γ . Thus, f is afunction that increases up to y = α , and decreases after α (see Figure 7). We have that f ( E/µ ) = 0and hence, if ∆ > ξ = E/µ , f ( ξ ) = 0, and the function f starts and ends in zero in Γ .Oppositely, if ∆ ≤ ξ = F and f ( ξ ) = − ∆ /µ > 0. These differences are depicted inFigure 7(B-C).The three cases (ii)(a)-(c) of the proposition follow from the determination of Γ. The idea isdepicted in Figure 8. For f ( Y ) ∈ Γ we require that f ( Y ) < ξ . The maximal value of f in Γ is f ( α ) which is strictly smaller than E/µ because µ Y + µ Y = µ Y + µ f ( Y ) = E . Thus, wehave • If f ( α ) < F , then for all Y ∈ Γ , we have f ( Y ) ≤ f ( α ) < ξ = min( F , E /µ ). Thus,Γ = [0 , ξ ), which corresponds to case (c) of the proposition. • If f ( α ) ≥ F , then the horizontal line Y = ξ = F intersects the graph of f over [0 , E/µ )in two points corresponding to Y -values α ≤ α ≤ α < E/µ ( α = α if ξ = f ( α )). Wealso have that α < ξ , since α ≤ α < ξ .In the latter case, if Y ∈ [0 , α ) ∪ ( α , E/µ ), then f ( Y ) ∈ Γ , while f ( Y ) lies outside Γ if Y ∈ [ α , α ]. Therefore, if f ( α ) ≥ F , we haveΓ = ([0 , α ) ∪ ( α , E/µ )) ∩ [0 , ξ ) = (cid:40) [0 , α ) if ξ ≤ α [0 , α ) ∪ ( α , ξ ) if ξ > α . The first case happens if and only if ξ = F and F ≤ α . Since α ≤ F , this condition isequivalent to ∆ ≤ − ∆ /µ = f ( F ) ≥ f ( α ) = F , or 0 ≤ E − ( µ F + µ F ). This proves( a ). The second case is equivalent to either ξ = E/µ or ξ = F > α . Proceeding as above, thisis equivalent to 0 < µ F + µ F − E , which proves case (b). Proposition 2. Assume that Λ > 0, that is Proposition 1(ii) applies. Let ¯ ϕ ( Y ) = ϕ ( Y , f ( Y )) sothat ϕ ( Y ) = ¯ ϕ ( Y ) + ϕ ( f ( Y )). The derivative of ϕ with respect to Y is ϕ (cid:48) ( Y ) = ∂ ¯ ϕ ∂Y ( Y ) + ∂ϕ ∂Y ( f ( Y )) ∂f∂Y ( Y ) with ∂ϕ ∂Y = (1 + µ ) + F δ ( F − Y ) . The term ∂ϕ ∂Y is the only one that depends on F . Note in addition that this term is decreasing in F for Y < F .The variables S , Y , X are all increasing functions of Y . As for S = µ Y η ( E − µ Y − µ f ( Y )) , wehave that the derivative of µ Y + µ f ( Y ) is δ η ( EF − µ Y ) + δ µ ( F − Y ) ( η Y + δ ( F − Y )) , (76)which is positive for Y < F , E/µ . Thus, S is also increasing in Y and hence ∂ϕ ∂Y ( Y ) > F = f ( α ) there exist(many) values of Y ∈ Γ (cid:48)(cid:48) for which ϕ (cid:48) ( Y ) < 0. Fix one such Y . Since ϕ (cid:48) is continuous in F , thenthere exists an (cid:15) > ϕ (cid:48) ( Y ) < F ∈ ( f ( α ) − (cid:15), f ( α ) + (cid:15) ). This ensures that forvalues of F ≥ f ( α ) close to f ( α ), we have multistationarity.On the other hand, the term ∂ϕ ∂Y ( f ( Y )) is positive for any Y ∈ Γ. Therefore, if ϕ (cid:48) vanishes forsome Y , it must be that Y > α , where f is decreasing. In the interval I = ( α, ξ ), − ∂f /∂Y is pos-itive, bounded from above and independent of F . It follows from the expression for ∂f /∂Y givenin the proof of Proposition 1. Further, for Y ∈ I , ∂ϕ /∂Y is positive (stated after equation (76)),bounded from below and independent of F . The latter two statements follow from Result 1 andequation (76) upon differentiation (where (76) is used to differentiate S ).Finally, when F tends to infinity, ∂ϕ /∂Y tends to zero uniformly for all Y ∈ I , since Y = f ( Y ) ≤ f ( α ) is independent of F . Putting it all together, we find that for F large ∂ϕ ∂Y ( Y ) > − ∂ϕ ∂Y ( f ( Y )) ∂f∂Y ( Y ) (77)for all Y ∈ I . Consequently the function ϕ is increasing in all Γ and there cannot be multistation-arity.Therefore, for values of F above f ( α ) but ‘close’ to f ( α ) there is multistationarity, and forlarge values of F , there is monostationarity. If F (cid:48) is such that ϕ ( Y ) is increasing for some value Y , then for any F ≥ F (cid:48) , ϕ ( Y ) must be increasing too (as ϕ (cid:48) ( Y ) continues being positive whenincreasing F , cf. (77)).This implies that the two regions (mono- versus multistationarity) are separated by a certainboundary value of F , say M , where M is the infimum of all F such that ϕ (cid:48) ( Y ) > Y ∈ Γ. For F = M , ϕ is also strictly increasing: The derivative ϕ (cid:48) might be zero (as M is theinfimum) but this can only happen in a finite number of points because the derivative is a non-zerorational function. Therefore, we see that for all F ≥ M , there is only one BMSS, while for all M > F > f ( α ) there is multistationarity. Proposition 3. Recall that σ = ( µ − µ )( µ δ η − µ δ η ) and ψ = ( ϕ , ϕ ). Let us compute det( J ψ ( Y , Y )). We have, writing C F = F − Y − Y and C E = E − µ Y − µ Y , ∂ϕ ∂Y = 1 + µ + Fδ C F − Y δ C F + µ Eη C E − µ µ Y η C E ∂ϕ ∂Y = Y δ C F + µ µ Y η C E ∂ϕ ∂Y = 1 + µ + Fδ C F − Y δ C F + µ Eη C E − µ µ Y η C E ∂ϕ ∂Y = Y δ C F + µ µ Y η C E Let C = ∂ϕ ∂Y − − µ > D = ∂ϕ ∂Y − − µ > 0. Then, the determinant is given asdet( J ψ ( Y , Y )) = (1 + µ )(1 + µ ) + (1 + µ ) C + (1 + µ ) D + Fδ δ C F + µ µ Eη η C E + A + B with A = µ ( E F − µ Y F − Y E ) δ η C E C F = µ δ η C E C F (cid:16) Y C F + µ Y C E + ( µ − µ ) Y Y C E C F (cid:17) B = µ ( E F − µ Y F − Y E ) δ η C E C F = µ δ η C E C F (cid:16) Y C F + µ Y C E + ( µ − µ ) Y Y C E C F (cid:17) where in the last equality of both equations, we substitute F by C F + Y + Y and E by C E + µ Y + µ Y . Hence, det( J ψ ( Y , Y )) is a sum of positive terms together with N = Y Y δ δ η η C E C F ( µ δ η ( µ − µ ) + µ δ η ( µ − µ ))Therefore, if σ = ( µ − µ )( µ δ η − µ δ η ) ≥ J ψ ( Y , Y )) > ϕ is always increasing. This proves (i).Let us assume now that σ < 0. We make the following observation: the only negative term N of det( J ψ ( Y , Y )) is also the term with denominator of highest degree in C E , C F . Further,all numerators in the expression of det( J ψ ( Y , Y )) above are bounded. Thus, by letting C E , C F simultaneously tend to zero, we can obtain a negative determinant.Since C F = F − Y − Y and C E = E − µ Y − µ Y , it is possible to make them smallsimultaneously by varying Y , Y , if and only if the two lines r F : F − Y − Y = 0 and r E : E − µ Y − µ Y = 0 intersect for valid values of Y , Y . The intersection of the two lines is the point( Y , Y ) with Y = µ F − Eµ − µ , Y = µ F − Eµ − µ . If µ − µ > 0, then for positive intersection values Y , Y we require µ F < E < µ F . If µ − µ > µ F < E < µ F .Hence, if these conditions above are satisfied, we are guaranteed the existence of values of Y , Y for which the determinant is negative and thus multistationarity occurs. Further, if C E , C F aresmall then S, P must be large.Assume finally that σ < (a) µ − µ > µ F ≥ E or E ≥ µ F .(b) µ − µ > µ F ≥ E or E ≥ µ F .In this case multistationarity cannot occur. Assume that (a) holds (case (b) is similar). Since µ − µ > 0, we have µ F > µ F and: • If µ F > µ F ≥ E , then using − Y k E ≥ − µ k Y k F , we have EF − µ Y F − Y E ≥ EF − µ Y F − µ Y F = F C E > EF − Y E − µ Y F ≥ EF − µ Y F − µ Y F = F C E > . • If E ≥ µ F > µ F , then using − µ k Y k F ≥ − Y k E , we have EF − µ Y F − Y E ≥ EF − Y E − Y E = EC F > EF − Y E − µ Y F ≥ EF − Y E − Y E = EC F > . Therefore, both A and B are positive and the proposition is proved. Proposition 4. Recall that f is a function of Y defined in equation (65) by elimination of Y from S = Φ( Y , Y ). Thus, we have that f (cid:48) = − Φ / Φ , whereΦ := ∂ Φ /∂Y = µ + Y δ ( F − Y − Y ) > , Φ := ∂ Φ /∂Y = 1 + µ + µ Eη ( E − µ Y ) + F − Y δ ( F − Y − Y ) > . Recall that ∂P ∂Y = µ δ η f ( f ( F − f − Y ) − f (cid:48) Y ( F − Y )) , and ∂P ∂Y = F − f + f (cid:48) Y δ ( F − f − Y ) . ∂P ∂Y < F − f + f (cid:48) Y < 0, that is0 > Φ ( F − Y ) − Φ Y = (1 + µ )( F − Y ) − µ Y + µ E ( F − Y ) η ( E − µ Y ) + Fδ ( F − Y − Y ) . (78)Since the last two summands are positive for valid values of Y , Y , a necessary condition for ∂P ∂Y < C µ )( F − Y ) − µ Y < 0. Note that if 1 + µ > µ , then the condition cannot be fulfilledand P is always an increasing function of Y for any E, F , S .Since f (cid:48) < 0, we require C F − Y − Y < ∂P ∂Y < C , C P , S are given, thefollowing are restrictions on Y , Y : Y + Y < F , µ Y < E, (1 + µ ) Y + µ Y < S, µ Y < P . Let σ = min( S, P ). We find that µ Y < σ and hence F − Y − Y > F − E/µ − σ/µ and (1 + µ )( F − Y ) − µ Y > (1 + µ )( F − E/µ ) − σ. Let σ µ = min(1 + µ , µ / F > E/µ and σ µ ( F − E/µ ) > σ , then we are guaranteed thatboth C , C > Proposition 5. Note that for Φ (i.e. S ) to tend to + ∞ for any fixed Y ∈ Γ = [0 , min( S/µ , F )) ⊆ [0 , F ), either Y must tend to E/µ or to F − Y . The first case arises if and only if Y ≤ F − E/µ .(i) Assume that Y > F − E/µ . Consider the functions Φ , Φ from the proof of Proposition 4and the expression for ∂P ∂Y and ∂P ∂Y : The derivative f (cid:48) = − Φ / Φ tends to − S tends to infinity(and Y tends to F − Y ). Likewise, it follows that f ( F − f − Y ) − f (cid:48) Y ( F − Y ) tends to zero,and hence ∂P ∂Y also tends to zero. From equation (78) in the previous proof it follows that ∂P ∂Y > Y approaches F − Y , that is as S becomes large (the last term is positive and dominates theother terms). Hence the derivative of ϕ ( Y ) becomes positive: ϕ (cid:48) ( Y ) = 1 + µ + ∂P ∂Y + ∂P ∂Y > = µ F − E > Y ≤ F − E/µ . As S tends to infinity, Y tends to E/µ . The limit curve ϕ ∞ is ϕ ∞ ( Y ) = (1 + µ ) Y + µ δ (∆ − Y ) Y η E/µ + Y δ (∆ − Y )for ∆ = ∆ /µ . Existence of values of Y for which ϕ (cid:48)∞ ( Y ) < p ( Y ) := δ (∆ − Y ) (1 + µ + β (∆ − Y )) + ∆ < β = µ µ δ /η E . The function p ( Y ) is a polynomial of degree 3. For Y = 0, it is positiveand when Y tends to infinity, the polynomial tends to −∞ . Therefore, it takes negative values in Y ∈ [0 , ∆ ) if and only if there is a root in this interval. At Y = ∆ it is positive, and hence, thereis at least one root after ∆ . The derivative with respect to Y is p (cid:48) ( Y ) = − δ (∆ − Y )(1 + µ + β (2∆ − Y )) . One zero of the derivative is Y = ∆ , while the other is γ = (1 + µ + 2 β ∆ ) / β . Note that at ∆ , p (∆ ) is positive. Therefore, for negative values of p to occur between 0 and ∆ , we need γ < ∆ and p ( γ ) < γ < ∆ if and only if ∆ > (1 + µ ) /β , in which case ∆ is a maximum and γ is aminimum. Evaluating p in γ gives the following condition:∆ = 27 µ µ δ η E ( µ F − E ) − δ ( δ µ µ F − ( µ δ + η (1 + µ )) E ) < , which concludes the proof of (ii). Proposition 6. It is convenient to write ϕ as a function of X , Y , i.e., without substituting X by f ( Y ). The derivative of ϕ with respect to Y takes the form ϕ (cid:48) = µ ( η E + δ η X ( µ F − E )) η η ( E − µ Y ) + 1 + µ + F δ ( F − Y ) + (cid:18) µ δ η ( F − Y ) η η ( E − µ Y ) + 1 (cid:19) f (cid:48) with only the first term susceptible of being negative since f is an increasing function. If we cancontrol f (cid:48) by varying some total amounts we might be able to make the derivative ϕ (cid:48) negative.We will see that there exist values of X , F , and E for which ϕ (cid:48) as a function of the two variables Y and X is negative. Then, by defining P = ϕ P ( Y , X ), we obtain regions of multistationarity.Note that this is equivalent to the approach taken in Motif (i), when we studied the Jacobian of( ϕ S , ϕ P ) rather than the derivative of ϕ .Using the implicit function theorem, the derivative of f with respect to Y is f (cid:48) = − ∂ϕ P /∂Y ∂ϕ P /∂X .We compute the partial derivatives of P , P , X , Y , X for each term in the numerator and thedenominator. First of all we repeat some definitions, X = Φ X ( Y , X ) = δ η X ( F − Y )( E − µ Y ) η Y + δ η X ( F − Y ) , and Y = Φ Y ( Y , X ) = µ − X + µ − Φ X . It follows that ∂ϕ P ∂X = ∂ ( P + P ) ∂X + 1 + µ − + (1 + µ − ) ∂ Φ X ∂X , ∂ϕ P ∂Y = ∂ ( P + P ) ∂Y + (1 + µ − ) ∂ Φ X ∂Y . For the denominator we have: ∂P ∂X = δ ( F − Y ) η Y , ∂ Φ X ∂X = δ η η Y ( F − Y )( E − µ Y )( η Y + δ η X ( F − Y )) ,∂P ∂X = F δ ( F − Φ Y ) · ( µ − + µ − ∂ Φ X ∂X ) . For the numerator, we have − ∂P ∂Y = δ F X η Y − ∂ Φ X ∂Y = δ η X µ ( F − Y ) η Y + δ η X ( F − Y ) + η F ( E − µ Y )( η Y + δ η X ( F − Y )) . − ∂P ∂Y = − F δ µ ( F − Φ Y ) ∂ Φ X ∂Y Now fix Y , F and let X , F , E tend to + ∞ . Assume that E goes slower than X , for exampleput X = E p for some p > 0. Let F − Φ Y = K be constant such that F = K + µ − X + µ − Φ X and X go at the same rate towards infinity. Note that we can vary X , F , E without changing Y , F or restricting their range of definition. With these assumptions the numerator of f (cid:48) , − ∂φ P ∂Y , goes to infinity proportionally to X . Alsothe denominator of f (cid:48) , ∂φ P ∂X , goes to infinity proportionally to X . Consequently, f (cid:48) convergestowards a constant that depends on F , Y and K (in addition to the rate constants) f (cid:48) → K δ δ µ F η Y + µ µ µ . Going back to the expression for ϕ (cid:48) , we find that under the same assumptions, ϕ (cid:48) ≈ − µ δ η X η η E + 1 + µ + F δ ( F − Y ) + f (cid:48) , where terms that eventually vanish are not shown. It follows that for large E , X = E p and F = K + µ − X + µ − Φ X , the derivative ϕ (cid:48) eventually becomes negative for any choice of Y and F1