Topological descriptors of the parameter region of multistationarity: Deciding upon connectivity

Switch-like responses arising from bistability have been linked to cell signaling processes and memory. Revealing the shape and properties of the set of parameters that lead to bistability is necessary to understand the underlying biological mechanisms, but is a complex mathematical problem. We present an efficient approach to address a basic topological property of the parameter region of multistationary, namely whether it is connected. The connectivity of this region can be interpreted in terms of the biological mechanisms underlying bistability and the switch-like patterns that the system can create. We provide an algorithm to assert that the parameter region of multistationarity is connected, targeting reaction networks with mass-action kinetics. We show that this is the case for numerous relevant cell signaling motifs, previously described to exhibit bistability. The method relies on linear programming and bypasses the expensive computational cost of direct and generic approaches to study parametric polynomial systems. This characteristic makes it suitable for mass-screening of reaction networks. Although the algorithm can only be used to certify connectivity, we illustrate that the ideas behind the algorithm can be adapted on a case-by-case basis to also decide that the region is not connected. In particular, we show that for a motif displaying a phosphorylation cycle with allosteric enzyme regulation, the region of multistationarity has two distinct connected components, corresponding to two different, but symmetric, biological mechanisms.


Introduction
Bistable switches are frequently observed and studied in living systems, and have been linked to cellular decision making and memory processes [1,2].These switches arise in different forms; one common form in parametric systems is that of hysteresis [3], that is, the system is monostable for small or large values of a parameter, and has two or more stable steady states for intermediate values.When the parameter changes slowly enough to allow the system to remain approximately at steady state, the resulting steady states depend on whether the parameter is increased or decreased.This is illustrated in Fig. 1(a), which displays a hypothetical system with three steady states.When the parameter increases from low to high, the steady state goes through a bifurcation at a critical parameter value τ max , after which the system discontinuously settles to another region of the output space.If the parameter value is decreased again, the system remains at the high steady state value, that is, it does not return to the low steady state value immediately.This will first happen if the parameter is decreased beyond another critical value τ min < τ max .This behavior confers the switch with robustness: after a change of level of steady state value takes place, small fluctuations in the parameter will not reverse the change.The larger the interval [τ min , τ max ] where the system has several steady states, the higher the degree of robustness of the change.More complicate switches can arise if, for example, the system has more than one steady state (is multistationary) in two intervals of the parameter, as illustrated in Fig. 1(b,c).Irreversible switches are obtained if the relevant interval is of the form (0, τ max ].Fig. 1(a) and Fig. 1(b,c) are qualitatively different in one important aspect: the parameter region where the system has multistationarity is connected in Fig. 1(a) and disconnected (has two disjoint pieces) in Fig. 1(b,c).This phenomenon appears also in higher dimensions, where instead of a curve of steady states, there is a steady state manifold with "bends", and instead of having one parameter being varied, a vector of parameters is changed along a curve.The shape of the parameter region of (a-c) Input-output curves for hypothetical systems.Input is thought to be a parameter of the system that is varied, and the output is the concentration of a species at steady state.Dashed lines correspond to unstable steady states, and solid lines to stable steady states.(a) displays a simple hysteresis switch; (b-c) show input-output curves for systems where bistability arises in two disjoint intervals of input.(d-e) For a system with two parameters, the system has more than one positive steady state in the orange regions, and one in the purple regions.In panel (d) the multistationarity region is connected, while in (e) it has two connected components.
multistationarity (or multistationarity region for short), and specifically the number of pathconnected components, modulates the type of switches that can arise.If the multistationarity region is path connected, as in Fig. 1(d), then any two parameter values in the region can be joined by a continuous path completely included in the region, typically giving rise to simple hysteresis switches (if the system has three steady states for parameters in the region).If, on the contrary, the region has two path-connected components, as in Fig. 1(e), then any path joining two parameter points in different regions, necessarily goes through parameter points where the system does not have several steady states, allowing for complex switches to arise.Mathematically, understanding the connectivity of a region is a basic topological property of a set, and the number of connected components is called the 0th Betti number of the set.Higher order Betti numbers describe the shape of the set in more detail, for instance, the first Betti number is the number of "holes" of the set.Tools from topological data analysis can infer the Betti numbers of a set from sample data points.By generating points in the multistationarity region, properties of the shape of the region have been explored for a specific dual phosphorylation system (Fig. 2(h)) in [4], where it has also been suggested that lack of connectivity may indicate that different biological mechanisms underlie multistationarity.
In this work, we address connectivity of the multistationarity region for polynomial systems describing the steady states of biochemical reaction networks.We achieve this by using exact symbolic tools and theoretical results relating the multistationarity region with the region where a polynomial attains negative values.Specifically, we work in the framework of chemical reaction network theory [5,6], where extensive work has been done to decide whether a network exhibits multistationarity, i.e. whether the multistationarity region is non-empty [7].More recent work focuses on understanding and finding the multistationarity region, but here progress is scarce and often restricted to special systems, e.g.[8][9][10][11][12][13][14].
Finding and studying the region where a polynomial system has more than one positive solution is a mathematical problem that belongs to the realm of semi-algebraic geometry and quantifier elimination [15,16].Although there are generic methods to address this question, these have high complexity, and fail for realistic networks, even of moderate size.To overcome these difficulties, methods targeting the specificities of reaction networks systems have been developed, and some partial results, for example describing the projection of the multistationarity region onto a subset of the parameters, have been developed [9][10][11].
Here we present a new algorithm to assert that the multistationarity region is connected without explicitly finding it, which bypasses the use of computationally expensive algorithms from semi-algebraic geometry and quantifier elimination.In fact, we reduce the problem to computing the determinant of a symbolic matrix, and finding a point in the feasible region of a system of linear inequalities.This makes the algorithm successful for networks of moderate size.
We apply our algorithm to numerous motifs in cell signaling, known to exhibit multistationarity, and conclude that these have connected multistationarity regions.These systems are shown in Fig. 2, where for all subfigures but (c) and (h), we confirm that the region is connected.For the system in Fig. 2(h), previously suggested to have a connected multistationarity region [4], our method is inconclusive.For system (c), modeling the enzymatic phosphorylation of a substrate S with allosteric regulation of the enzymes [17], the algorithm is also inconclusive.In this case, a detailed inspection of the system employing ideas similar to those behind our algorithm, allows us to assert that the region has exactly two path-connected components.These components are contained in the subset of parameters where κ 3 > κ 6 or κ 6 > κ 3 respectively.These two parameters are the catalytic constants of the phosphorylation and dephosphorylation processes, respectively.Therefore, if for example κ 6 increases from a small value to a value larger than κ 3 , then the multistationarity region will be crossed twice, and a non-simple hysteresis switch arises.
The fact that the remaining networks in Figure 2 have a connected multistationarity region, does not forbid complex switches, as a parameter path could still enter and exit the multistationarity region several times.However, in this scenario we can assert that for any pair of parameter values of multistationarity, there is a path connecting them and completely included in the multistationarity region.This implies that the conditions yielding to multistationarity vary continuously, and we cannot separate the multistationarity region into two disjoint sets, each corresponding to a distinct biological mechanism.
Our algorithm builds on two previous results.First, in [9], the authors associated with each reaction network a function, whose signs are closely related to the multistationarity region.Second, in [18], we developed a new criterion to determine the connectivity of a set described as the preimage of the negative real half-line by a polynomial map.We connect these two key ingredients in Theorem 2.2 to give a criterion for connectivity of the multistationarity region.Our results establish a stronger property, namely path-connectivity of the region, which in turn imply connectivity.
The criterion relies on computing the determinant of a matrix with symbolic entries, and this task might become unfeasible for large matrices with many variables.To bypass this, we show that the network can be reduced by removing some reactions and connectivity of the multistationarity region for the reduced network can be translated into the original network, see Theorem 2.4.This paper is organized as follows.We first establish the framework and background material, and present our algorithm for connectivity of the multistationarity region.We proceed to apply our algorithm to the networks in Fig. 2, and while doing so, we illustrate the strengths and limitations of the algorithm.To keep the exposition concise, we compile the proofs of the main theorems in Section 6 at the end.

2.1.
Theory.The results of this work are framed in the context of chemical reaction network theory, a formalism to study reaction networks that goes back to the 70s with the works of Horn, Reaction networks arising in cell signaling.The subindex 'p' indicates a phosphorylated site.When writing 'p' and 'pp' it is assumed the substrate has two phosphorylation sites, and phosphorylation/dephosphorylation is ordered.When writing '0p' for example, it means the substrate also has two sites numbered 1 and 2, and the second one is phosphorylated.All networks are known to be multistationary.For all networks but (c) and (h), the multistationarity region is path connected.For network (c), the multistationarity region has two path-connected components, while for network (h) our approach is inconclusive.
Jackson and Feinberg [5,6].To help the unfamiliar reader, and to fix the notation, we start with a brief introduction.This part ends with the main theoretical result to decide whether the set of parameters where multistationarity arises is path connected (and hence connected).To keep the exposition simple for non-experts, the proofs of the statements are given in Section 6 at the end.
Reaction networks.A reaction network (S, R) is a collection of reactions R = {R 1 , . . ., R r } between species in a set S = {X 1 , . . ., X n }.In our applications, species will be proteins, such as kinases and substrates.The reactions will encode events such as complex formation or posttranslational modifications.Formally, each reaction connects to linear combinations of species, that is, it has the form: where the coefficients a ij , b ij are non-negative integer numbers.The net production of each species when the reactions takes place is encoded in the stoichiometric matrix We do not consider reactions where the reactant and product are equal, hence N has no zero columns.
For illustration purposes, we consider a reaction network that is small enough to get a good feeling about the formal concepts but large enough to display the relevant features.To construct such a reaction network, the article [19] was particularly helpful.Realistic networks will be considered in Section 3 of this work.We refer to the following reaction network as the running example: Here, two species are related by three reactions, so the corresponding stoichiometric matrix has two rows and three columns, see Fig. 3.
Mathematical modeling offers us tools to get insights into the dynamics of the network and understand the temporal changes in the concentrations of the species of the network.By encoding the concentrations of X 1 , . . ., X n into the vector x = (x 1 , . . ., x n ) ∈ R n ≥0 , and under the assumption of mass-action kinetics, the evolution of the concentrations of the species over time is modeled by the ODE system: where Here κ = (κ 1 , . . ., κ r ) ∈ R r >0 is the vector of reaction rate constants, which are parameters of the system.Observe that each component of v κ (x) corresponds to one reaction, and it is obtained by considering the coefficients of the reactant of the reaction as exponents.The function v κ (x) and the associated ODE system of the running example are shown in Fig. 3.
The ODE system (2) is forward invariant on stoichiometric compatibility classes [20], that is, for any initial condition x 0 the dynamics takes place in the stoichiometric compatibility class of x 0 , which is the set (x 0 + S) ∩ R n ≥0 , where S denotes the vector space spanned by the columns of N .To work with stoichiometric compatibility classes it is more convenient to have equations for them.These are obtained by considering a full rank matrix W ∈ R (n−s)×n such that W N = 0, where s is the rank of the stoichiometric matrix N .Then for each parameter vector c ∈ R n−s , the associated stoichiometric compatibility class is the set This class is the set (x 0 + S) ∩ R n ≥0 for any initial condition x 0 satisfying W x 0 = c.Such a matrix W is called a matrix of conservation relations, any equation defining a class is a conservation relation, and the parameter vector c is called a vector of total concentrations.
For the running example, we have n = 2, s = 1, and hence W has one row, as given in Fig. 3.The figure depicts also the stoichiometric compatibility classes for c = 2, 3, 4, which are compact.In general, a reaction network is called conservative if each stoichiometric compatibility class is a compact subset of R n ≥0 , and this is the same as asking that there is a vector with all entries positive in the left-kernel of N [21].Under this assumption, all trajectories of the ODE system (2) are completely contained in a compact subset, and hence no component can go to infinity.For the purposes of this work, it will be enough to require a milder condition, namely that there is a compact set that all trajectories enter in finite time and do not leave again.In this case, the reaction network is said to be dissipative.See [9] for ways to verify that a network is dissipative.For simplicity, our applications are conservative networks.
We denote the set of non-negative steady states of the ODE system (2) by and call it the steady state variety.The positive steady state variety consists then of the points in V κ with all entries positive and is denoted by V κ,>0 .A steady state x ∈ V κ is a boundary steady state, if one of its coordinates equals zero.For our results, we will need that stoichiometric compatibility classes intersecting R n >0 do not contain boundary steady states.We call therefore a boundary steady state relevant if it belongs to a class P c with P c ∩ R n >0 = ∅.We say that a pair of parameters (κ, c) enables multistationarity if the intersection of the positive steady state variety V κ,>0 with the stoichiometric compatibility class P c contains more than one point.The set of parameters that enable multistationarity form the multistationarity region.We show in the right panel of Fig. 3 some stoichiometric compatibility classes, steady state varieties and their intersection, showing that multistationarity arises.
The running example: Stoichiometric matrix N , reaction rate function vκ(x) and a full rank matrix W such that W N = 0: Equation of the stoichiometric compatibility class Pc: Steady state varieties Vκ and stoichiometric compatibility classes Pc for different choices of parameters (κ, c): (κ, c) = ((18, 1, 10), 3) enables multistationarity, (κ, c) = ((18, 1, 10), 4) does not enable multistationarity.For κ = (1, 1, 1) there is no total concentration c such that ((1, 1, 1), c) enables multistationarity.We now recall the main theorem from [9] that is key to describe the multistationarity region.The theorem requires choosing a matrix of conservation relations W that is row reduced.Then, if i 1 < • • • < i n−s are the indices of the first non-zero coordinates of each row of W , we construct the matrix M κ (x) from the Jacobian of f κ (x) by replacing the i j th row by the jth row of W (for j = 1, . . ., n − s).
Our running example is dissipative, as it is conservative, and it has no relevant boundary steady states.We also find that This expression is negative for κ = (18, 1, 10) and x * ≈ (0.2448, 2.7552) ∈ V κ .Since W x * = 3, we can conclude using Theorem 2.1, that the intersection of V (18,1,10) and P 3 contains more than one point.This is exactly what Fig. 3 indicates.
Parametrizations.In Theorem 2.1, it is crucial to evaluate the determinant of M κ (x) at points in the incidence set: Therefore, we need to be able to describe the points in V in a useful way.This is done by considering parametrizations of V. Loosely speaking, a parametrization is a function whose image set consists precisely of the points in R n >0 × R r >0 that belong to V. Formally, we define a parametrization of V as a surjective analytic map In practice, D is the positive orthant of some R k and Φ is described by polynomials or quotients of polynomials such that their denominators do not vanish on D. Below, we discuss how to choose a parametrization and show that there is always at least one.
Using Theorem 2.1, one can show (see Lemma 6.4) that the multistationarity region is closely related to the preimage of the negative real half-line under the polynomial map g : V → R, (x, κ) → (−1) s det(M κ (x)), (7) which can be described using a parametrization Following [22], we call the function g • Φ a critical function, and observe that it depends on the choice of the parametrization.We are ready to present the main theoretical result of this work, namely a criterion for connectivity of the multistationarity region.The statement tells us that we can look at the number of path-connected components of (g • Φ) −1 (R <0 ), that is, of the set of values ξ where g • Φ is negative.The proof of the following theorem can be found in Section 6.
Theorem 2.2 (Deciding connectivity).Consider a reaction network that is dissipative and does not have relevant boundary steady states.Let g • Φ be a critical function as in (8) such that the closure of (g Then the number of path-connected components of the multistationarity region is at most the number of path-connected components of (g • Φ) −1 (R <0 ).

2.2.
Algorithm for checking path connectivity.Theorem 2.2 gives a theoretical criterion to decide upon connectivity, from which one can establish an algorithm for connectivity with the following steps: (Step 1) Check that the reaction network is dissipative and does not have relevant boundary steady states.(Step 2) Find a parametrization Φ of V and compute the critical function g • Φ.
The important point is that each of these steps can be addressed computationally, and hence the algorithm can be carried through without manual intervention, at least for networks of moderate size.We proceed to describe each of these steps in detail.
(Step 1) has already been described in detail in [9], as it consists of verifying that the conditions to apply Theorem 2.1 hold.Computable criteria that are sufficient to ensure that the properties hold are presented in [9].These are, however, not necessary and hence it might not always be possible to decide upon this step.
To verify dissipativity, the first attempt is to show that the reaction network is conservative by finding a row vector w ∈ R n >0 such that w N = 0 .This can be checked by solving the system of linear equalities: w N i = 0, for all i = 1, . . ., r and w j > 0 for all j = 1, . . ., n, (9) where N i denotes the ith column of N .We already noticed that the running example is conservative, by choosing A sufficient criterion to preclude the existence of relevant boundary steady states arises by using siphons, that is, subsets of species such that for all species in the set and all reactions producing them, there is a species in the reactant also in the set, see [23,Theorem 2], [25, Proposition 2], [26].In a nutshell, the criterion requires that for each minimal siphon it is possible to choose w ∈ R n ≥0 with wN = 0 and such that the positive entries of w correspond exactly to the species in the siphon.For more details, we refer to [9,23,25].Note that this criterion also relies on solving linear inequalities.Our running example has only one siphon, namely {X 1 , X 2 }.As the two entries of w in (10) are positive, the criterion holds and the network does not have relevant boundary steady states.
(Step 2) asks for the choice of a parametrization of V and the computation of the critical function.To find a parametrization systematically, we consider so-called convex parameters introduced by Clarke in [27].Since then, they have been applied to study reaction networks, for example to detect Hopf bifurcations and study bistability [28][29][30][31].
The idea behind convex parameters is the simple observation that the rate function v κ (x) has to be in the flux cone: It is easy to see that F is a convex polyhedral cone containing no lines.
Using software with packages for polyhedral sets (see Methods), one can compute a minimal collection of generators E 1 , . . ., E ∈ R r of F. These generators are often called extreme vectors of the cone.Their choice is unique up to multiplication by a positive number.Since the flux cone F does not contain lines, each of its elements can be written as a non-negative linear combination of extreme vectors [32, Corollary 18.5.2],that is for each v ∈ F there exists some λ = (λ 1 , . . ., λ where E ∈ R r× denotes the matrix with columns E 1 , . . ., E .We call E a matrix of extreme vectors.This gives rise to the following convex parametrization: where A 1 , . . ., A r denote the columns of the matrix A := a ij ∈ R n×r of the coefficients of the reactants of the reactions, h A j is short notation for h is the diagonal matrix with diagonal entries given by v, and 1/h is taken component-wise. In Corollary 6.2(a) in Section 6, we show that Ψ is surjective if E does not have a row where all the entries are equal to zero, and hence Ψ is a parametrization of V.This restriction is not relevant for our purposes: a zero row of E is equivalent to ker(N ) not having any positive vector, and hence there is no positive steady state of the ODE system (2), see Corollary 6.2(b).In particular, the reaction network cannot be multistationary.In the rest of the work, we assume that E does not have a zero row and when this holds, we say that the network is consistent [23] (consistent networks are called dynamically nontrivial in other works, e.g.[24]).
For the convex parametrization, the critical function g • Ψ can be represented in a direct way using the following observation.The Jacobian of f κ (x) evaluated at Ψ(h, λ) equals J(h, λ) := N diag(Eλ)A diag(h), (12) for each (h, λ) ∈ R n >0 × R >0 , see [28].We construct the matrix M (h, λ) from J(h, λ) as above: if W is row reduced and i 1 < • • • < i n−s are the indices of the first non-zero coordinates of each row, replace the i j throw of J(h, λ) by the jth row of W .Then, it holds From this equality, one computes g • Ψ directly using symbolic software.Since the entries of M (h, λ) are polynomials in (h, λ), so is g • Ψ.In the following, we call this polynomial the critical polynomial.
Let us find the critical polynomial g • Ψ for the running example.The flux cone F and its extreme vectors are displayed in Fig. 4. Now, all we have to do is to compute the matrix product in ( 12) and replace the first row by (1,1).After taking the determinant and multiplying by (−1) s = −1, we obtain the critical polynomial: The above discussion shows that we can always find a suitable parametrization and compute the critical polynomial.In some cases, other types of parametrizations arise by parametrizing each V κ,>0 separately.This is done by first trying to express some variables among x 1 , . . ., x n , κ 1 , . . ., κ r in terms of the others using the equations in (5).For finding these expressions, a computer algebra system such as SageMath [33] or Maple [34] can be useful.Once such an expression is found, one should check whether it gives a well-defined surjective analytic map, that is a parametrization.If all the components of the parametrization Φ are quotients of polynomials with positive denominators, then g • Φ is a quotient of polynomials too, and its denominator is positive.
Let us see how this works in practice for the running example.We see from the ODE system in Fig. 3 that positive steady states are characterized by This expression gives the parametrization Combining Φ with g from ( 7), we get the critical function: In general, there is no guarantee that such a parametrizations for each κ can be found.However, there are broad classes of reaction networks allowing such a parametrization, for example networks with toric steady states [35] and post-translational modification systems [36] to name a few.As we always can find a critical function using the convex parametrization, one Red dots correspond to positive exponents, the blue square corresponds to the only negative exponent.(b) The preimage of the negative real half-line under f .might wonder what the value of these other type of parametrizations is.The point is that with this type, the reaction rate constants are still present in the parametrization, and this is useful to get information about what parameter values yield to multistationarity.This was the theme of [9], and we will explore this advantage later in the application of our algorithm to the network with allosteric reciprocal enzyme regulation in Fig. 2(c).
Finally, we discuss how to address (Step 3), now that we know how to compute the critical polynomial/function.To check whether the preimage of the negative real half-line under a critical function is path connected is in general hard and depends strongly on the parametrization.As we discussed in (Step 2), critical functions are in practice polynomials or rational functions with positive denominator.In the latter case, we can restrict to the numerator of the rational function.To verify the conditions in Theorem 2.2, it is then enough to study the preimage of the negative real half-line under a polynomial function restricted to the positive orthant.
Recall that a polynomial function can be written as and σ(f ) ⊆ N k is a finite set, called the support of f .To determine whether the preimage of the negative real half-line is path connected, one can use methods from real algebraic geometry [15,Remark 11.19], [37,Section 3].These methods work well for polynomials in few variables, but they scale poorly.If the polynomial has many variables, the computation is unfeasible.
A hyperplane in R k is the set of solutions µ ∈ R k of a linear equation where v ∈ R k \ {0}, a ∈ R and v • µ denotes the Euclidean scalar product of two vectors.Each hyperplane has two sides, which are described by the linear inequalities A hyperplane is called strictly separating if the positive and negative exponents of f are on different sides of the hyperplane and not all the negative exponents are on this hyperplane.For a geometric interpretation, we refer to Fig. 5.
Strict separating hyperplanes can be used to decide upon path connectivity of the multistationarity region using Theorem 2.2 via the following theorem.
Theorem 2.3 (Preimage of negative real half-line).[18, Theorem 3.9] Let f : R k >0 → R be a polynomial function.If there exists a strict separating hyperplane of the support of f , then f −1 (R <0 ) is path connected and its closure equals f −1 (R ≤0 ).
For the running example, the supports of the two critical functions ( 14) and ( 15) form a quadrilateral.In both cases, there is only one negative exponent, which is at a corner (vertex) of the quadrilateral.Hence, one can easily find strict separating hyperplanes for the supports of each polynomial.To get a geometric intuition, we investigate the numerator of (15).Its support lives in a 2 dimensional subspace of R 4 , so we can project it onto the plane R 2 and find a strict separating hyperplane there.The projected support is precisely that depicted in Fig. 5(a).
Therefore, Theorem 2.3 holds for the running example and any of the two critical functions, and then, by Theorem 2.2 we conclude that the multistationarity region of the running example is path connected.
Note that a strict separating hyperplane of the support of f exists if the following system of linear inequalities has a solution (v, a) ∈ R k+1 : For the polynomial in Fig. 5, v = (2, 3), and a = 6 form a solution to the system.
In practice, we determine whether the system of linear inequalities ( 16)-( 18) has a solution as follows.First, we construct the polyhedral cone C ⊆ R k+1 defined by the inequalities ( 16)- (17).Second, we pick a point (v, a) in the relative interior of C. If there exists β ∈ σ − (f ) such that v • β > a, then (v, a) satisfies also inequality (18) and a strict separating hyperplane exists.If such β does not exist, then a simple argument gives that σ − (f ) is contained in the hyperplane defined by any (w, b) ∈ C, i.e. w • β = b for all β ∈ σ − (f ) and all (w, b) ∈ C. Therefore a strictly separating hyperplane does not exist.Theorem 2.3 gives a way to assert that the multistationarity region is path connected, but it is not informative if that is not the case.In [18] additional results are given to include two path-connected components.One of these results will be used to show that the multistationarity region of network in Fig. 2(c) has two path-connected components.
Model reduction for the simplification of the computations.Finding the critical function or the critical polynomial requires the computation of the determinant of a symbolic matrix, which can have a high computational cost if the matrix is large or the entries are long expressions in the symbolic variables.The next theorem shows that it is possible to remove certain reverse reactions from the network, and use a critical function for the reduced network to study the multistationarity region of the original network, thereby reducing (often dramatically) the computational cost.
Theorem 2.4 (Reduction and connectivity).Consider a conservative reaction network (S, R) without relevant boundary steady states.Assume that there exist species X 1 , . . ., X k ∈ S such that each X j participates in exactly 3 reactions of the form Let g • Φ be a critical function of the reduced network obtained by removing the reactions corresponding to κ 3j for j = 1, . . ., k. Assume that the closure of (g Then the number of path-connected components of the multistationarity region for both the reduced and the original reaction network is at most the number of path-connected components of (g In particular, if (g • Φ) −1 (R <0 ) is path connected, then the multistationarity region of the original network (S, R) is path connected.
The theorem might look a bit technical, but it is simply saying that it is enough to apply the algorithm to a smaller network obtained by removing the reverse reactions κ 3j , and the conclusions can be translated to the original network.Removal of reverse reactions contribute to the reduction of the computational cost as each of them gives an extreme vector to the flux cone (see Lemma 6.3(a) in Section 6).Making reversible reactions irreversible removes this extreme vector and thereby the matrix M (h, λ) depends on one less variable.
To illustrate Theorem 2.4, we consider the reaction network representing a signaling cascade with shared kinase in Fig. 2(i).This reaction network describes the phosphorylation of two substrates S and P with one phosphorylation site.The phosphorylation of S is catalyzed by a kinase E, while the phosphorylation of P is catalyzed both by E and by the phosphorylated form of S. The dephosphorylation processes are governed by two different phosphatases F 1 and F 2 [38].
Following Theorem 2.4, we remove the reactions corresponding to κ 3 , κ 6 , κ 9 , κ 12 , κ 15 .The reduced network has the form A matrix of extreme vectors is now The extreme vectors are the non-highlighted vectors in the matrix E above, with the entries corresponding to the reverse reactions removed (every third row of E).Computing the critical polynomial for the reduced network takes only 4 seconds.This critical polynomial is much simpler than the one for the full network.It has 15 variables and 204 terms.This example illustrates that the reduction in Theorem 2.4 might reduce the computational cost substantially.On one hand, the computation of the critical polynomial is faster and, on the other, the critical polynomial itself has less variables and terms, and therefore checking (Step 3) becomes faster as well.In the next section, we investigate networks where the benefit of applying network reduction and Theorem 2.4 is more dramatical.For example, for two of the networks, computing the critical polynomial for the full network turned out to be infeasible, but the computation became possible for the reduced network.By means of Theorem 2.4, we could assert connectivity of the multistationarity region for the full network (see Table 1 for more detail).This illustrates that Theorem 2.4 allows us to apply our approach to networks that were originally too large.
An important observation is that the existence of a strict separating hyperplane for a network or for a reduced version of it like in Theorem 2.4 are independent.That is, if we cannot find a strict separating hyperplane for the reduced network, it could still be that it exists for the original network.Also, the existence of this hyperplane depends on the choice of critical function, that is, of the parametrization.
Algorithm for path connectivity.We conclude this section by giving a procedure that checks a sufficient criterion for connectivity of the multistationarity region with no user intervention.Since most of the steps rely on solving linear inequalities, we implemented the algorithm using the computer algebra system SageMath [33].The code is given in the Supporting Information.We would like to emphasize that the multistationarity region could still be path connected, even if our algorithm terminates inconclusively.Output: 'The parameter region of multistationarity is path connected' or 'The algorithm is inconclusive'.

Investigating connectivity in relevant biochemical networks
We now demonstrate that Algorithm 2.5 is useful for realistic networks and that the number of connected components of the multistationarity region can be understood for several relevant networks in cell signaling of moderate size.
We start by going through the algorithm with two small networks: first, with the module regulating the cell cycle shown in Fig. 2(a), and then, with the simplified hybrid histidine kinase network in Fig. 2(b).The corresponding matrices and the critical polynomials become more complicated than for the running example, but still, are small enough to be displayed here.
Afterwards we analyze the rest of the networks in Fig. 2. Additionally, we consider the extensions of Fig. 2(e-f), with two phosphorylation sites, to several phosphorylation sites, and explore the strengths and weaknesses of the algorithm.When increasing the network size, the computation of the critical polynomial becomes unfeasible and we apply the reduction from Theorem 2.4.
We summarize the main properties of all the applications of the algorithm discussed in this work in Table 1.Table 1 shows the number of species, reactions, and extreme vectors of the reaction network.If the critical polynomial can be computed, it shows the number of positive and negative exponents of the critical polynomial and whether a strict separating hyperplane of the support exists.The same computations are repeated with the reduced network of Theorem 2.4, and we report the same data except the number of species and reactions.
3.1.Small networks.Cell cycle regulating module.We consider the model proposed in [39] for the second module that regulates the cell's transition from G2 phase to M-phase [40], which is shown in Fig. 2(a).This model has been analyzed for bistability in [30].With the order of species C, C + , M, M + , W, W + , the stoichiometric matrix, a matrix of conservation relations and a matrix of extreme vectors are: The sum of the three rows of W gives a positive vector and hence the network is conservative.We further verified that the network has no relevant boundary steady states.Furthermore, E has no zero row, so the network is consistent and the critical polynomial can be found using (13).
The matrix M (h, λ) is found by replacing the first, third and fifth rows of N diag(Eλ)A diag(h) by the rows of W : .
Hybrid histidine kinase.The hybrid histidine kinase network in Fig. 2(b) comprises a hybrid histidine kinase HK with the domain REC embedded, and separate histidine phospho-transfer domain Hpt.This reaction network has been studied in [41], where it was shown that the network displays multistationarity and a (complicate) description of the set of parameters with 3 steady states is given.It is further known that there is a choice of total concentrations such that the network is multistationary if and only if κ 3 > κ 1 , see [9] (this set is the projection of the multistationarity region onto the space of reaction rate constants).It was not known whether the full multistationarity region in κ and c is path connected.
With the order of species HK 00 , HK p0 , HK 0p , HK pp , Hpt, Hpt p , the stoichiometric matrix N , a matrix of conservation relations W , and a matrix E whose columns are a minimal set of extreme vectors are The network is conservative, consistent, and has no relevant boundary steady states.By replacing the first and fifth rows of N diag(Eλ)A diag(h) by the rows of W we find the matrix M (h, λ) and its determinant gives the critical polynomial: The polynomial has 8 variables and 17 positive and 2 negative coefficients.A strict separating hyperplane of its support is given by the equation (−5, −5, −1, 0, 5, 0, 0, 0) • µ = −3.
Using Theorem 2.2, it follows that the multistationarity region for the hybrid histidine kinase network is path connected.

Phosphorylation cycles.
We investigate models for phosphorylation and dephosphorylation of a substrate S with m binding sites with processes catalyzed by a kinase E and one or more phosphatases F .Sequential and distributive phosphorylation cycles.We first assume that phosphorylation and dephosphorylation occurs sequentially and distributively [42]: the kinase E catalyzes the phosphorylation one site at a time in a given order, while the phosphatase F dephosphorylates in the reverse order, also one site at a time.Under these assumptions, the network for m = 2 sites is shown in Fig. 2(e).
The dynamics of phosphorylation cycles have been intensively studied, e.g.[14,28,[42][43][44][45][46][47][48][49][50][51][52].In particular, it is known that they are multistationary for m ≥ 2, and there are choices of parameter values where they have m + 1 steady states for m even, and m for m odd [44], with half of them plus one being asymptotically stable [47].It has been conjectured that these networks can have up to 2m − 1 steady states, but this has only been established for small m [48].These networks are in the class of post-translational modification networks, which are conservative and consistent, and by the results in [25], since the underlying substrate network is strongly path connected, they do not have relevant boundary steady states.
The phosphorylation cycle with m = 2 binding sites has n = 9 species and the flux cone has = 6 extreme vectors.Then, the corresponding critical polynomial g • Ψ has 15 variables.It is a big polynomial with 288 positive and 112 negative exponents.Despite the large number of exponents, a strict separating hyperplane of its support can be found in less than a second.Our algorithm (and Theorem 2.2) can be applied to conclude that the multistationarity region for the sequential and distributive phosphorylation cycle with two binding sites is path connected.
By increasing the number m of binding sites, the reaction network size increases systematically.For example, for m = 3, the reaction network becomes S + E The critical polynomial g • Ψ has 2560 positive and 1536 negative exponents.The algorithm confirms that the multistationarity region is path connected in 96 seconds.
For m = 4, we could not compute the critical polynomial for the original network due to computer memory constraints.This shows that the computation of the critical polynomial is the bottleneck of the algorithm.After removing all reverse reactions as in Theorem 2.4, we could compute the critical polynomial, but it does not have a strict separating hyperplane.Therefore, for this family of networks we know that the multistationarity region is path connected for m = 2, 3, but it is unknown for m ≥ 4. Previous work has shown that the projection of the multistationarity region onto the reaction rate constants κ is path connected for all m ≥ 2, see [53], so we conjecture that the full multistationarity region is path connected for all m ≥ 2.
Phosphorylation cycles with different phosphatases.Different mechanisms for multisite phosphorylation have been observed, and in particular, phosphorylation and dephosphorylation of the different sites of a phosphate might not be catalyzed by the same kinase or phosphatase e.g.[54,55].If all steps are carried out by different enzymes, then multistationarity does not arise (see [38] for m = 2).Therefore, we consider the scenario where the phosphorylation occurs sequentially and the kinase acts in a distributive way, but we assume that each dephosphorylation step is governed by different phosphatases F 1 , . . ., F m (see Fig. 2(f) for m = 2).These networks are also conservative and do not have relevant boundary steady states for any m.
For m = 2 and m = 3, the algorithm finds a strict separating hyperplane, and hence the multistationarity region is path connected.For m = 4, the computation of the critical polynomial via the symbolic determinant in (Step 2) of the algorithm was too demanding, and the computer used for the tests ran out of memory.We employed the reduction approach given in Theorem 2.4 and removed eight reactions.The critical polynomial of the reduced network is significantly simpler: it has 22 variables and 178 monomials.Its support has a strict separating hyperplane.Thus, the multistationarity region of the original network is path connected by Theorem 2.4 and Theorem 2.3.
Weakly irreversible phosphorylation cycles.The two-site sequential and distributive phosphorylation network given in Fig. 2(e) assumes that each phosphorylation step proceeds via a Michaelis-Menten mechanism.This is referred to as strong irreversibility in [43].More plausible mechanisms have been argued to include the complex formation of the product with the enzyme, as for example a mechanism of this form would allow: A model incorporating this weak irreversibility at the dephosphorylation stage was proposed and analyzed for bistability in [56, Scheme 2].In the model, dephosphorylation of ERK by the phosphatase MKP3 proceeds as shown in Fig. 2(g) [57].For this network, our algorithm concludes that the multistationarity region is path connected.
A model with full weak irreversibility, that is, for both the phosphorylation and dephosphorylation processes, is shown in Fig. 2(h).The shape of the multistationarity region for this type of models was analyzed in [43], where it was concluded by means of a numerical approach that the multistationarity region in some aggregated steady state parameters is connected.Neither for the original network nor for the reduced network, a strict separating hyperplane exists, number of species; r = number of reactions; = number of extreme vectors of the flux cone; #σ ± (g • Ψ) = number of positive/negative exponents of the critical polynomial; existence of a strict separating hyperplane; ˜ = number of extreme vectors of the flux cone of the reduced network of Theorem 2.4; #σ ± (g • Ψ) = number of positive/negative exponents of the critical polynomial of the reduced network; existence of strict separating hyperplane for the reduced network.The number of variables of the critical polynomial is n + for the original network and n + ˜ for the reduced network.?? means that the computation could not be performed due to computer memory loss.The labels (a)-(k) refer to the networks in Fig. 2. and hence our algorithm is inconclusive.It remains thus open to be confirmed whether the multistationarity region is path connected.
Extracellular signal-regulated kinase (ERK) network.Dual-site phosphorylation and dephosphorylation of extracellular signal-regulated kinase has an important role in the regulation of many cellular activities [58], and a better knowledge of the dynamical properties of the ERK network might facilitate the prediction of this network's response to environmental changes or drug treatments [59].This network, analyzed in [60][61][62] and shown in Fig. 2(k), comprises as well phosphorylation of a substrate in two sites, but not in a distributive and sequential way.
Using Algorithm 2.5, we conclude that the multistationarity region for the ERK network is path connected.

Signaling cascades.
We next investigate two types of signaling cascades comprising phosphorylation cycles, which are known to be multistationary.Shared kinase.We consider first a two-layer signaling cascade with two single phosphorylation at each stage.The phosphorylated substrate of the first layer acts as the kinase of the second layer.We consider additionally that the kinase of the first layer also can act as kinase for the second layer, so the kinase is shared for the two layers, as shown in Fig. 2(i).Without this shared kinase, the cascade would not display multistationarity.
Algorithm 2.5 finds a strict separating hyperplane and we conclude that the multistationarity region for the cascade is path connected.Two-layer MAPK-cascade.Huang and Ferrell proposed a model for the MAPK cascade consisting of three layers, the first one being a single phosphorylation cycle, while the last two are dual phosphorylation cycles with phosphorylation and dephosphorylation proceeding sequentially and distributive [63].This network has bistability and also oscillations, and in fact for both properties only the first two layers of the cascade are required.
The network with two layers is shown in Fig. 2(j), and Algorithm 2.5 can be employed to conclude that the multistationarity region is path connected.
The full network with the three layers is large, with 22 species, 30 reactions, the rank of the stoichiometric matrix is 15, and the matrix E has 15 extreme vectors.Due to the computational cost, the computation of the critical polynomial was not possible for the original network, but was possible for the reduced network using the Julia package SymbolicCRN.jl[64].However, a strict separating hyperplane does not exist, so our algorithm is inconclusive.

3.4.
Reciprocal enzyme regulation.Finally, we consider two multistationary networks comprising single phosphorylation cycles, where the kinase and the phosphatase are subject to reciprocal regulation, both proposed and studied in [17].
Covalent regulation.We first consider the case where reciprocal regulation is via covalent modification catalyzed by the same enzyme, see Fig. 2(d).By means of Algorithm 2.5 we find that the multistationarity region for the the cascade is path connected.Allosteric regulation: An example with two path-connected components.The other mechanism of reciprocal regulation considered in [17] is via allosteric regulation: it is assumed that there is an allosteric effector L that binds both the phosphatase and the kinase, see Fig. 2(c).After performing a quasi-steady-state approximation, the authors of [17] show that a necessary condition for multistationarity is that κ 3 > κ 6 .This network is conservative, consistent, and has no relevant boundary steady states.
For all the applications we have seen so far, we could either conclude that the multistationarity region is path connected, or our approach was inconclusive.For this network our approach was inconclusive as well, however, by employing other theoretical results from [18] in conjunction with Theorem 2.4 and the approaches in [9], we were able to conclude that the multistationarity region has exactly two path-connected components, revealing two mechanisms underlying multistationarity.
The approach is as follows.On one hand, using the method to find parameter regions in κ from [9] relying on Theorem 2.1, we show in Section 6 that the two sets of parameters } both contain parameters that yield multistationarity, that is, these two regions both intersect the multistationarity region.To be precise, the condition κ 3 > κ 6 yields multistationarity if additionally the Michaelis-Menten constant for phosphorylation, K 1 = κ 2 +κ 3 κ 1 , is large enough relative to that for the dephosphorylation, K 2 = κ 5 +κ 6 κ 4 .Symmetrically, the condition κ 6 > κ 3 requires K 2 K 1 for multistationarity to arise.
We also show that if κ 3 = κ 6 , then the network cannot display multistationarity, no matter what the other parameters are.Therefore, any two points in each of the two sets above cannot be joined by a continuous path inside the multistationarity region, as any such path should cross at least one point where κ 3 = κ 6 .Hence, the full multistationarity region cannot be path connected: it has at least two path-connected components.
On the other hand, we show also in Section 6 that the multistationarity region of the reduced network has at most two path-connected components using [18].Hence, by Theorem 2.4, the original network in Fig. 2(c) has at most two path-connected components as well.
Putting the two pieces together, we conclude that the multistationarity region has precisely two connected components: one in which the catalytic rate of the phosphorylation step is larger than the catalytic rate of dephosphorylation, that is κ 3 > κ 6 , and the other where the inequality is reversed κ 6 > κ 3 .The second region was missed in [17] as it is outside the regime where the quasi-steady-state approximation employed there is valid.
To illustrate the type of switches that may arise from this system, we have considered the reduced model (for which there are also two connected components), and selected two parameter values α 1 , α 2 in different path-connected components of the multistationarity region.The two parameter values are identical except for the three parameters governing the dephosphorylation event: κ 4 , κ 6 and the total amount of P, see Fig. 6.We choose the path through the two points, Subfigures (a)-(c) show the bifurcation diagrams for the concentration of S p , K and L at steady state.In the intervals with three steady states, the one in the middle is unstable.Subfigure (d) shows the path in the three-dimensional space (K 2 , P tot , κ 6 ).The blue regions indicate the region of the path that belongs to the multistationarity region, and the displayed plane κ 6 = 1 separates the two regions.Subfigures (e)-(f) show the multistationarity regions when we fix κ 6 = 9.5 and κ 6 = 0.195 respectively.These are two slices of the two path-connected components, obtained by keeping only two free parameters.
which component-wise is given as 6(a-c) shows bifurcation diagrams, where t is the perturbed parameter, which perturbs simultaneously K 2 , κ 6 and the total amount of P, and we display the logarithm of the concentration of the phosphorylated substrate S p , kinase K and ligand L at steady state respectively.Note that by the choice of path we have made, t can take any real value, also negative.For the three concentrations, we obtain saddle-node bifurcations: a usual switch for larger values of t, while for smaller values of t, no hysteresis effect arises and the response curve has two components.In panel (d) of Fig. 6 a path in the three dimensional parameter space joining the two points is given.We see that the path enters and exists the multistationarity region twice, corresponding to the two path-connected components.The shape of these two components is displayed in Fig. 6(e-f) after slicing the three-dimensional space by fixing κ 6 for illustration purposes.

Discussion
Determining topological properties of semi-algebraic sets, that is, sets described by polynomial equalities and inequalities, is a highly complex problem that requires, for general sets, computationally expensive algorithms that scale poorly with the number of variables [15].For reaction networks with mass-action kinetics, the multistationarity region is a semi-algebraic set, and hence its description might not be straightforward.
Here we presented an approach to determine a basic topological property of a set, namely its connectivity.Non-connectivity of the multistationarity region may indicate different biological mechanisms underlying the existence of multiple steady states, and additionally, may give the cell the possibility to operate on complex switches as shown in Fig. 1.Our algorithm is to our knowledge the first to address the problem of connectivity in an effective and conclusive way.This is done by relying on linear programming and polyhedral geometry algorithms, rather than on semi-algebraic approaches, which reduces dramatically the computational cost.Additionally, our approach provides a symbolic proof of connectivity, and does not require numerical approaches, which unavoidably cannot explore the whole parameter space.
Although our algorithm might terminate inconclusively, even if the multistationarity region is path connected, we have shown that it is often applicable: for many motifs, the multistationarity region is connected because a strict separating hyperplane of the support of the critical polynomial exists.This came as a surprise to us, and might indicate a hidden feature present in realistic systems and brings up the question: What are the characteristics of the reaction networks from cell signaling that ensure that the support of the corresponding critical polynomial has a strict separating hyperplane?
It would certainly be relevant to understand the answer to this question, to bypass finding the critical polynomial, a step that is prohibitive for larger networks.This was illustrated with the networks of phosphorylation cycles with several phosphorylation sites, which revealed the computational boundaries of the algorithm.In several cases, the computation of the critical polynomial was not possible on a common computer.However, it was possible to compute the critical polynomial of the reduced network and still we were able to conclude that the multistationarity region is path connected.
Despite covering many networks, the algorithm remains inconclusive for some relevant networks, where we cannot decide whether the multistationarity region is connected, meaning that further investigations are required.For the m-site sequential and distributive systems, the projection of the multistationarity region onto the set of reaction rate constants is known to be path connected for all m [53], and this us believe that the same holds for the full region.In fact, based on the evidence gathered from the tested networks, we conjecture that if the projection onto the set of reaction rate constants κ is path connected, so is the multistationarity region.This would provide an additional strategy to study connectivity, which in particular might give a way to show that the fully weakly irreversible phosphorylation cycle studied in [4], see Fig. 2(h), is indeed path connected.
For the network in Fig. 2(c), where the algorithm was inconclusive, the network had two path-connected components.The strategy to show this was to combine knowledge about the projection of the multistationarity region onto the set of reaction rate constants, and a bound on the number of connected components of the reduced network found using ideas similar to those in the proof of [18,Theorem 3.9].This example opens up for new directions for counting path-connected components and understanding underlying features of reaction networks where the multistationarity region is disconnected.On one hand, it would be interesting to devise algorithms that can assert that the multistationarity region is disconnected, and ideally, count or give bounds on the number of path-connected components.On the other hand, one might wonder what network characteristics might give rise to disconnected multistationarity regions.For the reduced network of Fig. 2(c), the multistationarity region is no longer disconnected after deleting a reaction or a species, so this network can be viewed as a minimal motif with this property.A proper investigation of minimal networks with disconnected multistationarity region would require a better understanding on how the connectivity of the multistationarity region changes upon modifications on the network, in the spirit of Theorem 2.4.
We would like to point out that the proof of the key theorem, namely Theorem 2.2, is based on relating the multistationarity region to the preimage of the negative real half-line by the critical polynomial.When a strict separating hyperplane of the support of the critical polynomial exists, it is not only known that this preimage is path connected, but also that it is contractible [18,Theorem 3.9].This implies that all Betti numbers of the preimage set are zero.We conjecture that when this is the case, the multistationarity region is contractible, and hence topologically very simple (for example, it has no holes).However, this cannot be directly deduced by our arguments in the proof of Theorem 2.2.
To conclude, we propose the application of our work in the design of synthetic circuits displaying predefined switches.Indeed, by combining algorithms to determine multistationarity with the study of the connectivity of the multistationarity region, one can systematically study small networks and search for a desired input-output curve shape.This approach would adhere to related work already done in this direction, e.g.[65,66].

Methods
We implemented Algorithm 2.5 in SageMath 9.2 [33].A Jupyter notebook containing the code can be found in the Supporting Information.The computations for the networks in Table 1 were run on a Windows 10 computer with Intel Core i5-10310U CPU @ 1.70GHz 2.21 GHz processor and 8GB RAM.
In our implementation, each species is represented by a symbolic variable, and each reaction is represented by a list containing two symbolic expressions in the variables.For example, to run the code for the running example (1), one has to type: The output is given in the following format: The reaction network is conservative.There are no relevant boundary steady states.l = 2 Number of positive coefficients: 3 Number of negative coefficients: 1 The support set has a strict separating hyperplane.All the conditions are satisfied.We conclude that the parameter region of multistationarity is path connected.Although we used SageMath to compute the extreme vectors, the same computation can also be done with Polymake [67] (as a standalone [68] or as part of the Oscar project in Julia [69]).These programs compute extreme vectors of arbitrary pointed cones.For open cones of the form ker(N ) ∩ R n >0 , there exist specific algorithms designed in the context of stoichiometric network analysis and metabolic network analysis.See for example [70].
As discussed above, the bottleneck of Algorithm 2.5 is to compute the determinant of the symbolic matrix M (h, λ) in (13).In the tested examples, the Julia package SymbolicCRN.jl[64] is more efficient in computing the symbolic determinant than our implementation using SageMath.Using SymbolicCRN.jl,we were able to compute the critical polynomial for the reduced MAPK cascade with three layers, which was not possible with SageMath.For the 4-site phosphorylation networks in Table 1, we could not compute the critical polynomial using neither Julia nor SageMath.
The parameter regions of multistationarity displayed in Fig. 6(e-f) are found using cylindrical algebraic decomposition in Maple, using the command CellDecomposition from the packages RootFinding[Parametric] and RegularChains.

Proof of the results
6.1.Convex parameters.In this subsection we show basic results on the flux cone, and specially that positive combinations of extreme vectors of ker(N ) ∩ R r ≥0 parametrize the positive part of the cone.Proof.First, we observe that for (x, κ) ∈ V the vector v κ (x) lies in ker(N ) ∩ R r >0 by (3).Now, part (b) follows directly from Proposition 6.1(b).To show that Ψ is surjective, let (x, κ) ∈ V.By Proposition 6.1(a), there exist λ ∈ R >0 such that v κ (x) = Eλ.By letting h = 1/x, using the definition of v κ (x) in (3) and the definition of Ψ, one easily sees that Ψ(h, λ) = (x, κ), and hence, (x, κ) is in the image of Ψ, concluding the proof.Lemma 6.3.Let (S, R) be a reaction network such that the last two reactions R r−1 and R r are reverse to each other.Let ( S, R) be the reduced network obtained by removing the reaction R r .
(a) The vector (0, . . ., 0, 1, 1) ∈ R r ≥0 is an extreme vector of the flux cone of (S, R).(b) If v is an extreme vector of the flux cone of ( S, R), then (v, 0) is an extreme vector of the flux cone of (S, R).
In particular, the flux cone of (S, R) has more extreme vectors than the flux cone of ( S, R).
Proof.Before starting the proof, we recall some definitions and statements.The support of a vector u ∈ R r is the set of indices where the vector is non-zero: supp(u) := {i ∈ {1, . . ., r} | u i = 0}.A vector u ∈ ker(N ) ∩ R r ≥0 is said to have minimal support if for all w ∈ ker(N ) ∩ R r
(b) Let Ñ ∈ R n×(r−1) denote the stoichiometric matrix of the reduced network and hence we have that (v, 0) is contained in the flux cone of (S, R).To show that (v, 0) is an extreme vector, we show that it has minimal support.Let w ∈ ker(N ) ∩ R r ≥0 such that supp(w) ⊆ supp((v, 0)).Then w n = 0, and hence (w 1 . . ., w r−1 ) ∈ ker( Ñ ) ∩ R r−1 ≥0 .Since v is an extreme vector of ker( Ñ ) ∩ R r−1 ≥0 and supp((w 1 . . ., w r−1 )) ⊆ supp(v), it follows that supp(w) = supp((w 1 . . ., w r−1 )) = supp(v) = supp((v, 0)).Hence (v, 0) has minimal support and is therefore an extreme vector.6.2.Proof of Theorem 2.2.Theorem 2.2 is a direct consequence of two technical lemmas.To state the lemmas, we use the notation of the main text, and additionally we write Ω ⊆ R r >0 × R d for the parameter region of multistationarity: We will write Ω as the image of a subset of V under the map and then compare path-connected components.Recall that W is a matrix of conservation relations.
By Theorem 2.1, Ω is closely related to the preimage of the negative real half-line under the critical function given in (8) for a parametrization Φ : D → V: So we introduce the set: We summarize all the relevant functions that play a role in the proof of Theorem 2.2 in the following diagram: We denote by b 0 (X) the number of path-connected components of a set X ⊆ R k [75, Definition 3.3.7].Note that X is path connected if and only if b 0 (X) = 1.Lemma 6.4.For a dissipative reaction network without relevant boundary steady states, it holds that π(Θ) = Ω.
The second part follows from the fact that continuous images of path connected sets are path connected [75,Theorem 3.3.5].Lemma 6.5.Consider a dissipative reaction network without relevant boundary steady states.If the closure of (g Proof.For all (x, κ) ∈ V with g(x, κ) < 0, Theorem 2.1(B) gives that π(x, κ) ∈ Ω.Thus, by definition of Θ, it holds that By taking preimages under Φ, it follows that Using the Curve Selecting Lemma [76], there exists a continuous path γ : [0, 1] → R n such that γ(0) = η and γ (0, 1) ⊆ (g •Φ) −1 (R <0 ).Thus, there exists a continuous path between η and one of the path-connected components of (g Since Φ is surjective, it follows that Φ(Φ −1 (Θ)) = Θ.Since a continuous image of a path connected set is path connected [75,Theorem 3.3.5],we conclude that b 0 (Θ) ≤ b 0 (Φ −1 (Θ)).6.3.Proof of Theorem 2.4.To prove Theorem 2.4, we need to show that under the hypotheses of the theorem, we have b 0 (Ω) ≤ b 0 ((g • Φ) −1 (R <0 )).The proof is based on inductive application of the following statement.Proposition 6.6.Let (S, R) be a conservative reaction network without relevant boundary steady states.Assume that the species X n participates in exactly 3 reactions of the form Let ( S, R) denote the reduced network obtained by removing the reaction corresponding to κ r and Θ denote the set in (20) for ( S, R).It holds that b 0 (Θ) ≤ b 0 ( Θ).
Proof.For every object corresponding to a network, we write ˜to indicate that it corresponds to the reduced network, e.g. the reaction rate constants in the reduced network are denoted by κ1 , . . ., κr−1 .
As in the proof of Lemma 6.3, the stoichiometric matrices N and Ñ satisfy N r−1 = −N r , Ñj = N j for j = 1, . . ., r − 1, and rk(N ) = rk( Ñ ).Recall that W ∈ R d×n is a row reduced full rank matrix such that W N = 0.It follows that W Ñ = 0, so W is a matrix of conservation relations for ( S, R) and we use this matrix in the definition of π, analogously to (19).Using (9), we also have that ( S, R) is conservative if and only if (S, R) is conservative.
Following the proof of [77,Prop. 1] we introduce the maps: n−1 where a = (a 1 , . . ., a n−1 , 0).We start by relating the steady states of the two networks under the hypothesis that η(κ) = η(κ).Recall that (κ, c) ∈ Ω if and only if the equation system has at least two positive solutions.Redundant linear relations among the equations f κ (x) = 0 arise from linear dependencies of the rows of N .To remove these redundancies, we consider I = {i 1 , . . ., i d } to be the set of indices of the first non-zero coordinates of each row of W , i 1 < • • • < i d .For all j = 1, . . ., d, we replace the i j th row of f κ (x) by the jth row of W x − c.If we denote the resulting function by h κ,c (x), then x * ∈ R n >0 is a solution of ( 21) if and only if h κ,c (x * ) = 0 holds.Thus, the parameter pair (κ, c) enables multistationarity if and only if h κ,c (x) = 0 has at least two positive solutions.
Since x n appears linearly in h κ,c (x), there exist vectors z(κ), v(x, κ) such that where u(x, κ) is chosen such that (22) holds.
We define hκ,c (x) analogously for the reduced network and write hκ,c (x) = z(κ) −κ r−2 x n + ṽ(x, κ) −κ r−1 x a , which under the assumption that η(κ) = η(κ), we have It then holds that (23) For i ∈ I, it is straightforward to see that this equality holds.If i / ∈ I, then the equation reduces to the equality With this in place, consider the matrices , and B(κ where Id n−1 is the identity matrix of size n − 1.Using (23), we have Since the matrices B(κ) and B(κ) are invertible, it follows from (24) that every positive solution of h κ,c (x) = 0 is a solution of h κ,c (x) = 0.In particular, the reduced network does not have relevant boundary steady states.Additionally, since s = s, and the Jacobian of h κ,c (x) is precisely the matrix M κ (x) (and analogous for the reduced network), taking the Jacobian and determinant of both sides of (24) yields: Since κ r , κ r−2 , κr−2 > 0, it follows that g(x, κ) and g(x, κ) have the same sign if η(κ) = η(κ).
We aim at showing that the multistationarity region Ω of the network in Fig. 2(c) has exactly two path-connected components.We do this in two steps.First, we show that b 0 (Ω) ≥ 2, and then that b 0 (Ω) ≤ 2, giving the equality.
The coefficient of the other monomials of p κ (ξ) are sums of products of the parameters, and are positive.Thus, the value of the critical function (g • Φ)(ξ, κ) is positive for all ξ ∈ R 4 >0 if κ 3 = κ 6 .Now, one can apply part A of Theorem 2.1 to conclude that (κ, c) does not enable multistationarity if κ 3 = κ 6 .
Hence, for any ξ 4 > 0 large enough, this polynomial is negative, and so is q κ (ξ).Therefore, p κ (ξ) also attains negative values.If κ 3 − κ 6 < 0, all we need to do is to let K 1 be small enough, or K 2 large enough and repeat the argument.We conclude that b 0 (Ω) ≥ 2.
Showing that b 0 (Ω) ≤ 2. We now apply Theorem 2.4 to the reduced network:  A matrix of extreme vectors is
For each s ∈ [0, 1], we define the function f s : R >0 → h −1 (R <0 ) as c α γ(s) α t v•α , = h(γ(s)) + p s (t), where γ(s) • t v = (γ 1 (s)t v 1 , . . ., γ 12 (s)t v 12 ), and from (31) p s (t) is a sum of generalized monomials in t with all exponents negative and all coefficients positive.Hence the leading coefficient of f s (t) is h(γ(s)) < 0, all the other coefficients of f s (t) are positive, and f s (t) has a unique positive real root T s > 0. Since the roots of a polynomial depend continuously on the coefficients, T s depends continuously on s.Therefore, T := max s∈[0,1] T s exists.

Supporting information
An accompanying Jupyter notebook contains the code of the algorithm, written in SageMath 9.2 [33].

Figure 1 .
Figure 1.(a-c) Input-output curves for hypothetical systems.Input is thought to be

Figure 3 .
Figure 3. Illustration of the relevant objects for the running example given in (1).

Algorithm 2 . 5 .
Input: a reaction network (Step 1) Check that the reaction network is conservative and that it does not have relevant boundary steady states using siphons.(Step 2) Compute the convex parametrization map Ψ if the network is consistent, and the critical polynomial g • Ψ from (13).(Step 3) Decide whether a strict separating hyperplane of the support of g • Ψ exists.(Step 4) Eventually repeat Steps 1-3 with a reduced network as in Theorem 2.4.
(a) If E does not have any zero row, then the convex parametrization map Ψ : R n >0 ×R >0 → V from (11) is surjective.(b) If E has a zero row, then V = ∅.