Characterizing Multistationarity Regimes in Biochemical Reaction Networks

Switch like responses appear as common strategies in the regulation of cellular systems. Here we present a method to characterize bistable regimes in biochemical reaction networks that can be of use to both direct and reverse engineering of biological switches. In the design of a synthetic biological switch, it is important to study the capability for bistability of the underlying biochemical network structure. Chemical Reaction Network Theory (CRNT) may help at this level to decide whether a given network has the capacity for multiple positive equilibria, based on their structural properties. However, in order to build a working switch, we also need to ensure that the bistability property is robust, by studying the conditions leading to the existence of two different steady states. In the reverse engineering of biological switches, knowledge collected about the bistable regimes of the underlying potential model structures can contribute at the model identification stage to a drastic reduction of the feasible region in the parameter space of search. In this work, we make use and extend previous results of the CRNT, aiming not only to discriminate whether a biochemical reaction network can exhibit multiple steady states, but also to determine the regions within the whole space of parameters capable of producing multistationarity. To that purpose we present and justify a condition on the parameters of biochemical networks for the appearance of multistationarity, and propose an efficient and reliable computational method to check its satisfaction through the parameter space.


Introduction
Multistability is a nonlinear phenomenon characterized by the existence of two or more stable steady states, where a given dynamical system will evolve depending on its initial conditions. Important biological phenomena, like cellular decision processes, rely on multistable models, where the different functional phenotypic states or cell fates can be understood as discrete, stable and mutually exclusive stable states [1]. Experimental evidences for bistability have been found in numerous pathways involved in cell decision processes, such as the p42 MAPK/Cdc2 network governing the maturation of oocytes in Xenopus [2], the pheromone sensing MAPK pathway in S. cerevisiae [3], or the Rb-E2F pathway regulating proliferation in mammalian cells [4]. The analysis of mathematical models of the underlying multistable networks contributes to understand these biological phenomena from a systems perspective. In [5] for example, basic design principles for the control of the cell cycle are suggested based on the modeling and analysis of the gene circuit underlying the Rb-E2F switch, identified by the criterion of robustness.
The dynamics of biochemical reaction networks (i.e., the time evolution of the vector of species concentrations) can be described by models of coupled ordinary differential equations where the structure depends on the reaction connectivities, stoichiometry and kinetics, and the parameters are defined from kinetic rate constants. Modeling a biochemical system consists of inferring the structure and parameters from experimental data. For a given model structure, the corresponding parameters are typically estimated from experimental time course measurements of observables (usually linear combinations of some subsets of the species concentrations), by minimizing some measure of error between the experimental data and the model prediction [6]. Nontraditional methodologies for the determination of reaction mechanisms from kinetic data sets have been reviewed by [7,8]. In processes occurring within living cells, the access to quantitative information is often very limited, and this fact has severe implications on the development of mathematical models: it hampers model discrimination, often leads to poor parameter identifiability, and makes the parameter estimation task very challenging, since it entails to solve a nonconvex optimization problem in high dimensional search spaces [9], which cannot be reduced and/or constrained in absence of a priori knowledge about feasible parameter values.
One of the challenges of systems biology is to provide tools to overcome the lack of quantitative information, by exploring and exploiting connections between model structure and/or parameters with the expected dynamic behaviour [10][11][12]. In this context, methods to systematically detect multistationary regimes in biochemical systems will help modeling multistable systems, constraining the feasible parameter regions, for a given structure, based on the capability to produce multistationarity. These methods are also of great interest in the design of synthetic biological switches, where the robustness of the multistationarity property needs to be analyzed [13][14][15].
Current results in this direction are derived from different fields, from classical bifurcation theory [16] to monotone systems [17]. In particular, structural properties of reaction networks and their connection with multistationarity are at the core of the Chemical Reaction Network Theory, pioneered by Feinberg, Jackson and Horn [18,19] and subject to ongoing development [20,21] with special interest in the application to biological systems [22][23][24]. In the context of cell signaling, for example, CRNT has been used to discard kinetic mechanisms based on their capacity for multistationarity [25,26].
The deficiency one algorithm, the advanced deficiency theory, the deficiency zero and deficiency one theorems are part of the CRNT in which networks are classified by means of a nonnegative integer index called deficiency -a property of the graph of complexes of a network-and some structural conditions are evaluated to decide whether networks have the capacity for multiple positive equilibria [27].
As pointed out in [28], CRNT provides surprisingly strong results for reaction networks based only on the systems structure. For example, the deficiency zero theorem asserts that every weakly reversible network of zero deficiency has a unique equilibrium, for any choices of parameter values. However, when multistationarity cannot be ruled out, nothing can be said about how the parameters affect the qualitative behaviour of the solutions.
In a previous paper [29] we have introduced the parameters into the picture, providing a canonical expression for the equilibrium manifold in terms of the kinetic parameters and the so called deficiency parameters of the network. The concept of network layout, introduced in [29] as the difference between the deficiency of a network (d) and the dimension of the equilibrium manifold (l), allowed us to classify biochemical reaction networks in three groups: proper networks (d~l), overdimensioned networks (dwl) and underdimensioned networks (dvl). The analysis presented focused on proper networks, i.e. those networks where deficiency and manifold dimension coincide. For those networks, the qualitative behaviour of the manifold was evaluated through its derivative with respect to the deficiency parameters. Then a geometric intuitive idea was applied to find, under these restrictive assumptions (d~l), a condition on the parameters of the network giving room to multiple steady states. The condition was formulated as a (non convex and multimodal) optimization problem, to be solved by global optimization algorithms. Before the search, it was needed to partition the parameter space in regions with different qualitative behaviour. The parameter space was then characterized depending on whether the optimization algorithm could find a solution or not, i.e., depending on whether or not we could find a point (or a set of discrete points) in the parameter space fulfilling the multistationarity condition.
In this work, we formalize the geometric intuition in [29] to derive a sufficient condition for multistationarity which is general, i.e. applies to networks where the dimension of the equilibrium manifold is lower, equal or greater than the deficiency. The general condition is stated in a formal context, and a proof of its validity in arbitrary dimensional spaces is provided. In addition, we present a method to search the condition through the parameter space with several advantages over global optimization methods. On the one hand, it is capable of finding the regions where multistationarity condition is fulfilled, without requiring the a priori partitioning of the parameter space in areas of different qualitative behavior. On the other hand, the method allows the characterization of the multistationarity regimes of a biochemical reaction network in a reliable manner, i.e. those regions in the parameter-state space leading to multiple steady states are ensured to be enclosed by the solution set.

Fundamentals
Following CRNT classical description (see Feinberg's Lecture 3 in [30] for details), we consider a generic reaction network involving m species fS 1 , . . . ,S m g participating on a given set of irreversible reaction steps. Their concentrations c i are collected on a vector c defined on the space R m §0 we will refer to as the species space. In the following, we write x[R m w0 if x i w0 for all i~1,,m and Each reaction step will be represented by an arrow which connects two particular combination of species, thus indicating how a given set of reactants is transformed into a certain set of products. The set of species at both extremes of the arrows are known in CRNT as complexes. The set of all complexes connected by reaction steps conforms the reaction network which is represented by a directed graph (Cgraph or graph of complexes), where arrows (edges) indicate the reaction steps and the nodes correspond with the complexes (see Figure 2 for an example).
Let fC 1 , . . . ,C n g be the set of complexes of the network. To each complex C i we associate a set I i of integer elements which collects the indexes of those complexes that are directly reached from C i and a pair of vectors fy i ,e i g. The set I i can be formally defined as: §0 contains the molecularities of the species in complex i, and e i is a vector of the standard basis of R n such that for every i,j~1,, n: The complete set of edges in the graph is constructed by connecting C i ?C Ii for all complexes i~1, . . . ,n. Every edge in the graph (or reaction step) directly linking complex C i to complex C j , has its corresponding reaction rate, of the form: where k ij w0 is a constant parameter and y i a scalar function: In what follows we assume that the reaction rates are mass action. Thus each function y i (c) takes the form: Provided that c[R m w0 (i.e. it is a strictly positive vector), the expression (3) can also be written as: where the natural logarithm operator ln ( : ) acts on any vector element-wise.
The C-graph of a reaction network is composed by a number ' of ''isolated'' sub-graphs known in CRNT as linkage classes fL 1 , . . . ,L ' g, (see Lecture 3, page 14 in [30] for a complete discussion), each containing a number of complexes h k so that: The reaction network presented in Figure 2 consists of two linkage classes: one involving complexes fC 1 ,C 5 g, and the other involving complexes fC 2 ,C 3 ,C 4 g. Each linkage class L k is accompanied by a vector L k [R n §0 of the form: Complexes within a linkage class are linked by sequences of reaction steps (or equivalently by sequences of arrows) defining paths. Two complexes are strongly linked if they can be reached from each other by directed paths (every complex is strongly linked to itself). The maximal set of strongly linked complexes is a strong terminal linkage class if no other complex can be reached from any of its elements. A linkage class that is also a strong terminal linkage class is said to be weakly reversible. Weakly reversible networks are those composed of weakly reversible linkage classes. The network depicted in Figure 2 is an example of a weakly reversible network.
The dynamics of chemical reaction networks. Making use of the above definitions, the time evolution of species concentrations can be described by a set of ordinary differential equations of the form This system can be rearranged into the more familiar form, extensively employed in the context of CRNT: where vector y(c)[R n §0 contains as entries the scalar functions y i (2). In what follows, we will refer to the space R n §0 in which every y is defined as the space of complexes. The mapping A : R n §0 ?R n is such that: and Y in (7) is the m|n molecularity matrix, with columns being the molecularity vectors y 1 , . . . ,y n of complexes C 1 , . . . ,C n . The right hand side of equation (7) f( : ) : R m §0 ?R m is known in CRNT as the species formation function [18].
Trajectories of system (7) are constrained in the concentration space by some invariants of motion to lie on convex regions within the non-negative orthant known as reaction polyhedrons [31], or stoichiometric compatibility classes in the CRNT formalism (see Lecture 2, page 17, Definition 2.9 in [30]).
In order to characterize these regions, and suggested by the structure of the operator A (Eq. (8)), let us first consider the subspace spanned by vectors (e j {e i ), where the image of the operator A lies. The subspace is formally defined as: It must be noted that the number of independent vectors (e j {e i ) per linkage class is h k {1, and vectors from different linkage classes are independent, thus the dimension of D is n{'. Related to D there is another subspace we define as follows: The subspace S is known in CRNT as the stoichiometric subspace, where the species formation function f(c) lies, since from (7) and (8) we also have that: Let s be the dimension of S and B[R (m|m{s) a full rank matrix that spans column-wise the orthogonal complement S \ .
A linear vector-valued function of the form W (c)~B T c characterizes the invariants of motion, since it is constant along the trajectories defined by (7) for any initial condition c(0)~c 0 . This can be shown by differentiating W (c) along (7) so that: where the last equality follows because by construction B is orthogonal to the stoichiometric subspace S. Integrating this expression in time and taking into account the initial condition c 0 we get that: Each equation in the above expression is called a conservation law [32]. In this way, trajectories are constrained to regions that result from the intersection of the non-negative orthant R m §0 with any linear variety associated to the stoichiometric subspace. These regions we will refer to as reaction polyhedrons (equivalently stoichiometric compatibility classes), can be formally defined with respect to a reference concentration vector c 0 as: Mass conservation is a special case of invariance which leads to a particular class of polyhedrons, namely those which intersect every axis in the concentration space (in this case every species participate in at least one conservation law). Other polyhedrons are possible being parallel to some axis as it is the case of the example. In the extreme case (no species participating in conservation laws) the stoichiometric subspace spans the whole concentration space so that the orthogonal complement is the zero vector (the only column of B). In this case the reaction polyhedron is R m §0 . The nature of equilibrium points. Next we summarize some results from CRNT to be used in the contribution concerning possible equilibrium solutions of (7), namely vectors In what follows, we will concentrate on weakly reversible networks which in addition, for any initial condition c §0 produce equilibrium points in the interior of the positive orthant. Note that if the trajectories lie in the interior of the positive orthant the networks under study are persistent [32] (i.e. those which for any initial condition cw0 produce trajectories lying in R m w0 ). Any equilibrium point c Ã [R m w0 for (7) will satisfy either: or Related to condition (14) there exists a subspace D d : ker Y \ImA, we will refer to as the deficiency subspace, which plays a central role in CRNT. Its dimension is called the deficiency d and can be computed by making use of the standard relationship between the dimensions of the domain, kernel and image of a linear transformation (see [30], Lecture 4, page 23 in the proof of Proposition 4.7). For weakly reversible networks, the terms in the relationship coincide with the dimensions of subspaces D (see (9)), D d and S:ImYA, respectively. The final expression then reads: since from previous discussion dim(D)~n{' and dim(ImYA)~s, we finally get for dim( ker Y \ImA): It must be noted that as stated in the so-called deficiency zero theorem [30] (Theorem 5.1 in Lecture 5, page 2), any equilibrium point fulfilling (13) will be stable and unique in each compatibility class. Consequently, multistationarity, that is to say multiple equilibria, can only occur for c Ã satisfying (14). For this reason, any weakly reversible reaction network of zero deficiency possesses a unique and stable equilibrium point (associated to D 0 ) per stoichiometric compatibility class. This result remains valid independently of the values taken by the reaction rate constants.
In a previous paper [29], we have exploited the graph structure of biochemical networks to obtain an expression of the locus of equilibria -the set of points c Ã such that f(c Ã )~0 in (7)-in terms of as many parameters as the deficiency of the network. For a class of networks (the so called proper networks) we were able to partition (by continuation of variation parameters associated with the deficiency) the space of kinetic parameters in regions with different qualitative dynamic behavior. Exploiting the concept of equilibrium as a intersection between solutions satisfying (14) with the so called reaction polyhedron [31], conditions on the parameters of the network leading to multiple steady states were found for networks where the manifold dimension and deficiency coincide.
Here we start from this insight to provide a general condition for the existence of multistationarity, valid for weakly reversible mass action networks.

The Locus of Equilibria
We present in this section a canonical expression for the locus of equilibria, that is to say the set of all possible feasible equilibrium solutions in terms of the kinetic parameters of the network. In what follows, we will refer to this locus as the equilibrium manifold.
Mathematically it corresponds with an algebraic variety which results from the intersection of two other varieties: the family of solutions and the mass action manifold to be described below.
In preparing for the description and for simplicity, we assume that the molecularity matrix Y is full rank and mƒn. Furthermore, we assume that the m independent molecularity vectors y i are distributed among linkage classes so that each linkage has at least one independent vector. For such networks, and without loss of generality, let us number the first ' complexes so that each belongs to a different linkage class (note that by previous assumption m §'), and so that the first m columns of the matrix Y are linearly independent.
The family of solutions is a linear variety F defined in the space of complexes by: where a j are given real numbers, and vectors x k , f j [R n §0 are solutions of the following equations: A(f j )~w j j~1, . . . ,d: The set of vectors fw j g d j~1 that appears in (19), defines a basis for the deficiency subspace D d . As proved in [29], the elements of a basis for D d can be obtained from the left kernel of the matrix: where Y is the molecularity matrix and L is the n|' matrix with columns being the vectors L k defined in (5). Vectors x k [R n §0 for k~1,, ' constitute a basis for the kernel of A (Eq. (8)). As stated in Proposition 4.1 of [30] (Lecture 4, page 10), for weakly reversible networks the dimension of the kernel of A, and therefore the number of vectors of the basis, coincides with the number of linkage classes. Actually, the same holds for networks other than weakly reversible, provided that they have one terminal linkage class per linkage class. The same Proposition prescribes for each element of the basis a nonnegative vector. In particular, each vector x k associated to a linkage class L k , will have entries of the form: where parameters r ik w0 are functions of the reaction constants within the linkage class L k . The functional relation will be formally represented as r ik (K)w0 where K is the vector of kinetic rate constants. A proof of Proposition 4.1 is out of the scope of this contribution. It cannot be found in [30] either. However, two alternative proofs can be found in the literature: one based on Perron-Frobenius theorems applicable to weakly reversible networks [33]. The other more graph theoretically oriented has been proposed in [34].
Similarly, any entry f ij (K) of vectors f j is a function of kinetic parameters within the linkage class containing the complex C i . Explicit expressions relating parameters r ik and f ij (for i~1, . . . ,n, j~1, . . . ,d) with the original kinetic constants are obtained by solving Eqs. (18) and (19).
It must be noted that by construction, F in Eq. (17) under linear transformation A produces vectors: which correspond with elements of the deficiency subspace D d .
For the particular case of a i~0 for every i~1,, d, F characterizes elements of D 0 (i.e. the kernel of A). In this way, Eq. (17) provides a complete parametrization of complexes y, leading to equilibrium solutions.
The mass action manifold is a nonlinear algebraic variety M defined in the space of complexes by: where matrix Q is of the form: with Y 1 being a matrix containing the first m columns of the molecularity matrix Y . Note that each element M i can be written as: M i~P m j~1 y q ij j : Expression (23) is the equivalent of (4), but defined in the space of complexes instead of the concentration space. The relationship between both spaces is given by the bijective mapping: The equilibrium manifold is the algebraic variety which results from the intersection in the space of complexes of the family of solutions F and the mass action manifold M. Formally the intersection can be expressed as: where a is the vector that collects all the a i parameters for i~1, . . . ,d. For a given parameter vector K, the equilibrium manifold can be written as: where On occasions it may be convenient to transform the algebraic variety by means of (24) to its equivalent in the concentration space, namely: where: It should be noted that given a rate constant vector K, function H s (or equivalently H c ) is continuous and differentiable. In addition the dimension of the expression (28) (or equivalently (26)), either in the space of complexes or in the species space is l~m{s. This is so since l~mzd{(n{'), and by (16) d~n{'{s.

Condition for Multistationarity
This section contains the main result of the contribution namely a condition for a given network to have multiple (positive) equilibria within the same stoichiometric compatibility class. The condition is geometric in essence and makes use of the equilibrium manifold, expressed either in the space of complexes (26) or in the species space (28), and the linear variety that defines the set of reaction polyhedrons (12).
The underlying idea behind the condition relates to the question of whether or not a given set of equations can accept more than one solution. A formal statement and discussion of this question is presented in the Appendix S1 for a general class of functions defined in R n . There we present the necessary mathematical background as well, but first let us illustrate the basic concept on a simple two dimensional case depicted in Figure 1.
The example consists of a nonlinear manifold (the continuous curve) and two possible families of linear varieties represented by dashed lines with gradient vectors B and B' respectively. Vector P, perpendicular to the curve at x Ã , is the gradient of the curve at x Ã . This vector also defines the tangent subspace to the curve at that point. As it can be observed in Figure 1 A, all linear varieties associated to vector B intersect the curve in just one point so no multiple solutions are expected.
On the other hand, some linear varieties associated to B' in Figure 1 B intersect the curve in two points y 1 and y 2 what corresponds with two different solutions. What differs between Figures 1 A and 1 B is the relative orientations of the curve and linear variety gradients. Thereby vector alignment (or linear dependency) is what seems to be at stake to determine the number of solutions. In this way, multiple solutions are expected to appear whenever vector alignment takes place. This is the notion we take advantage of and extend to higher dimensional manifolds (i.e. the equilibrium manifold).
As discussed in the Appendix S1, essentially all is needed is the equilibrium manifold to be locally smooth. If at a given point in the space of concentrations this is the case, hyperplanes secant to the equilibrium manifold can be constructed in a small neighborhood of the point by parallel translation of the corresponding hyperplanes tangent at that point (see Figure S1). Multiple solutions are then possible if the hyperplanes coincide with a given reaction polyhedron.
In the remaining of the section the results discussed in the Appendix S1 will be adapted to detect multiple equilibria associated to the equilibrium manifold. To that purpose a mzd (equivalently n{'zm{s) dimensional space will be employed, which includes the variables: (c,a)~(c 1 ,c m ,a 1 ,,a d ): It must be pointed out that since the map (24) is bijective, the condition can be established either in the species space (c variables) or in the space of complexes (y m variables). However in the species space the reaction polyhedron is linear what simplifies the derivation of the condition. Thus for convenience, manifold (28) will be employed in first place. Nevertheless, some comments will be made at the end of the section on the condition expressed in the the space of complexes.
Firstly let us note that the Jacobian of H s reads: where D c H s [R (n{')|m and D a H s [R (n{')|d denote the Jacobians of H s with respect to c and a, respectively. In the space described by the variables (c,a) we express the linear variety associated to the reaction polyhedron (12) as: where function W : R m w0 |R d ?R m{s is of the form: Note that because D c W~B T and D a W~0 (i.e. a zero (m{s)|d matrix), its Jacobian can be written as: DW is full rank by construction, since as discussed in the Fundamentals, B is a basis for the orthogonal complement of the stoichiometric subspace S. For a given vector of rate constants K, let H s (c,a; K)~0 be continuous in the vicinity of a point (c Ã ,a Ã ). Note that by the implicit function theorem, this implies that DH s (c Ã ,a Ã ; K) is full rank. Then we are under the conditions of Proposition A1 (see Appendix S1) where H s corresponds with F(x), x:(c,a) and P~DH T s (c Ã ,a Ã ; K). Furthermore, DW T takes the place of C in matrix G in Corollaries A1 and A2. The corresponding n{'zm{s square G matrix then becomes: We are now in the position to formally state the geometric condition. Proposition 1 Consider a reaction network with a given vector of rate constants K, and let H s (c,a; K)~0 be continuous in its domain. If for any (c,a) satisfying H s (c,a; K)~0 matrix G is full rank, the reaction network for K has at most one positive equilibrium solution per stoichiometric compatibility class.
Proof: the result follows directly from Corollary A1 (see Appendix S1), applied to the domain where the equilibrium manifold (29) is defined (i.e. positive concentration space).
Proposition 2 Given a vector of rate constants K, a sufficient condition for the reaction network to exhibit multiple (positive) steady states within the same stoichiometric compatibility class is that for at least some (c Ã ,a Ã ) such that H s (c Ã ,a Ã ; K)~0, matrix G(c Ã ,a Ã ; K) is rank deficient. Equivalently, a sufficient condition for the reaction network to exhibit multistationarity is that: for some (c Ã ,a Ã ). Proof: The result follows directly from Corollary A2 (see Appendix S1), applied to the domain where the equilibrium manifold (29) is defined (i.e. positive concentration space). Rank deficiency can be checked through expression (35).
The condition for multistationarity remains valid in the space of complexes, since the map (24) is bijective. In this space however, W (y m ; c 0 )~0 is nonlinear, although continuous for every y m [R m w0 (namely in the interior of the space of complexes). This can be shown by using (24) to compute its Jacobian with respect to y m , so that: with s T being of the form: where diag(v) and diag {1 (v) are the diagonal and inverse diagonal matrices operating over the vector v, respectively. The Jacobian D y m W is full rank since B is full rank and s[R m|m invertible. Continuity of W (y m ; c 0 )~0 then follows from the implicit function theorem (see Appendix S1). In the space of complexes, Condition (35) from Proposition 2 should be checked on the matrix:

Interval Based Search
In order to find the regions in the parameter space fulfilling the multistationarity condition we formulate a so called continuous constraint satisfaction problem (CSP) [35,36] and solve it numerically by using interval methods [37]. Methods based on interval analysis allow mathematical operations to be carried out over real intervals instead of real numbers, and thus to represent a continuum of solutions to a given CSP by a finite number of intervals or boxes, the union of which encloses the solution set.
Let z[R p be the space where a set of constraints is defined. A domain in that space is constructed by interval variables z i (for i~1, . . . ,p) defined on closed real intervals z i l ,z i u ½ , where subindexes l and u stand for lower and upper bounds, respectively. For p variables, the cartesian product of einterval domains z~z 1 ||z p is called a box.
Following [36], a constraint satisfaction problem consists in finding an interval domain where a set of constraints G~g 1 , . . . ,g q hold. This can be formally stated as: where constraints involve nonlinear analytic expressions and the symbol stands for either equality or inequality constraints, that is to say [f~,ƒ, §g. A solution of a CSP is an element of the search space which fulfills all the equalities and inequalities simultaneously. In our case, the search space involves the elements of the rate constant vector K and the independent variables which characterize the equilibrium manifold (which equals l~m{s, the manifold dimension). In this way, for K[R r w0 , the number of variables is p~rzm{s.
Regarding constraints, the equality ones correspond with the equilibrium manifold plus condition (35) which adds up to n{'z1 equations. On the other hand, inequality constraints must be imposed on the dependent variables that describe the equilibrium manifold, to ensure positivity of variables representing concentration and to search on non-zero a values.
Characterizing multistationarity regimes boils down to find the regions in the parameter-state space composed of the rate constants and independent variables which fulfill the constraints, i.e. computing all real feasible solutions of the corresponding CSP.
Solutions in interval methods are approximated by subpavings which consist of unions of boxes. Formally, a subpaving S for a given region R of the search space is defined as the union of nonoverlapping boxes approximating R. If we construct subpavings S and S such that:

S5R5S,
the region R is bracketed between inner and outer approximations. The outer approximation is reliable [38], since it is guaranteed that the solution region is contained within S.
Simple algorithms as SIVIA (Set Inverter Via Interval Analysis) proposed in [37] can be used to compute inner and outer subpavings by successive bisections and selections. Implementing interval algorithms requires environments supporting interval arithmetics such us the free available software package INTLAB, which provides an interactive environment within Matlab [39,40]. Efficiency can be gained by branch-and-prune algorithms where a set of boxes that contains all the solutions of the CSP is computed, and then each box is reduced and split. The free software package REALPAVER [36] provides a modeling language and a generic branch and prune algorithm which combine different splitting strategies and pruning techniques.

Example: Multistationary Regimes in an Autocatalytic Network
As a proof of concept, we make use of a simple autocatalytic network introduced by Edelstein [41], and already used as a case study for the bifurcation analysis in the context of biochemical systems [42]. Bistability has been shown for very simple autocatalytic systems like phosphorylation-dephosphorylation cycles with autocatalytic kinase [43].
The Edelstein network involves three species fA,B,Cg, distributed into five complexes, fA,B,C,AzB,2Ag. The C-graph for the Edelstein network is depicted in Figure 2, where the following index sets indicate complex interconnections: I 1~f 5g, I 2~f 3g, I 3~f 2,4g, I 4~f 3g, I 5~f 1g. The reaction steps that correspond with the set of edges in the graph are presented in Table 1. The network has two different linkage classes: L 1 containing the complexes (C 1 ,C 5 ), and L 2 containing the complexes (C 2 ,C 3 ,C 4 ). Their corresponding vectors L 1 and L 2 (see (5)) read: From the C-graph structure it follows that the network is weakly reversible since each linkage class is weakly reversible (in fact the network is reversible).  Let us consider the set of kinetic constants given in Table 3. For b~0:6451 they satisfy the condition described by (42). The resulting values for c 1 and c 2 are: c 1~2 :7416, c 2~1 2:4230: In Figure 3, the equilibrium manifold for these values of the kinetic constants is depicted, together with the reaction polyhedron corresponding to ½B 0 z½C 0~3 0. The manifold is one dimensional, and intersects the reaction polyhedron in three points, corresponding to three different equilibria. The points fulfilling the rank deficiency condition, corresponding to a~{15:787 and a~{9:079 are also indicated. As it can be deduced from the figure, three steady states will exist for a range of the sum of initial concentrations ½B 0 z½C 0 . In fact, performing a continuation of the curve of equilibria by varying the values of ½B 0 z½C 0 we obtain the curve shown in Figure 4, where two limit points or saddle node bifurcations appear for ½B 0 z½C 0~2 9:7768 and ½B 0 z½C 0~3 0:6949. Within these values, corresponding to different positions of the reaction polyhedron, three steady states will exist. Note that the points fulfilling the rank deficiency condition indicated in Figure 3 correspond precisely with the bifurcation points in Figure 4.
Alternatively, interval methods can be employed to search for the condition given by (35). Here it is important to remark that the method allows searching for parameters and/or steady state values of species concentrations within multistationary regimes, provided some other parameter values and/or steady state concentrations fixed. To illustrate this, let us assume we are interested in the ranges of parameters allowing for multistationarity, and the steady state values of c 1 , for a given steady state concentration of the species B (i.e. ½B~c 2 ). The variables in the constraint satisfaction problem are thus the free kinetic parameters and the steady state concentration of the species A.
A three dimensional plot of the result is given in Figure 5 for c 2~1 2:4230, where the x and y axis represent two of the kinetic parameters (k 23 and k 51 ), and the z axis represents a function of the steady state concentration of the species A in the steady state. We represent b in the z axis in order to facilitate the comparison with the results derived analytically. The corresponding values of c 1 can be computed as c 1~b r 51 =2.  . Equilibrium curve for the Edelstein network. The curve is obtained by varying ½B 0 z½C 0 , using the software Cl Matcont [46]. The kinetic parameters are kept fixed with the values shown in Table 3. doi:10.1371/journal.pone.0039194.g004  Table 3. Stars indicate steady states of the system for ½B 0 z½C 0~3 0. Dots indicate those points in the locus of equilibria where the rank deficiency condition is fulfilled. doi:10.1371/journal.pone.0039194.g003

Discussion
In this work we present a method to compute the regimes of multiple steady states in biochemical reaction networks, i.e. the regions in the parameter space, or in the state-parameter space of network models leading to multistationarity.
The main result of the paper consists of a sufficient condition for multistationarity, demonstrated to be valid for weakly reversible networks of arbitrary dimension and deficiency.
The idea behind is based on the fact that one steady state is an intersection between the locus of equilibria, or equilibrium manifold, and the reaction simplex, or stoichiometric compatibility class. In search for the existence of multiple steady states in reaction networks, we explored the requirements for multiple intersections between these two varieties.
In a previous work of the authors [29] it was shown that, for networks with d~l, (or equivalently with m~n{'), and fixed a kinetic rate vector K, the equilibrium manifold could be continued by the variation of the deficiency parameters of the network, and its qualitative behaviour evaluated through the derivative of the manifold with respect to the deficiency parameters. After partitioning of the parameter space in regions of different qualitative behaviour, global optimization algorithms [44] were used to check, within every region, whether multiple intersections between the manifold and the simplex were possible.
Here we use and extend this geometric insight to state a general condition of multistationarity for networks where the deficiency might be lower, equal or greater than the dimension of the equilibrium manifold. In this regard, [29] deals with a particular case of the general condition presented (and proved) here.
The evaluation of the multistationarity condition presented boils down to check the rank of a matrix which depends on the kinetic parameters, the concentrations, and the deficiency parameters. This matrix is systematically derived as indicated in the Analysis section, from the equations of the manifold and the mass conservation laws. In order to check the condition through the (state-)parameter space, we reformulate the problem as a constraint satisfaction one to be solved by interval methods.
For the purpose of characterizing the (state-)parameter space in terms of the capability for multiple steady states, methods based on interval arithmetics presents several advantages over classical global optimization methods [44]. On the one hand, they allow identifying regions, and not sets of discrete points, in the (state-)parameter space. In this way, there is no need to partition a priori the (state-)parameter space in regions with different qualitative behaviour. On the other hand, they ensure reliability of the solution [38], i.e. they guarantee that all the multistationarity regimes are enclosed by the solution regions.
As it has been commented in the Introduction, the CRNT provides particularly strong results to rule out multiple steady states based on the network structure irrespective of the network parameters. In this regard, this paper extends the results of the CRNT by introducing the kinetic parameters into the picture, giving a general condition for the appearance of multiple steady states.
In terms of applicability and performance, the method presented here is valid independently of the value of the network deficiency and dimension. The efficiency of the search will depend on the computational cost of the algorithm which increases in high dimensional spaces. In case we are only interested in finding points in the state-parameter space leading to multistationarity (for example if the goal is to decide only whether the network can exhibit multiple steady states or not), global optimization methods [44,45] would perform much faster.
Once the condition for multistationarity is derived, and depending on the particular scenario we want to explore, the search can be performed considering some of the parameters and/ or states to be fixed. In this way we obtain bifurcation diagrams in the desired low dimensional projections of the state-parameter spaces much more efficiently than using standard continuation techniques. For example, if we are interested in computing all the parameter sets giving multiple steady states for a particular value of the equilibrium concentration of one the species, we can perform the search keeping this concentration fixed.
To end up, we would like to stress the applicability of our method to standard problems appearing in systems and synthetic biology. The method can be particularly convenient, for example, to evaluate the robustness of the switching property for a given network. As remarked by [13], building a synthetic switch does not only require to start from a network structure allowing bistability, but such bistability must be sufficiently robust, meaning that the range of parameters leading to bistability is wide enough to enable a successful practical implementation. In this context, our method will efficiently characterize the regions of interest within the stateparameter space. Figure S1 Representation of the tangent and secant hyperplanes associated to a nonlinear manifold, in this case a curve, at x Ã ; e 1 and e 2 are unit vectors centered at x Ã ; y 1 and y 2 are two different points of the manifold.

Acknowledgments
We thank the anonymous reviewers for their helpful comments and suggestions.   Table 3. doi:10.1371/journal.pone.0039194.g005