Identifying parameter regions for multistationarity

Mathematical modelling has become an established tool for studying the dynamics of biological systems. Current applications range from building models that reproduce quantitative data to identifying systems with predefined qualitative features, such as switching behaviour, bistability or oscillations. Mathematically, the latter question amounts to identifying parameter values associated with a given qualitative feature. We introduce a procedure to partition the parameter space of a parameterized system of ordinary differential equations into regions for which the system has a unique or multiple equilibria. The procedure is based on the computation of the Brouwer degree, and it creates a multivariate polynomial with parameter depending coefficients. The signs of the coefficients determine parameter regions with and without multistationarity. A particular strength of the procedure is the avoidance of numerical analysis and parameter sampling. The procedure consists of a number of steps. Each of these steps might be addressed algorithmically using various computer programs and available software, or manually. We demonstrate our procedure on several models of gene transcription and cell signalling, and show that in many cases we obtain a complete partitioning of the parameter space with respect to multistationarity.


Introduction
Mathematical models in the form of parameterized systems of ordinary differential equations (ODEs) are valuable tools in biology. Often, qualitative properties of the ODEs are associated with macroscopic biological properties and biological functions [1][2][3][4]. It is therefore important that we are able to analyse mathematical models with respect to their qualitative features and to understand when these properties arise in models. With the growing adaptation of differential equations in biology, an automated screening of ODE models for parameter dependent properties and discrimination of parameter regions with different properties would be a very useful tool for biology, and perhaps even more for synthetic biology [5]. Even though it is currently not conceivable how and if this task can be efficiently formalized, we view the procedure presented here as a first step in this direction.
Multistationarity, that is, the capacity of the system to rest in different positive equilibria depending on the initial state of the system, is an important qualitative property. Biologically, multistationarity is linked to cellular decision making and 'memory'-related on/off responses to graded input [2][3][4]. Consequently the existence of multiple equilibria is often a design objective in synthetic biology [6,7]. Various mathematical methods, developed in the context of reaction network theory, can be applied to decide whether multistationarity exists for some parameter values or not at all, or to pinpoint specific values for which it does occur [8][9][10][11][12][13][14][15][16][17]. Some of these methods are freely available as software tools [18,19].
It is a hard mathematical problem to delimit parameter regions for which multistationarity occurs. Often it is solved by numerical investigations and parameter sampling, guided by biological intuition or by case-by-case mathematical approaches. A general approach, in part numerical, is based on a certain bifurcation condition [15,16,20,21]. Alternatively, for polynomial ODEs, a decomposition of the parameter space into regions with different numbers of equilibria could be achieved by Cylindrical Algebraic Decomposition (a version of quantifier elimination) [22]. This method, however, scales very poorly and is thus only of limited help in biology, where models tend to be large in terms of the number of variables and parameters.
Here we present two new theoretical results pertaining to multistationarity (Theorem 1 and Corollary 2). The results are in the context of reaction network theory and generalize ideas in [23,24]. We consider a parameterized ODE system defined by a reaction network and compute a single polynomial in the species concentrations with coefficients depending on the parameters of the system. The theoretical results relate the capacity for multiple equilibria or a single equilibrium to the signs of the polynomial as a function of the parameters and the variables (concentrations).
The theoretical results apply to dissipative reaction networks (networks for which all trajectories eventually remain in a compact set) without boundary equilibria in stoichiometric compatibility classes with non-empty interior. These conditions are met in many reaction network models of molecular systems. We show by example that the results allow us to identify regions of the parameter space for which multiple equilibria exist and regions for which only one equilibrium exists. Subsequently this leads to the formulation of a general procedure for detecting regions of mono-and multistationarity. The procedure verifies the conditions of the theoretical results and further, calculates the before-mentioned polynomial. A key ingredient is the existence of a positive parameterization of the set of positive equilibria. Such a parameterization is known to exist for many classes of reaction networks, for example, systems with toric steady states [11] and post-translational modification systems [25,26].
The conditions of the procedure might be verified manually or algorithmically according to computational criteria. The algorithmic criteria are, however, only sufficient for the conditions to hold. For example, a basic condition is that of dissipativity. To our knowledge there is not a sufficient and necessary computational criterion for dissipativity, but several sufficient ones. If these fail, then the reaction network might still be dissipative, which might be verified by other means. By collecting the algorithmic criteria, the procedure can be formulated as a fully automated procedure (an algorithm) that partitions the parameter space without any manual intervention. The algorithm might however terminate indecisively if some of the criteria are not met. Table 1 shows two examples of reaction network motifs that occur frequently in intracellular signalling: a two-site protein modification by a kinase-phosphatase pair and a one-site modification of two proteins by the same kinase-phosphatase pair. These reaction networks are in the domain of the automated procedure and conditions for mono-and multistationarity can be found without any manual intervention. The conditions discriminating between a unique and multiple equilibria highlight a delicate relationship between the catalytic and Michaelis-Menten constants of the kinase and the phosphatase with the modified protein as a substrate (the k cand k M -values). If the condition for multiple equilibria is met, then multiple equilibria occur provided the total concentrations of kinase, phosphatase and substrate are in suitable ranges (values thereof can be computed as part of the procedure).
The paper has three main sections: a theoretical section, a section about the procedure and an application section. We close the paper with two brief sections discussing computational limitations, related work and future directions. In the theoretical section we first introduce notation and mathematical background material. We then give the theorem and the corollary that links the number of equilibria to the sign of the determinant of the Jacobian of a certain function, which is derived from the ODE system associated with a reaction network. In the Multiple: b1(κ) < 0, Unique: b1(κ) ≥ 0 and b2(κ) ≥ 0  second section we state the procedure, derive the algorithm and comment on the feasibility and verifiability of the conditions. Finally, in the application section we apply the procedure to several examples. The Supporting Information has six sections. All proofs are relegated to §1-4 together with background material. In §5 we elaborate further on how the conditions of the procedure/algorithm can be verified. In §6 we provide details of the algorithmic analysis of the examples in Table 1. Also we include a further monostationary example for illustration of the algorithm.

Theory
In this part of the manuscript we present the theoretical results. We start by introducing the basic formalism of reaction networks. Theorem 1, Corollary 1 and 2 below apply to dissipative networks without boundary equilibriaand concern the (non)existence of multiple equilibria in some stoichiometric compatibility class. Corollary 2 assumes the existence of a positive parameterization of the set of positive equilibria. Before stating the results these five concepts are formally defined.

Reaction networks.
A reaction network, or simply a network, consists of a set of species {X 1 , . . . , X n } and a set of reactions of the form: where α ij , β ij are non-negative integers. The left hand side is called the reactant, while the right hand side is called the product. We let N = (N ij ) ∈ R n× be the stoichiometric matrix of the network, defined as N ij = β ij − α ij , that is, the (i, j)-th entry encodes the net production of species X i in reaction R j . We refer to the 'running example' in Fig 1 for an illustration of the definitions. The concentrations of the species X 1 , . . . , X n are denoted by lower-case letters x 1 , . . . , x n and we let x = (x 1 , . . . , x n ). We denote by R n >0 (R n ≥0 ), the positive (non-negative) orthant in R n . The evolution of the concentrations with respect to time is modeled as an ODE system derived from a set of reaction rate functions. A reaction rate function for reaction R j is a C 1 -function Reaction network: Stoichiometric matrix N , vector of mass-action reaction rate functions v(x) and a matrix W such that W N = 0: The set V from Eq (4) for κ1 = 5, κ2 = 1, κ3 = 3 and three stoichiometric compatibility classes for c = 1, 2, 3: Figure 1: Running example. Example network with two species, X 1 and X 2 , and three reactions with mass-action kinetics.
v j : R n ≥0 → R ≥0 that models the (non-negative) speed of the reaction. We further assume that v j (x) = 0 ⇔ x i = 0 for some i such that α ij > 0, (2) that is, the reaction only takes place in the presence of all reactant species. We refer to the set of reaction rate functions as the kinetics.
A particular important example of a kinetics is that of mass-action kinetics. In this case the reaction rate functions are given by where κ j is a positive number called the reaction rate constant and we assume 0 0 = 1.
Other important examples are Michaelis-Menten kinetics and Hill kinetics. All three types of kinetics fulfil the assumption in Eq (2). For a choice of reaction rate functions v = (v 1 , . . . , v ), the ODE system modelling the species concentrations over time with initial condition x(0) = x 0 , iṡ Under assumption Eq (2), the orthants R n >0 and R n ≥0 are forward-invariant under f in Eq (21) [27, Theorem 5.6], [28,Section 16]. Forward-invariance implies that the solutions to the ODE system stays in R n >0 (resp. R n ≥0 ) for all positive times if the initial condition is in R n >0 (resp. R n ≥0 ). The trajectories of the ODEs in Eq (21) are confined to the so-called stoichiometric compatibility classes, which are defined as follows. Let s = rank(N ) be the rank of the network and d = n − s be the corank. Further, let W ∈ R d×n be any matrix of full rank d such that W N = 0, see Fig 1 for an example. This matrix is zero-dimensional if N has full rank n. For each c ∈ R d , there is an associated stoichiometric compatibility class defined as This set is empty if c / ∈ W (R n ≥0 ). The positive stoichiometric compatibility class is defined as the relative interior of P c , that is, the intersection of P c with the positive orthant: The sets P + c and P c are convex. Since by construction W x is conserved over time and determined by the initial condition, then P + c and P c are also forward-invariant. An equation of the form ω · x = c for some ω ∈ im(N ) ⊥ and c ∈ R is called a conservation relation. In particular, W x = c forms a system of d conservation relations.
For the running example in Fig 1, the rank of the network is s = 1 and the corank is d = 1. The matrix W in the figure leads to the conservation relation x 1 +x 2 = c. Here the stoichiometric compatibility class P c has non-empty interior, that is, P + c = ∅, if and only if c > 0. In the following, to ease the notation, we implicitly assume a reaction network comes with a kinetics (a set of reaction rate functions) and the associated ODE system.
Dissipative and conservative reaction networks. A reaction network is dissipative if, for all stoichiometric compatibility classes P c , there exists a compact set where the trajectories of P c eventually enter (see §3.2 in the Supporting Information). A reaction network is conservative if there exists a conservation relation with only positive coefficients, or, equivalently, if for all species X i there is a conservation relation such that the coefficient of x i is positive and all other coefficients are non-negative. This is equivalent to the stoichiometric compatibility classes being compact sets [29]. Hence, in particular, a conservative reaction network is dissipative because we can choose the attracting compact set to be the stoichiometric compatibility class itself. Because of the conservation relation x 1 + x 2 = c, the reaction network of the running example is conservative.

Equilibria.
Given the ODE in Eq (21), the set of non-negative equilibria is the set of points for which f (x) vanishes: We are interested in the positive equilibria in each stoichiometric compatibility class, that is, in the set V ∩ P + c . Generically, this set consists of isolated points obtained as the simultaneous positive solutions to the equations We introduce some definitions: a network admits multiple equilibria (or is multistationary) if there exists c ∈ R d such that V ∩ P + c contains at least two points, that is, the system in Eq (5) has at least two positive solutions. Equilibria belonging to V ∩ P c but not to V ∩ P + c for some c are boundary equilibria. A boundary equilibrium has at least one coordinate equal to zero.
The function ϕ c (x).
Some of the n equations in the system f (x) = 0 might be redundant. Indeed, every vector ω ∈ im(N ) ⊥ fulfils ω · f (x) = 0, and hence gives a linear relation among the entries of f (x). As a consequence, there are (at least) as many independent linear relations as rows of W , that is, d, and there are at most s = n − d linearly independent equations in the system f (x) = 0. Thus d of the equations are redundant. By removing these from f (x) = 0, the system in Eq (5) becomes a system of n equations in n variables.
In order to systematically choose d equations to remove, we proceed as follows. We choose the matrix of conservation relations W ∈ R d×n to be row reduced and let i 1 , . . . , i d be the indices of the first non-zero coordinate of each row. Then the scalar product of the j-th row of W with f (x) can be used to express f i j (x) as a linear combination of the entries of f (x) with indices different from i 1 , . . . , i d . It follows that the equations For the running example in Fig 1 the matrix W is already row reduced with i 1 = 1. Hence ϕ c is obtained by replacing f 1 (x) with x 1 + x 2 − c: As the function ϕ c (x) is obtained by replacing redundant equations in f (x) = 0 with equations defining P c , we have Consequently, a network admits multiple equilibria if the equation ϕ c (x) = 0 has at least two positive solutions for some c ∈ R d .
A theorem for unique and multiple equilibria.
Let M (x) ∈ R n×n be the Jacobian matrix of ϕ c (x), that is, the matrix with (i, j)-th entry equal to the partial derivative of ϕ c,i (x) with respect to x j . The matrix M (x) does not depend on c, see Eq (6).
We say that an equilibrium x * ∈ V ∩ P c is non-degenerate if the Jacobian of ϕ c at x * , M (x * ), is non-singular, that is, if det(M (x * )) = 0 [10].
Theorem 1 (Unique and multiple equilibria). Assume the reaction rate functions fulfil Eq (2), let s = rank(N ) and let P c be a stoichiometric compatibility class such that P + c = ∅, where c ∈ R d . Further, assume that (i) The network is dissipative. (ii) There are no boundary equilibria in P c . Then the following holds. (A') Uniqueness of equilibria. If sign(det(M (x))) = (−1) s for all positiveequilibria x ∈ V ∩ P + c , then there is exactly one positive equilibrium in P c . Further, this equilibrium is non-degenerate.
then there are at least two positive equilibria in P c , at least one of which is non-degenerate. If all positive equilibria in P c are non-degenerate, then there are at least three and always an odd number.
The proof of Theorem 1 is based on relating det(M (x)) to the Brouwer degree of ϕ c at 0 (see §1- §4 in the Supporting Information). Note that the only situation that is not covered by Theorem 1 is when sign(det(M (x))) takes the value 0 for some x, but never the value (−1) s+1 . The determinant of M (x) is the same as the core determinant in [30,Lemma 3.7]. See also [10,Remark 9.27].
To check whether the sign conditions in part (A') or (B') hold requires information about the equilibria in P + c . As such, these conditions are difficult to check. If sign(det(M (x))) is constant for all x in a set containing the positive equilibria, then the condition in (A') is always fulfilled. In particular, this is the case for injective networks, where sign(det(M (x))) = (−1) s for all x ∈ R n >0 [10] (see also [31][32][33][34][35][36] for related work on injective networks). The latter might be verified or falsified without any knowledge about the equilibria of the system (see the comments to Step 5 and Step 7 in the section "Procedure for finding parameter regions for mono-and multistationarity").
Corollary 1 (Unique equilibria). Assume that the assumptions of Theorem 1 hold and that sign(det(M (x))) = (−1) s for all x ∈ R n >0 . Then there is exactly one positive equilibrium in each stoichiometric compatibility class. Further, this equilibrium is non-degenerate.
The conclusions of Theorem 1 refer specifically to non-degenerate equilibria. Non-degenerate equilibria are always isolated from each other within a given stoichiometric compatibility class, as det(M (x)) = 0 ensures M (x) is locally invertible. In some situations we might be able to "lift" non-degenerate equilibria of a reaction network to another reaction network that in some sense is larger, thereby proving lower bounds on the number of non-degenerate equilibria of the larger reaction network. This is for example the case if the smaller network is embedded in the larger [37,38], if the smaller network is without inflows/outflows while the larger has all inflows/outflows [39], or if the smaller is obtained by elimination of intermediate species [40]. Conditions for the existence of degenerate equilibria, where det(M (x)) is expected to change sign, are also known [20,41].
Positive parameterizations and a corollary. Verifying condition (A') or (B') is considerably easier if there exists a positive parameterization of the set V ∩ R n >0 of all positive equilibria. In this subsection we define such a parameterization and restate Theorem 1 as Corollary 2 in this situation. In the following sections this corollary will become the foundation for the procedure to partition the parameter space into regions with different equilibrium properties.
By a positive parameterization of the set of positive equilibria we mean a surjective function for some m < n, such thatx ∈ R m >0 is the vector of free variables. In other words, a positive parameterization implies that x 1 , . . . , x n are expressed at equilibrium as functions ofx: such that x 1 , . . . , x n are positive providedx is positive. Thus Typically, the number of free variables equals the corank of the network, that is m = d = n − s. We say that a parameterization is algebraic if the components Φ i (x) are polynomials or rational functions (quotients of polynomials) and can be given such that the denominator is positive for allx. See Fig 2 (Step 6) for an application to the running example. Note that the parameterizations considered here do not make use of the conservation relations.
A positive equilibrium Φ(x),x ∈ R m >0 , belongs to the stoichiometric compatibility class P c where c := W Φ(x).
Combining Eq (8) and Eq (9), it follows that the positive solutions to Eq (5) for a given c are in one-to-one correspondence with the positive solutions to Eq (9), that is, In order to restate Theorem 1 using the parameterization Φ, we consider the determinant of M (x) evaluated at Φ(x), Corollary 2 (Positive parameterization). Assume the reaction rate functions fulfil Eq (2) and let s = rank(N ). Further, assume that (i) The network is dissipative.
(ii) There are no boundary equilibria in P c , for all c ∈ R d such that P + c = ∅.
(iii) The set of positive equilibria admits a positive parameterization as in Eq (27).
Then the following holds.
(A) Uniqueness of equilibria. If sign(a(x)) = (−1) s for allx ∈ R m >0 , then there is exactly one positive equilibrium in each P c with P + c = ∅. Further, this equilibrium is non-degenerate.
then there are at least two positive equilibria in the stoichiometric compatibility class P c where c := W Φ(x). Further, at least one of the equilibria is non-degenerate. If all positive equilibria in P c are non-degenerate, then there are at least three equilibria and always an odd number.
Note that, contrary to Theorem 1, the stoichiometric compatibility class P c is not fixed in the corollary.
In the next section we formulate a procedure based on Corollary 1 and Corollary 2 to find regions of mono-and multistationarity. Before that we end this section with an application to the running example. The analysis is divided into seven steps which prelude the steps of the procedure.
Application of Corollary 1 and Corollary 2 to the running example.
We start with the setup given in Fig 1 and first check whether the sign condition of Corollary 1 is fulfilled, in which case there is a single equilibrium in all stoichiometric compatibility classes. The steps of the analysis are illustrated in Fig 2. The assumptions of the corollary are easily verified in this case. As we are assuming massaction kinetics, Eq (2) is fulfilled (Fig 2, Step 1). Further the network is conservative, hence dissipative and (i) is fulfilled (Fig 2, Step 2). It is easily seen that there are no boundary equilibria in any stoichiometric compatibility class with non-empty interior (Fig 2, Step 3). Hence (ii) is fulfilled. We then construct ϕ c (x) and calculate the determinant of M (x). It is a polynomial (in fact, a linear function) in x 1 , x 2 with coefficients containing both positive and negative terms (Fig 2, Step 4). By choosing (x 1 , x 2 ) ∈ R 2 >0 with x 1 large enough, the determinant of M (x) is positive (Fig 2, Step 5). Therefore Corollary 1 cannot be applied as s = 1. We note that this conclusion is independent of the specific choice of the parameter vector κ = (κ 1 , κ 2 , κ 3 ), so in fact it holds for all parameter values. Corollary 2 has the same assumptions as Corollary 1. We find a positive parameterization by solving the equilibrium equation for x 1 . That is, we treat it as an equation in x 1 , while x 2 (=x) is treated as a parameter. The function a(x 2 ) obtained by substituting x by Φ(x 2 ) in the determinant of M (x) is given in Fig 2, Step 6. It is clear from the expression of a(x 2 ) that it takes the sign −1 for all x 2 > 0. Also this conclusion does not depend on the specific value of κ. By application of Corollary 2(A) with s = 1, we conclude that there exists a unique positive non-degenerate equilibrium in each stoichiometric compatibility class with c > 0, for all values of the reaction rate constants (Fig 2, Step 7). The possibility of multiple equilibria is therefore excluded. In this particular example the existence of a positive parameterization is essential to draw the conclusion.
To illustrate how Corollary 2 can be used to find parameter regions for multistationarity, we consider the polynomial a( x) given in the first row of Fig 4, where n = 6 and s = 4 (the example is worked out in detail below): Input: N and v(x) from Fig. 1. We have s = 1.
Step 1: Kinetics is mass-action, thus (2) holds; f (x) and W are given in Fig 1. W is row reduced.
Step 2: The network is conservative and hence dissipative.
Step 4: Since i1 = 1 in W , we have The Jacobian M (x) of ϕc and its determinant are Step 5: One term of det(M (x)) has sign +1 = (−1) s+1 for all x ∈ R 2 >0 . We proceed to the next step.
Step 6 establishes a positive parameterization and finds the polynomial a(x).
Step 7 is similar to step 5, but for a(x).
Only one of the coefficients of the polynomial a(x) in x 4 , x 5 can be negative. If κ 3 ≤ κ 1 , then sign(a(x)) = (−1) 4 = 1 for all positivex. Corollary 2(A) implies that there is a unique positive non-degenerate equilibrium in each stoichiometric compatibility class with non-empty positive part.
Oppositely, we show that for κ 3 > κ 1 , Corollary 2(B) applies. For that, let x 4 = T and x 5 = T . Then a(x) becomes a polynomial in T with negative leading coefficient of degree 3. For T large enough, a( x) is negative, and we conclude that there existsx such that sign(a(x)) = (−1) 5 = −1. Corollary 2(B) implies that there exists a stoichiometric compatibility class P c that contains at least two positive equilibria. In summary, the region of the parameter space for which multistationarity exists is completely characterized by the inequality κ 3 > κ 1 .
Step 5 and 7 are sign analyses of det(M (x)) and a(x), respectively. These are crucial steps and essential for determining parameter regions with mono-and multistationarity. In general, the sign of a polynomial might be studied by studying the signs of the coefficients of the monomials in the polynomial. If all coefficients have the same sign, then the polynomial is either positive or negative for all x ∈ R n >0 , respectively,x ∈ R m >0 , depending on the sign, and Corollary 1, respectively, Corollary 2(A) applies. If this is not the case, then Corollary 2(B) might be applicable if we can show that the polynomial has the sign (−1) s+1 for somex ∈ R m >0 . In the two examples discussed here, the signs of det(M (x)) and a(x) are straightforward to analyse. However, this is not always the case, see the section "Checking the steps of the procedure".

Procedure for finding parameter regions for multistationarity
In the previous subsection we applied Corollary 1 and Corollary 2 to the running example by going through a number of steps corresponding to the conditions of the statements and the calculation of the determinant. In this section we outline the steps formally. Afterwards we discuss the steps and how they can be verified either manually or algorithmically, that is, without user intervention. Finally we devise an algorithm to conclude uniqueness of equilibria or to find regions in the parameter space where multistationarity occurs. We conclude this section with some extra examples that follow the steps of the procedure.
We assume the reaction rate functions v(x) depend on some parameters κ. The reaction rate functions are further assumed to be polynomials (as for mass-action kinetics) or quotients of polynomials (as for Michaelis-Menten and Hill kinetics with integer exponents).
The input to the procedure is v(x) and N (the stoichiometric matrix) and the output is parameter regions for which the network admits multistationarity or uniqueness of equilibria.
Procedure (Identification of parameter regions for multistationarity) Input: N and v(x) depending on κ.
1. Find f (x), a row reduced matrix W of size d × n such that W N = 0, and check that v(x) vanishes in the absence of one of the reactant species, that is, check that it satisfies Eq (2). 2. Check that the network is dissipative. 3. Check for boundary equilibria in P c for P + c = ∅ and c ∈ R d . 4. Construct ϕ c (x), M (x) and compute det(M (x)). 5. Analyze the sign of det(M (x)). Find conditions on the parameters κ such that sign(det(M (x))) = (−1) s for all x ∈ R n >0 , in which case Corollary 1 holds. If Corollary 1 does not hold for all κ, continue to the next step. 6. Obtain an algebraic parameterization Φ(x) of the set of positive equilibria for all κ, as in Eq (27), such that the coefficients of the numerator and the denominator of each Φ i (x) possibly depend on κ. Compute a(x) = det(M (Φ(x))). By hypothesis, a(x) can be written as the quotient of two polynomials inx with coefficients depending on κ, whose denominator takes positive values. 7. Analyze the sign of the numerator of a(x).
7a. Identify coefficients with sign (−1) s+1 and coefficients that can have different signs depending on the parameters.
7b. Use the terms corresponding to identified coefficients to construct parameter inequalities such that, whenever these inequalities hold, one has either sign(a(x)) = (−1) s for allx ∈ R m >0 or sign(a(x)) = (−1) s+1 for at least onex ∈ R m >0 , in which case either Corollary 2(A) or (B) holds. There is no guarantee that all steps of the procedure can be carried out successfully, let alone automatically. While step 1 and 4 usually are straightforward (only computational issues might arise for large networks), step 2, 3, 5, 6 and 7 might in particular require case specific approaches. However, there exist computationally feasible sufficient criteria that guarantee the conditions in each step can be checked efficiently.
Checking the steps of the procedure.
If the network is not dissipative, then at least one concentration grows to infinity over time. This is typically not the case for realistic networks, but it needs to be ruled out in order to apply the procedure.
We start by checking whether the network is conservative. This implies solving the linear system ω t N = 0 with the constraint ω > 0. Alternatively, conservation relations are often easily established by inspection of the reactions. For example, in many signalling networks, the total concentration of enzyme (free and bounded) and of substrate (phosphoforms) are conserved.
If the network is not conservative, then we check whether it is strongly endotactic [42,43]. Strongly endotactic reaction networks are in particular permanent, that is, dissipative and the compact set can be chosen such that it does not intersect the boundary of R n >0 , see [42][43][44] for details.
If the network is neither conservative nor strongly endotactic, then we can use the following proposition to decide on dissipativity (see the §3.2 in the Supporting Information).
Proposition 1 (Dissipative network). Let || · || be a norm in R n . Assume that for each c with P + c = ∅, there exists a vector ω c ∈ R n >0 and a number R > 0 such that ω c · f (x) < 0 for all x ∈ P c with ||x|| > R. Then the network is dissipative.
Thus, we look for vectors ω c with all coordinates positive and such that ω c ·f (x) < 0 for large x. To avoid restricting the parameter values, this computation should be done symbolically.
Step 3: absence of boundary equilibria.
For systems of moderate size it is often possible to establish nonexistence of boundary equilibria by arguments similar to those employed in the analysis of the running example: for each i, assume x i = 0, and show that it leads to a contradiction.
A systematic procedure to check for the existence of boundary equilibria relies on computing the so-called minimal siphons of the network [45]. A siphon is a set of species Z ⊆ {X 1 , . . . , X n } fulfilling the following closure property: if X i ∈ Z and X i is produced in reaction R j (that is, β ij > 0), then there exists X k ∈ Z such that X k is consumed in the same reaction (that is, α kj > 0). A minimal siphon is a siphon that does not properly contain any other siphon.
Proposition 2 (Siphons) ( [46,47]) If for every minimal siphon Z there exists a subset {X i 1 , . . . , X i k } ⊆ Z, and a conservation relation λ 1 x i 1 + · · · + λ k x i k = c for some positive λ 1 , . . . , λ k , then the network has no boundary equilibria in any stoichiometric compatibility class P c with P + c = ∅. The hypothesis of the proposition can be summarised by saying that each minimal siphon contains the support of a positive conservation relation.
For example, the running example has only one minimal siphon, namely {X 1 , X 2 }. The conservation relation x 1 + x 2 = c fulfils the requirement of Proposition 2, and hence the network has no boundary equilibria in any P c with c > 0.
More information about using siphons to preclude boundary equilibria is given in the section "Computational issues" below and in §5.1 of the Supporting Information.
Step 5: determining the sign of det(M (x)). If the kinetics is mass-action, then det(M (x)) is a polynomial in x. In general, if the reaction rate functions are rational functions in x, then so is det(M (x)). In the latter case, if the jth reaction rate function fulfils We determine conditions on the parameters such that all coefficients of p(x) have sign (−1) s . Then the sign of det(M (x)) is also (−1) s for all x ∈ R n >0 and Corollary 1 holds.
Step 6: finding an algebraic positive parameterization.
Computer algebra systems like Maple or Mathematica can be used to find a parameterization. One strategy is to solve the equations f i (x) = 0, i / ∈ {i 1 , . . . , i d }, for some subset of (at most) s variables, treating the remaining (at least) d variables as coefficients of the system. If a parameterization found in this way exists but is not positive, another set of variables should be tried out. This can be systematically addressed by trying out all possible subsets of variables. It requires computation and analysis of at most ( n d ) parameterizations. Alternatively, one can compute the circuits of degree one of the matroid associated with the equilibrium equations [48].
In some cases, the network structure implies that a positive parameterization of the set of equilibria exists. A set, say {X k+1 , . . . , X n } with n − k elements for some k, is non-interacting if two species never appear on the same side of a reaction and they have coefficient at most one in all reactions. In this case the equilibrium equations f k+1 (x) = · · · = f n (x) = 0 form a linear system in the variables {x k+1 , . . . , x n }. Provided that the determinant of the coefficient matrix of the linear system is not identically zero, this system can be solved and we obtain a positive parameterization of the non-interacting variables x k+1 , . . . , x n at equilibrium in terms of the remaining variables x 1 , . . . , x k [49,50]. A necessary condition for the determinant of the coefficient matrix not being identically zero is that there is no conservation relation of the form x i 1 + · · · + x i l with i 1 , . . . , i l ∈ {k + 1, . . . , n}. If a non-interacting set with k = d exists, that is, with s = n − d elements, then this guarantees the existence of the desired parameterization. In the running example there is not a non-interacting set because both species have coefficient 2 in the reaction 2X 1 → 2X 2 .
The non-interacting condition can be relaxed in some cases by requiring that none of the species in {X k+1 , . . . , X n } appear together in a reactant (these sets are called reactant-noninteracting [51]). Proceeding as above, provided that the determinant of the coefficient matrix is not identically zero, x k+1 , . . . , x n can be expressed at equilibrium in terms of x 1 , . . . , x k . Conditions that ensure this is a positive parameterization are given in [51]. In the running example, species X 1 is a reactant-non-interacting set and we can obtain a positive parameterization of x 1 in terms of x 2 , see Fig 1 and Fig 2. If the network admits so-called toric steady states, then a positive parameterization also exists [11].
Step 7: the sign of a(x) and the Newton polytope. This is perhaps the hardest step of all. We write a(x) = p(x)/q(x) with q(x) positive for allx and would like to determine the sign of p(x). We first look for conditions that ensure uniqueness of positive equilibria by imposing that all coefficients of p(x) as a polynomial inx have sign (−1) s . We next identify the monomials of p(x), where the sign of the coefficient, say β, is (−1) s+1 for some parameter values. For each of these monomials we check whether the monomial can "dominate" the sign of p(x). That is to say, if sign(β) = (−1) s+1 , then we determine whether there is anx such that p(x) also has the sign (−1) s+1 . If it is the case, then the condition sign(β) = (−1) s+1 is a sufficient condition for multiple equilibria according to Corollary 2(B).
Given a coefficient of a monomial with sign (−1) s+1 , it might not be straightforward to decide if the polynomial p(x) has the same sign for some value ofx. (For example, the polynomial x 2 −2xy +y 2 = (x−y) 2 has one monomial with negative sign, but the polynomial itself can never be negative.) When the number of variables is small, one can attempt to decide the sign as we did in the examples above. Otherwise, our strategy is to determine whether the monomial of interest corresponds to a vertex of the Newton polytope. If that is the case, then the monomial can dominate the sign of p(x) (see §5.2 in the Supporting Information). The Newton polytope of p(x) is defined as the convex hull of the exponent vectors α = (α 1 , . . . , α m ) ∈ R m corresponding to the monomialsx α 1 1 · . . . ·x αm m of p(x). If α is a vertex of the Newton polytope, then there existsx ∈ R m >0 such that the sign of p(x) agrees with the sign of the coefficient of the monomial (see §5.2 in the Supporting Information).

An algorithm.
In the previous subsection we have outlined computational criteria that might be used to verify the conditions of the steps in the procedure. These computational criteria are only sufficient, that is, even if they fail the procedure might still work on the given network. For example, a sufficient computational criterion for the absence of boundary equilibria is based on Proposition 2. However, it might happen that Proposition 2 cannot be applied, but that the network nonetheless has no boundary equilibria in stoichiometric compatibility classes with non-empty interior.
We have collected sufficient computational criteria that guarantee the conditions of the procedure are fulfilled. In this way the procedure is formulated as an algorithm with decision diagram shown in Fig 3. If one step of the algorithm fails, then we say that the algorithm ends indecisively. In that case we might check whether the step can be verified by other means.
For simplicity, we have restricted to mass-action kinetics. Under this assumption, det(M (x)) is a polynomial in x and the parameters κ, and a(x) is a rational function inx and κ because the parameterization is assumed to be algebraic.

Applications to selected examples
To illustrate several aspects of the algorithm we provide a detailed step-by-step analysis of a collection of examples.

Two-component system
We have chosen this example to illustrate the situation where an algebraic parameterization is not required, as already det(M (x)) is of constant sign. The algorithm therefore stops successfully at Step 5 (and consequently skips Step 6 and 7). Step 1 Step 2 Step 3 Step 4 Step 5 Step 6 Step 7 If that is the case, the corresponding condition might still be verified manually and the algorithm resumed from the next following step.

Yes
We consider a simple version of a two-component system consisting of a histidine kinase HK that autophosphorylates and transfers the phosphate group to a response regulator RR, which undergoes autodephosphorylation. The reactions of the network are We let X 1 = HK, X 2 = HK p , X 3 = RR and X 4 = RR p . The stoichiometric matrix N and a row reduced matrix W such that W N = 0 are The matrix W gives rise to the conservation relations x 1 + x 2 = c 1 and x 3 + x 4 = c 2 . With mass-action kinetics, the vector of reaction rates is v(x) = (κ 1 x 1 , κ 2 x 2 x 3 , κ 3 x 4 ), and the function We apply the algorithm to this network.
Step 1. Mass-action kinetics fulfills assumption in Eq (2) on the vanishing of reaction rate functions. The function f (x) and W are given above. The matrix W is row reduced.
Step 3. The minimal siphons of the network are {X 1 , X 2 } and {X 3 , X 4 }. These two sets are the supports of the conservation relations. By Proposition 2, there are no boundary equilibria in any P c as long as P + c = ∅.
Step 4. With our choice of W , we have i 1 = 1, i 2 = 3. Hence ϕ c is obtained by replacing the components f 1 (x), f 3 (x) of f (x) by the expressions derived from the two conservation relations: The Jacobian M (x) of ϕ c and its determinant are Step 5. All terms of det(M (x)) have sign +1 = (−1) s , since s = 2, and thus the conclusion of Corollary 1 holds. The network admits exactly one non-degenerate equilibrium point in every stoichiometric compatibility class with non-empty positive part.

Hybrid histidine kinase
This example has been analysed in [52]. Taken with mass-action kinetics the network is known to be multistationary for specific choices of reaction rate constants. We have chosen this example to illustrate how the algorithm can be used to sharpen known results: not only does it establish multistationarity for some parameter values, it provides precise conditions for when it occurs and allows a complete partition of the parameter space into regions with and without multistationarity. It also illustrates the use of an algebraic parameterization, which can be obtained by identifying sets of reactant-non-interacting species, and the use of the Newton polytope in Step 7. This reaction network is an extension of the two-component system discussed above and it is given in the first row of Fig 4. Specifically, the histidine kinase is assumed to be hybrid, that is, it has two ordered phosphorylation sites [52]. Whenever the second phosphorylation site is occupied, the phosphate group can be transferred to a response protein.
Using the notation X 1 = HK 00 , X 2 = HK p0 , X 3 = HK 0p , X 4 = HK pp , X 5 = RR and X 6 = RR p , the stoichiometric matrix N and a row reduced matrix W such that W N = 0 are The matrix W gives rise to the conservation relations We assume mass-action kinetics multiple equilibria: unique equilibria: and the function is We apply the algorithm to this network.
Step 1. Mass-action kinetics fulfills assumption in Eq (2) on the vanishing of reaction rate functions. The function f (x) and W are given above. The matrix W is row reduced.
The network has two minimal siphons {X 1 , X 2 , X 3 , X 4 } and {X 5 , X 6 }, which are respectively the supports of the two conservation relations. We apply Proposition 2 to conclude that there are no boundary equilibria in any P c as long as P + c = ∅.
Since i 1 = 1, i 2 = 5, the function ϕ c is obtained by replacing the components f 1 (x), f 5 (x) of f (x) by the expressions derived from the two conservation relations: The Jacobian M (x) of ϕ c (x) and its determinant are Step 5. The sign of the first coefficient of det(M (x)) depends on the parameters. If κ 1 ≥ κ 3 , then the sign is positive and det(M (x)) has sign +1 = (−1) 4 (s = 4) as the remaining terms are positive. According to Corollary 1, there is a single non-degenerate equilibrium in each stoichiometric compatibility class with non-empty positive part. If κ 1 < κ 3 , then Corollary 1 cannot be applied. We proceed to the next step to investigate the parameter space further.
Step 6. The set {X 1 , X 2 , X 3 , X 6 } is reactant-non-interacting and consists of s = 4 elements. We solve the equilibrium equations f 1 = f 2 = f 3 = f 6 = 0 for x 1 , x 2 , x 3 , x 6 . This gives the following algebraic parameterization Φ : The function a(x), which is det(M (x)) evaluated at Φ(x 4 , x 5 ), is the polynomial given in the first row of Fig 4. Step 7. We assume κ 1 < κ 3 , as the case κ 1 ≥ κ 3 is analysed in Step 5. Only one coefficient of a(x) has sign −1 = (−1) s+1 = (−1) 5 . The monomial associated with this term is x 4 x 2 5 . As the point (1, 2) (the degrees of the monomial) is a vertex of the Newton polytope (see Fig 4), then there existsx ∈ R m >0 such that the sign of a(x) is −1. Corollary 2(B) implies that there exists c = (c 1 , c 2 ) such that P c contains at least two positive equilibria.
Multistationarity is thus completely characterized by the inequality κ 3 > κ 1 . This condition states that the reaction rate constant for phosphorylation of the first site of the hybrid kinase is larger if the second site is phosphorylated than if it is not.

Gene transcription network
We consider the gene transcription network given in row 2 of Fig 4. This example has been studied in [53]. The particularities of this example are that the network is dissipative but not conservative, and that it displays multistationarity for all parameters κ. Further, this network illustrates the situation where the algorithm stops inconclusively at some step, but can be resumed after successful manual verification.
The network represents a gene transcription motif with two proteins P 1 , P 2 , produced by their respective genes X 1 , X 2 , and such that P 2 dimerises [53]. Further, the proteins cross regulate each other as depicted in Fig 4. Using the notation X 1 = X 1 , X 2 = X 2 , X 3 = P 1 , X 4 = P 2 , X 5 = X 2 P 1 , X 6 = P 2 P 2 , and X 7 = X 1 P 2 P 2 the stoichiometric matrix N and a row reduced matrix W such that W N = 0 are From W we find the conservation relations x 1 + x 7 = c 1 and x 2 + x 5 = c 2 . Here s = 5. We consider mass-action kinetics such that We apply the algorithm to this network: Step 1. Mass-action kinetics fulfills assumption in Eq (2) on the vanishing of reaction rate functions. The function f (x) and W are given above. The matrix W is row reduced.
The network is neither conservative nor strongly endotactic. Thus the algorithm terminates inconclusive. We take a manual approach: we pick ω c = (1, 1, 1, 1, 2, 2, 3) ∈ R 7 >0 and observe Note that x 1 , x 2 are bounded (due to the conservation relations) while x 3 , x 4 can be arbitrarily large. Then, for x 3 , x 4 large enough, ω c ·f (x) < 0 and the network is dissipative by Proposition 1 (as has been shown in [53] by other means).
Step 3. This network has two minimal siphons: {X 1 , X 7 } and {X 2 , X 5 }, which are the supports of the two conservation relations. Therefore, by Proposition 2, there are no boundary steady states in stoichiometric compatibility class with non-empty positive part.
In section §5.1 in the Supporting Information we illustrate how to apply a simplification technique, based on the removal of so-called intermediates and catalysts, to check whether Proposition 2 holds for this network.
Step 5. One coefficient of det(M (x)) has sign (−1) s+1 = 1 for all values of κ. Thus we proceed to the next step.
Step 6. There is not a set of non-interacting species nor reactant-non-interacting with s = 5 elements. Thus the algorithm terminates inconclusively. We take a manual approach and solve the equilibrium equations f 3 = f 4 = f 5 = f 6 = f 7 = 0 for x 1 , x 2 , x 3 , x 6 , x 7 . This gives the following algebraic parameterization Φ : R 2 >0 → R 7 >0 of the set of equilibria in terms of x = (x 4 , x 5 ): Evaluating det(M (x)) at Φ(x 4 , x 5 ) we obtain the polynomial Step 7. The coefficient of the monomial x 2 4 x 5 of the numerator of a(x 4 , x 5 ) has sign (−1) s+1 = (−1) 6 = 1. Since the monomial x 2 4 x 5 is a vertex of the associated Newton polytope (see Fig 4), there exists (x 4 , x 5 ) ∈ R 2 >0 such that the sign of a(x 4 , x 5 ) is 1. We conclude from Corollary 2(B) that for all κ i > 0 there exists c such that P c contains at least two positive equilibria.

Special classes of networks
There are several classes of networks for which some of the steps of the procedure are automatically fulfilled. We review some of them here.
Post-Translational Modification (PTM) networks consist of enzymes (E i ), substrates (S i ) and intermediate species (Y i ) [25]. Allowed reactions are of the form All intermediates are assumed to be the reactant, respectively, the product of some reaction. These networks are conservative (hence dissipative) and boundary equilibria are precluded provided the underlying substrate network obtained by ignoring enzymes and intermediates is strongly connected [47], see also §5.1 in the Supporting Information. When equipped with mass-action kinetics, these networks have a non-interacting set with d elements consisting of all enzymes and some of the substrates, namely one per (minimal) conservation relation involving the substrates [25]. Thus, a positive parameterization can always be found under the conditions stated above in Step 6. The class of PTM networks is contained in the class of cascades of PTM networks. Also this class admits a positive parameterization in terms of the concentrations of the enzymes and some of the substrate forms [26]. Cascades of PTM networks might further be generalized to so-called MESSI networks [54]. These networks are all conservative. Easy-to-check conditions for the absence of boundary equilibria and to decide whether the network admits toric steady states (and hence a positive parameterization) are given in [54].
A class of networks that cannot have boundary equilibria in any stoichiometric compatibility class with non-empty interior is given in [55].
The two examples in Table 1 are both PTM networks. Hence they are conservative and positive parameterizations exist. The underlying substrate network is strongly connected (they pass the criterion based on minimal siphons). For both networks the conditions shown in Table  1 are obtained by the algorithm. See §6.1 and §6.2 in the Supporting Information. For illustration purposes, we apply the algorithm in §6.3 of the Supporting Information to an additional network and show that it is monostationary.

Computational issues
The computational complexity of some of the steps in the procedure are demanding. Some conditions can be checked using linear algebra and do not depend on parameter values, others depend on parameter values and require symbolic manipulations. In some situations, the calculation can be done for even large networks at the cost of time, while in other situations symbolic software (like Mathematica and Maple) have inherent limits to what it can process. We offer here a few remarks about computational strategies and time complexity.
(1) Dissipativity. There are efficient algorithms to check whether the network is conservative and strongly endotactic, using linear algebra or mixed-integer linear programming [18,44]. We are not aware of a systematic way to check if Proposition 1 is fulfilled or not.
(2) Finding the minimal siphons of a network requires in general exponential time and there might be exponentially many of these [56]. Different algorithms developed in Petri Net theory can be applied to find the minimal siphons; see for example [45,46,56] and references therein. The complexity of this computation can often be substantially reduced by removing so-called intermediates and catalysts from the network [47] (see §5.1 in the Supporting Information for details).
(3) Finding all non-interacting and reactant-non-interacting sets requires in general exponential time. One strategy is the following. We first remove all species S i for which α ij > 1 or β ij > 1 for some reaction R j (the latter constraint is omitted if we are looking for reactantnon-interacting sets only). Then we build non-interacting (reactant-non-interacting) sets by adding new species recursively until no more species can be added without having an interacting pair of species in the set.
(4) Calculation of the symbolic determinant of the matrix M (x), and hence also of a(x), often fails in our experience for networks with more than 20 variables on common laptops [57]. However, this clearly depends on the sparsity of the matrix M (x), that is, on the number and order of the reactions. Strategies to reduce the complexity of the computation by expanding the determinant along the non-symbolic rows (conservation relations) were inspected in [57]. (6) Finding the vertices of the Newton polytope can be done with existing symbolic software, for example Polymake [58] or Maple, as we demonstrate in the Supporting Information.
We stress that it is always beneficial to guide the procedure/algorithm whenever possible in the sense that, if something is known for the network, there is no reason to go through many possibilities.

Discussion
The main result of this paper, the procedure to identify parameter regions for unique and multiple equilibria, combines Brouwer degree theory and algebraic geometry. In particular, under the assumptions of Corollary 2, we show that there exist stoichiometric compatibility classes with at least two equilibria if, and only if, a certain multivariate polynomial can attain a specific sign.
Discriminating regions of the parameter space where multistationarity occurs is a hard mathematical problem, theoretically addressable by computationally expensive means [22]. Our approach beautifully overcomes these difficulties by building on a simple idea, the computation of the Brouwer degree of a function related to a dissipative network. Additionally, not only closed-form expressions in the parameters are obtained, but, as illustrated in examples, these expressions are often interpretable in biochemical terms, providing an explanation of why multistationarity occurs.
The procedure applies theoretically to any choice of algebraic reaction rate functions. However, in practice, the procedure works well with mass-action kinetics. For example, we have considered the two-site phosphorylation cycle depicted in the second row of Table 1, but now modelled with Michaelis-Menten kinetics instead of mass-action kinetics. This network is known to be multistationary [59], and the conditions to apply Corollary 1 and Corollary 2 are valid. However, a positive algebraic parameterization does not exist, and hence our approach cannot be used to find parameter conditions for multistationarity.
However, Corollary 1 might be used with rational reaction rate functions for monostationary networks. This is the case for example for the one-site phosphorylation cycle S −− −− S p with Michaelis-Menten kinetics [59]. This network has two species and rank one. The sign of det(M (x)) is −1 for all parameter values and all x ∈ R 2 >0 . By Corollary 1, the network admits exactly one positive equilibrium in every stoichiometric compatibility class P c with P + c = ∅ for all parameter values.
If a reaction network does not have any conservation relation, then the set of equilibria consists typically of a finite number of points. In this case an algebraic parameterization is an algebraic expression of the equilibria in terms of the parameters of the system. Since m = 0, then R m consists of a single point and it follows directly that there is a unique equilibrium. Such an expression rarely exists. Therefore the procedure applies mainly to reaction networks with conservation relations. In particular, this rules out reaction networks where each species is produced and degraded.
Several natural questions remain outside the reach of our procedure. Firstly one would like to determine the particular stoichiometric compatibility classes for which there are multiple equilibria. As stated in Corollary 2, if sign(a(x)) = (−1) s+1 , then c := W Φ(x) defines a stoichiometric compatibility class with multiple equilibria. However, this only establishes c indirectly througĥ x. In some situations, it might be possible to find a positive parametrization that uses some of the conservation relations (ideally, all but one) and the stoichiometric compatibility classes with multiple/single equilibria would be determined up to a single parameter.
Secondly, one could ask for parameter regions that differentiate between the precise number of equilibria (that is, 0, 1, 2, . . .). This question should be seen in conjunction with the previous question: in typical examples, when there are two equilibria in a particular stoichiometric compatibility class, then there exists another class for which there are three. Hence the number of equilibria cannot be separated from the stoichiometric compatibility classes.
A third question concerns the stability of the equilibria, which cannot be obtained from our procedure. It is, however, known that if the sign of the Jacobian evaluated at an equilibrium is (−1) s+1 , then it is unstable [31]. This is in particular the case for an equilibrium fulfilling the condition in Corollary 2(B).
We have shown that for some reaction networks our procedure can be formulated as an algorithm. We consider therefore our research a step in the direction of providing 'black box tools' to analyse complex dynamical systems. Such tools would easily find their use in systems and synthetic biology, where it is commonplace to consider (many) competing models. A particular problem is to exclude models that cannot explain observed qualitative features, such as multistationarity.
Supporting Information: Proof of mathematical statements and examples. In this document we first prove the claims of the main text. Next, we provide details on how to check the steps of the procedure. Finally, we give details of the examples in Table 1 and include an extra example which is a PTM network.

Identifying parameter regions for multistationarity Supplementary Information
Carsten Conradi, Elisenda Feliu, Maya Mincheva, Carsten Wiuf September 11, 2018 In this document we prove the claims of the main text.
Sections 1 to 4 focus on the proofs of the theorem and its corollaries in the main text. We start by introducing some preliminaries before recapitulating the main facts about Brouwer degree theory. Then we compute the Brouwer degree for a special class of functions (Theorem 2). We proceed to introduce the necessary background on reaction networks and to state and prove a key result regarding the Brouwer degree of a reaction network with a dissipative semiflow (Theorem 3). In Section 4 we use Theorem 3 to prove Theorem 1 of the main text. The first four sections of the document are self-contained and do not require parallel reading of the main text. For this reason some parts of the main text are repeated here for convenience.
Subsequently in Section 5, we provide details on how to check the steps in the procedure of the main text. In Section 6 we give details of the examples in the main text and apply the algorithm to an extra network that is monostationary. A set B is convex if the following holds: Let B ⊆ R n be a convex set. We say that v ∈ R n points inwards B at x ∈ bd(B) if x + v ∈ cl(B) for all > 0 small enough. In particular, v = 0 points inwards B at all x ∈ bd(B). If v points inwards B at x ∈ bd(B), then it also points inwards cl(B) at x ∈ bd(B). The vector v points outwards B at x ∈ bd(B), if it does not point inwards B at x ∈ bd(B).
We will use the following facts about convex sets.
Lemma 1. Let B ⊆ R n be a convex set. Then the following holds: be the half-closed line segment between x 1 and x 2 . Then [x 1 , x 2 ) ∈ B.

Functions
Given an open set B ⊆ R n , we let C k (B, R m ) denote the set of C k -functions from B to R m . If B is open and bounded, then we let C k (cl(B), R m ) denote the subset of C k (B, R m )-functions f whose j-th derivative d j f , j = 0, . . . , k, extends continuously to the boundary of B. Equivalently, d j f is uniformly continuous in B for j = 0, . . . , k, since cl(B) is compact.
For f ∈ C 1 (B, R n ) and x * ∈ B, we let J f (x * ) ∈ R n×n be the Jacobian of f evaluated at x * , that is, J f (x * ) is the matrix with (i, j)-entry ∂f i (x * )/∂x j . We say that y ∈ R n is a regular value for f if J f (x) is non-singular for all x ∈ B such that y = f (x). If this is not the case, then we say that y is a critical value for f . If B ⊆ R n is open and bounded, f ∈ C 1 (cl(B), R n ) and y is a regular value for f such that y / ∈ f (bd(B)), then the set {x ∈ B|f (x) = y}. 2 Brouwer degree and a theorem

Brouwer degree
We first recall basic facts about the Brouwer degree. We refer to Section 14.2 in [62] for background and fundamental properties of the Brouwer degree. See also the lecture notes by Vandervorst [61].
In this section we let B ⊆ R n be an open bounded set. We use the symbol deg(f, B, y) to denote the Brouwer degree (which is an integer number) of a function f ∈ C 0 (cl(B), R n ) with respect to (B, y), y ∈ R n \ f (bd(B)).
A main property of the Brouwer degree is that if y / ∈ f (cl(B)), then deg(f, B, y) = 0 (but not vice versa) and if deg(f, B, y) = 0, then there exists at least one x ∈ B such that y = f (x). In particular, the Brouwer degree can be used to study the number of solutions to the equation provided y / ∈ f (bd(B)) and f ∈ C 0 (cl(B), R n ). The Brouwer degree deg(f, B, y) is characterized by the following properties: where the sum over an empty set is defined to be zero.

The Brouwer degree for a special class of functions
In this section we use Theorem 1 and the homotopy invariance of the Brouwer degree (A3) to compute the Brouwer degree of certain functions. Specifically, we are concerned with C 1 - and matrices W ∈ R d×n of maximal rank d. A priori there is no relationship between f and W . Assume that W is row reduced and let i 1 , . . . , i d be the indices of the first non-zero coordinate of each row, i 1 < . . . < i d . Let c ∈ R d and define the C 1 -function We say that ϕ c is constructed from f and W . The dependence of ϕ c on f and W is omitted in the notation. We will make use of this construction with different choices of f and W . Define the positive closed and open level sets of W by It follows readily that the two set are convex. By reordering the columns of W , the vector (x 1 , . . . , x n ) and the coordinates of f simultaneously, if necessary, we can assume without loss of generality that {i 1 , . . . , i d } = {1, . . . , d}. In this case, W has the block form where W ∈ R d×s , s := n − d, and I d is the identity matrix of size d. The last s coordinates of the function ϕ c come from f . Assuming this reordering, let π : R n → R s be the projection onto the last s coordinates. Using (15), it follows that In particular, for x, y ∈ R n fulfilling W x = W y, we have that x = y if and only if π(x) = π(y).
If W f (x) = 0, then it follows from (17) that f (x) = 0 if and only if π(f (x)) = 0. Our first result concerns the Brouwer degree of ϕ c . The proof of the theorem is adapted from the proof of Lemma 2 in [63] in order to account for the reduction in dimension introduced by P c .
Theorem 2. Let f : R n ≥0 → R m be a C 1 -function and W ∈ R d×n a matrix of rank d. Let s := n − d, c ∈ R d , P c as in (14) and ϕ c as in (13). Let B c be an open, bounded and convex subset of R n >0 such that (ii) f (x) = 0 and W f (x) = 0 for x ∈ bd(B c ) ∩ P c .
(iii) for every x ∈ bd(B c ) ∩ P c , the vector f (x) points inwards B c at x.
Proof. Without loss of generality, we might assume that W has the block form in (15). Choose an arbitrary pointx ∈ B c ∩ P c , which exists by assumption (i), and consider the continuous function G : cl(B c ) → R n defined by where π is the projection map onto the last s coordinates of R n . By (15), the Jacobian of G has the block form Therefore, det(J G (x)) = (−1) s for all x. In particular, 0 is a regular value for G. Furthermore, if G(x) = 0, then x ∈ P c since W x = c and π(x) = π(x). Using (17), we conclude thatx = x. Sincex / ∈ bd(B c ), it follows that G does not vanish on the boundary. We apply Theorem 1 to compute the degree of G for 0: deg(G, B c , 0) = sign(det(J G (x))) = (−1) s .
Consider now the following homotopy between the functions ϕ c and G: Clearly, H is continuous. To apply (A3) to find the degree of ϕ c , we need to show that H(x, t) = 0 implies that W x = c and hence x ∈ P c . Thus, we need to show that For t = 1, (18) follows from (17) using that f (x) = 0 and W f (x) = 0 for x ∈ bd(B c ) ∩ P c by assumption (ii). For t = 0, we have already shown that G does not vanish on the boundary of B c . Assume now that for t ∈ (0, 1), (18) does not hold. That is, there exists Since t−1 t < 0,x ∈ B c and x ∈ bd(B c ), it follows from Lemma 1(iii) that f (x ) points outwards B c at x , contradicting assumption (iii). Therefore

Setting
Consider a chemical reaction network with species set {X 1 , . . . , X n } and reactions: where α ij , β ij are non-negative integers. The left hand side is called the reactant complex and the right hand side the product complex.
The ODE system associated with the chemical reaction network G (as described in the main text) takes the formẋ where N ∈ R n× is the stoichiometric matrix and v(x) is the vector of rate functions, which are assumed to be C 1 -functions (e.g. mass-action monomials). We say that the network has rank s if the rank of the stoichiometric matrix is s and define d = n − s to be the corank of the network. The stoichiometric compatibility classes are the convex sets P c defined in (14), where W is a matrix such that the rows form a basis of im(N ) ⊥ . By construction, a trajectory of (21) is confined to the stoichiometric compatibility class where its initial condition belongs to. The positive stoichiometric compatibility classes P + c are defined accordingly.
The positive solutions to the system of equations ϕ c (x) = 0 with ϕ c as in (13), are precisely the positive equilibria of the network in the stoichiometric compatibility class P c .
Let φ(x, t) denote the flow of the ODE system and let the semiflow of the ODE system be the restriction of the flow to t ≥ 0. It is assumed that the choice of rate functions v(x) is such that v j (x) = 0 if x i = 0 for some i with α ij > 0. (22) In particular, mass-action kinetics fulfil this condition. Under this assumption, the non-negative and the positive orthants, R n ≥0 and R n >0 , are forward invariant under the ODE system (21), cf. [28,Section 16]. That is, if x 0 ∈ R n ≥0 (resp. R n >0 ), then the solution to the ODE system (21) with initial condition x 0 is confined to R n ≥0 (resp. R n >0 ): Forward invariance implies that the semiflow φ(x, t) maps R n ≥0 to itself for any fixed t ≥ 0 for which the solution is defined.
Since the dynamics is confined to the stoichiometric compatibility classes, this implies that for a point x at the relative boundary of P c , the vector f (x ) points inwards P c . Further, both P c and P + c are also forward invariant sets. Recall that these are convex sets.

Conservative and dissipative networks
A network is conservative if and only if the stoichiometric compatibility classes P c are compact subsets of R n ≥0 [29].
Definition 2. Consider a network with associated ODE systemẋ = N v(x). The semiflow of the network is dissipative if for all c ∈ R d such that P + c = ∅, there exists a compact set K c ⊆ P c such that φ(x, t) ∈ K c for all x ∈ P c and t ≥ t(x), for some t(x) ≥ 0. That is, all trajectories in P c enter K c at some point.
The set K c is called attracting (and sometimes absorbing) [64]. Equivalently, the semiflow of a network is dissipative if all trajectories are eventually uniformly bounded, that is, there exists a constant k > 0 such that lim sup t→+∞ x i (t) ≤ k for all i = 1, . . . , n and all initial conditions in P c , provided that P + c = ∅ for c arbitrary. If the semiflow of the network is dissipative, then the unique solution to the ODE system (21) for a given initial condition is defined for all t ≥ 0, in which case the semiflow is said to be forward complete.
Lemma 2. Consider a network with a dissipative semiflow and let c ∈ R d such that P + c = ∅. Then the following holds: (i) An attracting set K c can be chosen such that K c ∩R n >0 is non-empty, that is, K c ⊆ bd(R n ≥0 ).
(ii) All ω-limit points in P c of the system are contained in K c . In particular, all positive equilibria in P c belong to K c .
(iii) There exists an attracting set K c such that K c ⊆ K c , K c is forward invariant and all ω-limit points outside the boundary of R n ≥0 are in the interior of K c (relatively to P c ).
Proof. (i) Consider an attracting set K c ⊆ P c and assume that K c ⊆ bd(R n ≥0 ). Since P + c = ∅, there exists a compact set K c ⊆ P c that includes K c and such that K c ∩ R n >0 is non-empty. This set is also an attracting set.
(ii) If it were not the case, there would exist an ω-limit point x ∈ P c \K c , a trajectory φ(x, t) and a sequence of time points t i such that lim i→∞ t i = ∞ and lim i→∞ φ(x, t i ) = x . As K c is closed, there exists an open ball B (x ) in R n such that B (x ) ∩ K c = ∅ and φ(x, t) ∈ B (x ) for arbitrary many time points. However, this contradicts that K c is an attracting set.
(iii) By (ii) and choosing K c potentially larger, all ω-limit points outside the boundary of R n ≥0 are in the interior of K c (relatively to P c ). The existence of an attracting set K c that includes K c and is forward-invariant is proven in the first part of the proof of Lemma 2 in [64]. In the notation of [64], K c = K + .
The semiflow of a conservative network is dissipative. Indeed, it is sufficient to take K c = P c , since P c is compact. If the network is not conservative, then the semiflow associated with the network might still be dissipative (see Example "Gene transcription network" in the main text). However, in general, it is not straightforward to show that. In some cases it is possible to prove dissipativity by constructing a suitable Lyapunov function. It is the idea underlying the proof of the next proposition.

Proposition 1.
Assume that for all c ∈ R d such that P + c = ∅, there exists a vector ω ∈ R n >0 and a real number r > 0 such that ω · f (x) < 0 for all x ∈ P c with ||x|| ≥ r, where || · || is any norm. (Note that ω and r might depend on c.) Then the semiflow of the network is dissipative.
Proof. Let c ∈ R d with P + c = ∅ and let ω be as given in the statement. Define The function V (x) satisfies V (0) = 0 and V (x) > 0 for all x ∈ R n ≥0 , different from 0. Further, for ||x|| ≥ r and x ∈ P c ,V (x) = ∇V · f (x) = ω · f (x) < 0 by assumption. Thus, V (x) is a strict Lyapunov function and V (φ(x, t)) is strictly decreasing along trajectories φ(x, t) in P c as long as ||φ(x, t)|| ≥ r. Choose R > 0 such that The set K c is compact by construction and forward invariant sinceV (x) < 0 for all ||x|| ≥ r. Further, all trajectories eventually enter K c within finite time, that is, K c is attracting. Indeed, if this were not the case, then there would exist x ∈ P c , x / ∈ K c (hence ||x|| > r) such that V (φ(x, t)) is decreasing for all t ≥ 0 in the interval of definition and bounded below by R. As a consequence, the trajectory is defined for all t ≥ 0 and ( * ) lim t→∞ V (φ(x, t)) = R ≥ R for some R . Hence φ(x, t) is in B := {x | V (x) ≤ R + } for large t (and any > 0). Since B is compact it follows that the semiflow φ(x, t) has at least one ω-limit point in B . By virtue of ( * ), all ω-limit points x of φ(x, t) must fulfil V (x ) = R . Further, the set of ω-limit points is forward invariant and since V (x ) = R it must be thaṫ V (x ) = 0. This contradicts the assumption thatV (x) < 0 for all x with ||x|| ≥ r. We conclude that there exists t(x) ≥ 0 such that φ(x, t) ∈ K c for all x ∈ P c and t ≥ t(x). Hence, the semiflow is dissipative.

Degree for dissipative semiflows
The main results to establish a characterization of regions of multistationarity (Theorem 4) are Theorem 1 and the theorem below. The proof of the theorem relies on Theorem 2 and ideas developed in [64].
Theorem 3. Consider a network of rank s with an associated ODE systemẋ = f (x) where f (x) = N v(x) as in (21). Assume (22) holds on the rate functions and let W ∈ R d×n , d = n − s, be a row reduced matrix such that the rows of W form a basis of im(N ) ⊥ . Let c ∈ R d such that P + c = ∅. Further, assume that: Step (B). The open set U1 ⊆ R n (in green) is chosen such that Kc ⊆ U1 ⊆ B. In the C 1 -partition of unit, the support of ψ1 is in U1 and that of ψ2 is in R n \ Kc.
• The semiflow of the network is dissipative, and that That is, there are no boundary equilibria in P c .
Then there exists an open bounded and convex set B c ⊆ R n >0 that contains all positive equilibria of the network in the stoichiometric compatibility class P c , and such that where ϕ c is defined in (13) from f and W .
Proof. The idea of the proof is to construct a function g defined on R n ≥0 and a set B c ⊆ R n >0 such that the conditions of Theorem 2 are fulfilled for g, W and B c . If we let ϕ g c be the function ϕ c in (13) constructed from the function g and W , this will imply that deg(ϕ g c , B c , 0) = (−1) s . Subsequently, we will use homotopy invariance to conclude that also deg(ϕ c , B c , 0) = (−1) s .
The function g will be defined as where φ(x, t) is the semiflow ofẋ = f (x), K c is a suitably chosen attracting set, T is the maximum entrance time into K c from a specific set, and ρ(x) is an auxiliary function with certain useful properties (see below). The proof is divided into four steps. In step (A) we define the set B c , choose K c and find basic properties of B c and K c . In step (B), we construct the function ρ. In step (C), we properly define g and show that g, B c and W have the required properties to apply Theorem 2. In step (D) we show that ϕ g c and ϕ c are homotopy equivalent and conclude the proof of the theorem using the homotopy invariance of the Brouwer degree.
(A) Let K c ⊆ P c be as in Definition 2, that is, a compact attracting set of all trajectories with initial condition in P c . According to Lemma 2, K c can be chosen such that K c is forward invariant, K c ∩ R n >0 = ∅, and all ω-limit points in P + c are interior points of K c (relatively to P c ). Let B ⊆ R n be an open, bounded and convex set containing K c , that is, Then B c is also open, bounded and convex. Since K c ⊆ R n ≥0 ∩ B, then B c contains all points in K c except those on the boundary K c ∩ bd(R n ≥0 ). Further, see Figure 5.
Since f (x) = 0 for all x ∈ bd(R n ≥0 ) ∩ P c by assumption and K c contains all zeros of f in P c , then B c contains all zeros of f in P c , that is (B) The function ρ : R n ≥0 → R n in the definition of g is defined such that it has the following properties: We first construct two other functions ρ and ψ 1 , and subsequently define ρ : R n ≥0 → R n as the product ρ = ρ ψ 1 . Let x ∈ K c ∩ R n >0 and define ρ : be an open set containing K c (which exists since B is open), see Figure 5. Consider the open cover of R n given by U 1 and U 2 = R n \K c , such that U 1 ∩U 2 = ∅ and U 1 ∪U 2 = R n . Choose a C 1 -partition of unit ψ 1 , ψ 2 : R n → [0, 1] associated with this open cover. This implies in particular that the support of ψ i is included in U i and ψ 1 (x) + ψ 2 (x) = 1 for all x.
Define ρ : R n ≥0 → R n by ρ(x) = ψ 1 (x) ρ(x), x ∈ R n ≥0 (note the restriction to R n ≥0 ). This function fulfils properties (i)-(iv) above. Property (i): Follows by definition of ρ(x) = ψ 1 (x)( x − x), ψ 1 (x) ≥ 0 and Lemma 1(iii), using that x ∈ B c and x ∈ bd(B c ). Property (ii): Since the support of ψ 1 is contained in (C) Let T be defined as the maximum of the entry times to K c from any x ∈ cl(B c ) ∩ P c . The time T is finite because cl(B c ) ∩ P c is compact and the semiflow is dissipative with respect to K c . Note that once a trajectory is in K c , it stays there since K c is forward invariant Redefine T to be any positive number if T = 0.
We define g : R n ≥0 → R n , g(x) := Observe that W g(x) = 0 for all x ∈ bd(B c ) ∩ P c , using property (iv) in step (B) and that φ(x, T ), x ∈ P c . By definition of T , φ(x, T ) ∈ cl(B c ) ∩ P c if x ∈ cl(B c ) ∩ P c and hence 1 T (φ(x, T ) − x) points inwards B c at x ∈ bd(B c ) ∩ P c by convexity of cl(B c ). Also T ρ(x) points inwards B c at x by property (i) in step (B). Hence, g(x) points inwards B c at x ∈ bd(B c ) ∩ P c by convexity again.
Therefore, the function g, together with B c and W , fulfil the conditions of Theorem 2. By letting ϕ g c be the function ϕ c in (13) constructed from g and W , we conclude that deg(ϕ g c , B c , 0) = (−1) s .

(D)
We define a homotopy between ϕ c and ϕ g c on cl(B c ) × [0, T ] by The function H(x, t) is continuous since φ(x, t) is differentiable and is the semiflow ofẋ = f (x). Note that H(x, 0) = ϕ c (x) and H(x, T ) = (W x − c, π(g(x))) = ϕ g c (x). Thus H(x, t) is a homotopy between ϕ c (x) and ϕ g c (x). We need to show that H(x, t) does not vanish on the boundary bd(B c ).
If H(x, 0) = ϕ c (x) = 0, then x ∈ P c is an equilibrium of the ODE system. Hence H(x, 0) does not vanish on bd(B c ) since B c contains all zeros of f in P c , see (24). Now let x ∈ bd(B c ) and assume that H(x , t) = 0 for some t ∈ (0, T ]. It follows that x ∈ P c , hence x ∈ bd(B c ) ∩ P c , and π(φ(x , t) − x ) = −t 2 π(ρ(x )). (25) Using (17) and property (iv) in step (B) we have that By construction of K c , all fixed points and periodic orbits are contained in K c . If ρ(x ) = 0, then (26) implies x ∈ K c ∩ bd(B c ) as x ∈ bd(B c ) by assumption. However, this contradicts property (iii) in step (B). Hence, it must be the case that ρ(x ) = 0.
Using that x ∈ bd(B c ) ∩ P c from (25) and ρ(x ) = 0, we conclude that x ∈ bd(R n ≥0 ) by property (ii) in step (B), since x ∈ R n >0 ∩ bd(B c ) ∩ P c . It follows that there exists i such that x i = 0 and we have Here we have used that ψ 1 (x ) = 0, since ρ(x ) = ψ 1 (x ) ρ(x ) = 0 and ψ 1 (x ) is a scalar. Now, by the inequality above, φ(x , t) does not belong to R n ≥0 . However, this contradicts the forward invariance of R n ≥0 with respect to the flow. Therefore, H does not vanish on bd(B c ) × [0, T ]. With this in place, homotopy invariance of the Brouwer degree implies that Remark 1. The statement and proof of the theorem focus exclusively on one stoichiometric compatibility class, that is, on a fixed value c ∈ R d . Therefore, if a semiflow admits an attracting set in one specific stoichiometric compatibility class (and not necessarily in all), then the theorem and computation of the Brouwer degree holds for this specific class.

Multistationarity in dissipative networks
In this section we prove the theorem and corollaries stated in the main text, which are consequences of Theorem 3 from the previous section. Consider the Jacobian of the map ϕ c (x). Because ϕ c (x) is independent of c, we denote the Jacobian by M (x). The i-th row of this matrix is given as where W i is the i-th row of W . That is, one can think of M (x) as being the matrix obtained from the Jacobian of f (x), with the i j -th row, j = 1, . . . , d, replaced by the j-th row of W .
An equilibrium x * is said to be non-degenerate if M (x * ) has rank n, that is, if det(M (x * )) = 0.
Theorem 4. Assume the reaction rate functions fulfil (22), let s = rank(N ) and let P c be a stoichiometric compatibility class such that P + c = ∅ (where c ∈ R d ). Further, assume that (i) The semiflow of the network is dissipative. (ii) There are no boundary equilibria in P c .
Then the following holds: then there is exactly one positive equilibrium in P c . Further, this equilibrium is non-degenerate.
(B') Multiple equilibria. If sign(det(M (x))) = (−1) s+1 for some equilibrium x ∈ V ∩ P + c , then there are at least two positive equilibria in P c , at least one of which is non-degenerate. If all positive equilibria in P c are non-degenerate, then there are at least three and always an odd number.
Proof. The hypotheses ensure that we can apply Theorem 3. Therefore choose an open bounded convex set B c ⊂ R n >0 that contains all positive equilibria of the network in the stoichiometric compatibility class P c and such that Let V c be the set of positive equilibria in the stoichiometric compatibility class P c . Note that (A') Since sign(det(M (x))) = (−1) s = 0 for all equilibria in P c , 0 is a regular value for ϕ c . We can therefore apply Theorem 1 and obtain where #V c is the cardinality of V c . We conclude that #V c = 1 and therefore that there exists a unique positive equilibrium in the stoichiometric compatibility class. Furthermore, since sign(det(M (x))) = 0 for all equilibria, the equilibrium is non-degenerate. implies that there must exist at least two other points x , x ∈ V c , such that sign(det(M (x ))) = sign(det(M (x )) = (−1) s , that is, there are at least three positive equilibria in P c , all of which are non-degenerate. In this case by Corollary 1, there is an odd number of equilibria and they are all non-degenerate.
Assume now that 0 is not a regular value for ϕ c . Then there must exist another positive equilibrium x in P c for which the Jacobian of ϕ c (x ) is singular. This implies that there are at least two positive equilibria in P c , x * and x , one of which is non-degenerate.
In typical applications we find an odd number of equilibria (≥ 3), all of which are nondegenerate. Observe that the hypothesis for Part (A) holds if the sign of det(M (x)) is (−1) s for all x in a set containing the positive equilibria. In particular, this is the case if det(M (x)) = (−1) s for all x ∈ R n >0 . We assume now that the positive solutions to the system f (x) = 0 (with f (x) as in (21)) admit a parameterization x = (x 1 , . . . ,x m ) → (Φ 1 (x), . . . , Φ n (x)), for some m < n. That is, we assume that we can express x 1 , . . . , x n at equilibrium as functions ofx 1 , . . . ,x m : Proof. Given c, note that Φ induces a bijection between the sets V c and S c : An element of S c corresponds to a positive equilibrium in the stoichiometric compatibility class P c .
(A) Consider a stoichiometric compatibility class P c defined by c such that P + c = ∅. Let x ∈ V c andx such that x = Φ(x). Then det(M (x)) = det(M (Φ(x))) = a( x).
By hypothesis sign(det(M (x))) = sign(a( x)) = (−1) s . Since this holds for all equilibria in V c , Theorem 4(A') gives that there is exactly one positive equilibrium in P c , which is non-degenerate.
(B) Let x be such that sign(a( x)) = (−1) s+1 and let c be defined as in the statement of the theorem. Then x = Φ( x) is a positive equilibrium in V c for which the sign of det(M (x)) is (−1) s+1 . Theorem 4(B') gives the desired conclusion.

Details on the steps of the procedure
In this section we expand further on how to check step 3 and 7 of the algorithm.

On siphons and boundary equilibria
A proof of Proposition 2 in the main text for mass-action kinetics can be found in [46], where strategies to find siphons are also detailed. The proof in [46] is however valid for general kinetics fulfilling assumption (22) (see [47,Prop. 2]). Different algorithms developed in Petri Net theory can be applied to find the siphons of a reaction network.
For large networks, the task of finding the siphons can be daunting. A way to reduce the complexity of the computation is by the removal of intermediate species and catalysts [47]. We explain the key aspects of this reduction method here. The method is used in the examples below.
The first reduction concerns removal of intermediates. Intermediates are species in the network that do not appear interacting with any other species, are produced in at least one reaction, and consumed in at least one reaction. For example the species ES 0 in the reaction network is an intermediate. Given a network, we obtain a reduced network by "removing" some intermediates, one at a time. This is done in the following way. Say we want to remove an intermediate Y from the network. We remove all reactions of the original network that involve Y and add a reaction y → y whenever y → Y → y with y = y belongs to the original network. Here y and y are the reactant complex of a reaction y → Y and product complex of a reaction Y → y , respectively.
To illustrate this, we consider the removal of the intermediate ES 0 in the network (29). The reactions of the reduced network are obtained by considering all length 2 paths of the original network that go through ES 0 . We have two such paths: By "collapsing" these paths we obtain the reactions Clearly the process could be repeated now by choosing other intermediates of the network (if any). In this way we can obtain reduced networks by removing several intermediates.
The second reduction concerns removal of catalysts. Catalysts are species that whenever they appear in a reaction, then they appear at both sides and with the same stoichiometric coefficient. For example, E in the reaction network (30) is a catalyst. Catalysts are actually defined in more generality in [47], but we restrict to this scenario to keep the discussion simple. Catalysts are removed from a network by literally removing them from the reactions where they appear. Removal of E in the reaction network (30) yields the reaction network This network has one minimal siphon, namely {S 0 , S 1 }, and s 0 + s 1 = c is a conservation relation. By Proposition 2 in the main text it does not admit boundary equilibria in stoichiometric compatibility classes with non-empty positive part. The next proposition allows us to conclude that the original network in (29) neither admits boundary equilibria in stoichiometric compatibility classes with non-empty positive part.
Proposition 2 (Theorems 1 and 2 in [47]). Let G be a network and G be a network obtained after iterative removal of intermediates or catalysts from G as described above. Each minimal siphon of G contains the support of a positive conservation relation if and only if this is the case for G .
In several cases, removal of intermediates and catalysts yields a so-called monomolecular network. That is, a network whose complexes agree with some species or the complex zero. For example, the network in (31) is monomolecular. In this case, checking the hypothesis of Proposition 2 in the main text is straightforward, in view of the next lemma.
Lemma 3 (Proposition 3 in [47]). Let G be a monomolecular network. Each minimal siphon of G contains the support of a positive conservation relation if and only if all connected components of G are strongly connected.
The network in (31) is clearly strongly connected. Thus, we do not need to find the siphons of the network to conclude that each of its minimal siphons contains the support of a positive conservation relation and thereby conclude that (29) does not admit boundary equilibria in stoichiometric compatibility classes with non-empty positive part.
Corollary 3. Let G be a network and G be a network obtained after iterative removal of intermediates or catalysts from G as described above. If G is a monomolecular network with all connected components strongly connected, then G has no boundary equilibria in any stoichiometric compatibility class P c such that P + c = ∅.
For example, consider the "gene transcription network" of the main text. The species X 5 and X 7 are intermediates. The reaction network obtained upon their removal is: For this network, X 6 is an intermediate and X 1 , X 2 are catalysts. Removal of these three species yields the reaction network This is a strongly connected monomolecular network. By Corollary 3, there are no boundary equilibria in any P c as long as P + c = ∅. We have reached the same conclusion as in the main text without the need of finding the minimal siphons of the network.

Newton polytope
We write a multivariate polynomial f (x) ∈ R[x 1 , . . . , x n ] as a sum of monomials: . . x αn n and c α ∈ R, for which only a finite number are non-zero. The Newton polytope of f (x), denoted by N (f ), is a closed convex set in R n , defined as the convex hull of the exponents α ∈ N n for which c α = 0 (See [60, Section 2] for a definition of convex hull). The set of vertices of N (f ) is a subset of the set of points α for which c α = 0.
The following is a well-known fact about the Newton polytope of a polynomial. The proof of the fact is constructive and provides an explicit way to find x in Corollary 2(B). Thus it offers a way to find stoichiometric compatibility classes (i.e. values of c) for which multistationarity exists.
Proposition 3. Let f (x) = α∈N n c α x α and let α be a vertex of N (f ). Then there exists x ∈ R n >0 such that sign(f (x )) = sign(c α ).
Proof. By hypothesis c α = 0. Since α is a vertex in a bounded convex polytope, there exists a separating hyperplane ω · x = T that intersects the polytope only in α and such that ω · x < T for any other point x of the polytope (see e.g. Definition 3.5 and Theorem 3.8 in [65]). In particular, ω · α < ω · α for all vertices α = α . For Now f (y) is a well defined function for t ∈ R >0 , which tends to +∞ for t tending to infinity and c α > 0 and to −∞ for c α < 0 (by assumption c α = 0). Hence, by letting t be large enough, the sign of f (y) agrees with the sign of c α .
Finding the vertices in practice.
In the examples below, we find the vertices of the Newton polytope of the polynomial of interest as follows. We use Maple (version 2015). We construct first the polytope using the command PolyhedralSet and subsequently use the command VerticesAndRays, from the package PolyhedralSets, to find the vertices.
6 Details on the examples in the main text 6

.1 Phosphorylation of two substrates
In this subsection we consider the network in the first row of Table 1 in the main text.
We consider a system in which two substrates can be either unphosphorylated, A, B or phosphorylated A p , B p . Phosphorylation of both substrates is catalyzed by the same kinase K and dephosphorylation of A p , B p is catalyzed by the same phosphatase F . That is, the system consists of two futile cycles sharing kinase and phosphatase.
The reactions of the system are: This network is a PTM network with substrates A, B, A p , B p , enzymes K, F and intermediates AK, BK, A p F, B p F . It was shown in [66] that this network with mass-action kinetics is multistationary. Here we find the necessary and sufficient condition on the reaction rate constants for having multistationarity in some stoichiometric compatibility class. We let The stoichiometric matrix N of the network and a row reduced matrix W whose rows from a basis of im(N ) ⊥ are The rank of N is s = 6. The matrix W gives rise to the conservation relations where c 1 , c 2 , c 3 , c 4 correspond to the total amounts of kinase, phosphatase, substrate A and substrate B, respectively. With mass-action kinetics, the vector of reaction rates is v(x) = (κ 1 x 1 x 3 , κ 2 x 7 , κ 3 x 7 , κ 4 x 2 x 4 , κ 5 x 9 , κ 6 x 9 , κ 7 x 1 x 5 , κ 8 x 8 , κ 9 x 8 , κ 10 x 2 x 6 , κ 11 x 10 , κ 12 x 10 ).
We apply the algorithm to this network with the matrix N and the vector v(x).
Mass-action kinetics fulfils assumption (22). The function f (x) and W are given above and the matrix W is row reduced.
Step 2. The network is a PTM network, hence it is conservative and thus dissipative.
Step 3. We apply the reduction technique from Section 5.1. The network has four intermediates AK, A p F, BK, B p F . After their elimination, we are left with the reaction network This network has two catalysts: K, F . Their elimination yields the reaction network (the socalled underlying substrate network in the main text) This is a monomolecular network with two strongly connected components. By Corollary 3, there are no boundary equilibria in any P c for which P + c = ∅.
Step 4. For our choice of W , we have i 1 = 1, i 2 = 2, i 3 = 3, i 4 = 5. The function ϕ c (x) is thus The determinant of M (x) is a large polynomial. We omit it here.
Note that κ 3 , κ 6 , κ 9 , κ 12 are the catalytic constants of phosphorylation/dephosphorylation of A and B (k c1 , k c2 , k c3 , k c4 in the main text), and are the inverses of the Michaelis-Menten constants of K and F for each substrate. Therefore, the necessary and sufficient condition for multistationarity can be written in terms of the catalytic constants and the Michaelis-Menten constants, This proves the condition for multiple and unique equilibria given in the first row of Table 1 in the main text, by letting In particular, we have that • If κ 3 κ 12 > κ 6 κ 9 , then we need κ 3 κ 12 k M 1 k M 4 < κ 6 κ 9 k M 2 k M 3 for multiple equilibria to occur. • If κ 3 κ 12 < κ 6 κ 9 , then we need κ 3 κ 12 k M 1 k M 4 > κ 6 κ 9 k M 2 k M 3 for multiple equilibria to occur.

Two-site phosphorylation system
In this subsection we consider the network in the second row of Table 1 in the main text. The conditions given here were also found in [23], the paper that lay the foundations of this algorithm. In this work we consider a direct route using the function ϕ c and avoiding changes of variables. We explain here how to find the conditions using the algorithm in the main text. We consider a system in which one substrate undergoes sequential and distributive phosphorylation by a kinase K and sequential and distributive dephosphorylation by a phosphatase F . The three phosphoforms of the substrate are A, A p , A pp . The reactions of the system are: This network is a PTM network with substrates A, A p , A pp , enzymes K, F and intermediates AK, A p K, A p F, A pp F . We let The stoichiometric matrix N of the network and a row reduced matrix W whose rows from a basis of im(N ) ⊥ are The rank of N is s = 6. The matrix W gives rise to the conservation relations where c 1 , c 2 , c 3 correspond to the total amounts of kinase, phosphatase and substrate A, respectively.
We apply the algorithm to this network with the matrix N and the vector v(x).
Mass-action kinetics fulfils assumption (22). The function f (x) and W are given above and the matrix W is row reduced.
Step 2. This network is dissipative since it is a PTM network.
Step 3. The network has four intermediates AK, A p K, A p F, A pp F . After their elimination, we are left with a reaction network with two catalysts: K, F . Their elimination yields the following underlying substrate network This is a monomolecular network with two strongly connected components. By Corollary 3, there are no boundary equilibria in any P c for which P + c = ∅.
The Jacobian matrix M (x) = J ϕc (x) is The determinant of M (x) is a large polynomial. We omit it here.
Step 5. The determinant of M (x) has terms of sign (−1) s+1 = −1. We postpone the discussion of the conditions on the reaction rate constants for which all terms have sign (−1) s to Step 7. We proceed to the next step.
Step 6. This network is a PTM network and has a non-interacting set with s = 6 species: {X 4 , X 5 , X 6 , X 7 , X 8 , X 9 } = {A p , A pp , AK, A p F, A p K, A pp F }.
The free variables of this parameterization are the concentrations of the two enzymes and one of the substrates. We substitute x 4 , . . . , x 9 with their expressions in the parameterization in det(M (x)) to find a( x). The function a( x) is a large rational function with positive denominator which we do not include here.
Step 7. The numerator of a( x), the polynomial p( x), determines therefore the sign of a( x).
The coefficients are polynomials in κ 1 , . . . , κ 10 . The polynomial has 15 terms, 9 of which are positive for all values of the reaction rate constants. The remaining 6 coefficients are polynomials in κ 1 , . . . , κ 10 that can either be positive or negative.
The vertex highlighted in bold corresponds to the monomial x 3 1 x 2 2 x 3 , whose sign depends on b 1 (κ). By Proposition 3, if b 1 (κ) < 0, then there exists x such that p( x) is negative. Corollary 2(B) (Corollary 2(B) in the main text) gives that there is a stoichiometric compatibility class that admits positive multiple equilibria. This proves the condition for multistationarity given in the second row of Table 1 in the main text.
The exponent vector of the monomial corresponding to the coefficient α(κ), (2, 2, 1), is not a vertex of the Newton polytope. In this case it is uncertain whether the condition α(κ) < 0 is sufficient for multistationarity.

Two-substrate enzyme catalysis
This section contains an additional example to illustrate the application of the algorithm to a monostationary network for which a parameterization is required to reach the conclusion.
We consider a mechanism in which an enzyme E binds two substrates, S 1 , S 2 , in an unordered manner in order to catalyze the reversible conversion to the product P . A variation of this system was considered in [67]. The reactions of the system are: We let X 1 = E, X 2 = S 1 , X 3 = ES 1 , X 4 = S 2 , X 5 = ES 2 , X 6 = ES 1 S 2 , X 7 = P.
The stoichiometric matrix N of the network and a row reduced matrix W whose rows from a basis of im(N ) ⊥ are The rank of N is s = 4. The matrix W gives rise to the conservation relations where c 1 , c 2 , c 3 , c 4 correspond to the total amounts of kinase, substrate S 1 and substrate S 2 , respectively. With mass-action kinetics, the vector of reaction rates is v(x) = (κ 1 x 1 x 2 , κ 2 x 3 , κ 3 x 1 x 4 , κ 4 x 5 , κ 5 x 4 x 3 , κ 6 x 6 , κ 7 x 6 , κ 8 x 2 x 5 , κ 9 x 6 , κ 10 x 1 x 7 ).
We apply the algorithm to this network with the matrix N and the vector v(x).
Mass-action kinetics fulfils assumption (22). The function f (x) and W are given above and the matrix W is row reduced.
Step 2. The network is conservative since the concentration of every species is in the support of a conservation relation with positive coefficients. Therefore the network is dissipative.
This network has only one intermediate ES 1 S 2 . Its removal yields the reaction network The conservation relations of this new network are (with the notation above): The minimal siphons of the network are These siphons contain the support of the conservation relations for c 1 , c 2 , c 3 respectively. Thus, by Proposition 2 in the main text and Proposition 2, the original network does not have boundary equilibria in any stoichiometric compatibility class that intersects the positive orthant.
The Jacobian matrix M (x) = J ϕc (x) is The determinant of M (x) is a large polynomial. We omit it here.
Step 5. The determinant of M (x) has terms of sign (−1) s+1 = −1. We postpone the discussion of the conditions on the reaction rate constants for which all terms have sign (−1) s to Step 7. We proceed to the next step.