Mathematical modelling has become an established tool for studying the dynamics of biological systems. Current applications range from building models that reproduce quantitative data to identifying systems with predefined qualitative features, such as switching behaviour, bistability or oscillations. Mathematically, the latter question amounts to identifying parameter values associated with a given qualitative feature. We introduce a procedure to partition the parameter space of a parameterized system of ordinary differential equations into regions for which the system has a unique or multiple equilibria. The procedure is based on the computation of the Brouwer degree, and it creates a multivariate polynomial with parameter depending coefficients. The signs of the coefficients determine parameter regions with and without multistationarity. A particular strength of the procedure is the avoidance of numerical analysis and parameter sampling. The procedure consists of a number of steps. Each of these steps might be addressed algorithmically using various computer programs and available software, or manually. We demonstrate our procedure on several models of gene transcription and cell signalling, and show that in many cases we obtain a complete partitioning of the parameter space with respect to multistationarity.
Mathematical models in the form of parameterized systems of ordinary differential equations (ODEs) are valuable tools in biology. Often, qualitative properties of the ODEs are associated with macroscopic biological properties and biological functions [1–4]. It is therefore important that we are able to analyse mathematical models with respect to their qualitative features and to understand when these properties arise in models. With the growing adaptation of differential equations in biology, an automated screening of ODE models for parameter dependent properties and discrimination of parameter regions with different properties would be a very useful tool for biology, and perhaps even more for synthetic biology [5]. Even though it is currently not conceivable how and if this task can be efficiently formalized, we view the procedure presented here as a first step in this direction. Multistationarity, that is, the capacity of the system to rest in different positive equilibria depending on the initial state of the system, is an important qualitative property. Biologically, multistationarity is linked to cellular decision making and ‘memory’-related on/off responses to graded input [2–4]. It has been suggested that different stable equilibria of a cell represent different cell types [6,7]. Whole-cell modelling provides an opportunity to understand the number and type of the stable equilibria of the cell and could potentially give insight into the different cell types that a particular cell can differentiate into and transition between. Currently, this is an important open question in biology [8]. Moreover, the existence of multiple equilibria is often a design objective in synthetic biology [9,10]. Various mathematical methods, developed in the context of reaction network theory, can be applied to decide whether multistationarity exists for some parameter values or not at all, or to pinpoint specific values for which it does occur [11–20]. Some of these methods are freely available as software tools [21,22]. It is a hard mathematical problem to delimit parameter regions for which multistationarity occurs. Often it is solved by numerical investigations and parameter sampling, guided by biological intuition or by case-by-case mathematical approaches. A general approach, in part numerical, is based on a certain bifurcation condition [18,19,23,24]. Alternatively, for polynomial ODEs, a decomposition of the parameter space into regions with different numbers of equilibria could be achieved by Cylindrical Algebraic Decomposition (a version of quantifier elimination) [25]. This method, however, scales very poorly and is thus only of limited help in biology, where models tend to be large in terms of the number of variables and parameters. Here we present two new theoretical results pertaining to multistationarity (Theorem 1 and Corollary 2). The results are in the context of reaction network theory and generalize ideas in [26,27]. We consider a parameterized ODE system defined by a reaction network and compute a single polynomial in the species concentrations with coefficients depending on the parameters of the system. The theoretical results relate the capacity for multiple equilibria or a single equilibrium to the signs of the polynomial as a function of the parameters and the variables (concentrations). The theoretical results apply to dissipative reaction networks (networks for which all trajectories eventually remain in a compact set) without boundary equilibria in stoichiometric compatibility classes with non-empty interior. These conditions are met in many reaction network models of molecular systems. We show by example that the results allow us to identify regions of the parameter space for which multiple equilibria exist and regions for which only one equilibrium exists. Subsequently this leads to the formulation of a general procedure for detecting regions of mono- and multistationarity. The procedure verifies the conditions of the theoretical results and further, calculates the before-mentioned polynomial. A key ingredient is the existence of a positive parameterization of the set of positive equilibria. Such a parameterization is known to exist for many classes of reaction networks, for example, systems with toric steady states [14] and post-translational modification systems [28,29]. The conditions of the procedure might be verified manually or algorithmically according to computational criteria. The algorithmic criteria are, however, only sufficient for the conditions to hold. For example, a basic condition is that of dissipativity. To our knowledge there is not a sufficient and necessary computational criterion for dissipativity, but several sufficient ones. If these fail, then the reaction network might still be dissipative, which might be verified by other means. By collecting the algorithmic criteria, the procedure can be formulated as a fully automated procedure (an algorithm) that partitions the parameter space without any manual intervention. The algorithm might however terminate indecisively if some of the criteria are not met. Table 1 shows two examples of reaction network motifs that occur frequently in intracellular signalling: a two-site protein modification by a kinase–phosphatase pair and a one-site modification of two proteins by the same kinase–phosphatase pair. These reaction networks are in the domain of the automated procedure and conditions for mono- and multistationarity can be found without any manual intervention. The conditions discriminating between a unique and multiple equilibria highlight a delicate relationship between the catalytic and Michaelis-Menten constants of the kinase and the phosphatase with the modified protein as a substrate (the kc- and kM-values). If the condition for multiple equilibria is met, then multiple equilibria occur provided the total concentrations of kinase, phosphatase and substrate are in suitable ranges (values thereof can be computed as part of the procedure). The paper has three main sections: a theoretical section, a section about the procedure and an application section. We close the paper with two brief sections discussing computational limitations, related work and future directions. In the theoretical section we first introduce notation and mathematical background material. We then give the theorem and the corollary that links the number of equilibria to the sign of the determinant of the Jacobian of a certain function, which is derived from the ODE system associated with a reaction network. In the second section we state the procedure, derive the algorithm and comment on the feasibility and verifiability of the conditions. Finally, in the application section we apply the procedure to several examples. The S1 File has six sections. All proofs are relegated to §1–4 together with background material. In §5 we elaborate further on how the conditions of the procedure/algorithm can be verified. In §6 we provide details of the algorithmic analysis of the examples in Table 1. Also we include a further monostationary example for illustration of the algorithm.
In this part of the manuscript we present the theoretical results. We start by introducing the basic formalism of reaction networks. Theorem 1, Corollary 1 and 2 below apply to dissipative networks without boundary equilibria and concern the (non) existence of multiple equilibria in some stoichiometric compatibility class. Corollary 2 assumes the existence of a positive parameterization of the set of positive equilibria. Before stating the results these five concepts are formally defined. In the previous subsection we applied Corollary 1 and Corollary 2 to the running example by going through a number of steps corresponding to the conditions of the statements and the calculation of the determinant. In this section we outline the steps formally. Afterwards we discuss the steps and how they can be verified either manually or algorithmically, that is, without user intervention. Finally we devise an algorithm to conclude uniqueness of equilibria or to find regions in the parameter space where multistationarity occurs. We conclude this section with some extra examples that follow the steps of the procedure. We assume the reaction rate functions v (x) depend on some parameters κ. The reaction rate functions are further assumed to be polynomials (as for mass-action kinetics) or quotients of polynomials (as for Michaelis-Menten and Hill kinetics with integer exponents). The input to the procedure is v (x) and N (the stoichiometric matrix) and the output is parameter regions for which the network admits multistationarity or uniqueness of equilibria. Procedure (Identification of parameter regions for multistationarity) Input: N and v (x) depending on κ. 1. Find f (x), a row reduced matrix W of size d × n such that WN = 0, and check that v (x) vanishes in the absence of one of the reactant species, that is, check that it satisfies Eq (2). 2. Check that the network is dissipative. 3. Check for boundary equilibria in P c for P c + ≠ ∅ and c ∈ R d. 4. Construct φc (x), M (x) and compute det (M (x) ). 5. Analyze the sign of det (M (x) ). Find conditions on the parameters κ such that sign (det (M (x) ) ) = (−1) s for all x ∈ R > 0 n, in which case Corollary 1 holds. If Corollary 1 does not hold for all κ, continue to the next step. 6. Obtain an algebraic parameterization Φ (x ^) of the set of positive equilibria for all κ, as in Eq (7), such that the coefficients of the numerator and the denominator of each Φ i (x ^) possibly depend on κ. Compute a (x ^) = det (M (Φ (x ^) ) ). By hypothesis, a (x ^) can be written as the quotient of two polynomials in x ^ with coefficients depending on κ, whose denominator takes positive values. 7. Analyze the sign of the numerator of a (x ^). 7a. Identify coefficients with sign (−1) s+1 and coefficients that can have different signs depending on the parameters. 7b. Use the terms corresponding to identified coefficients to construct parameter inequalities such that, whenever these inequalities hold, one has either sign (a (x ^) ) = (- 1) s for all x ^ ∈ R > 0 m or sign (a (x ^) ) = (- 1) s + 1 for at least one x ^ ∈ R > 0 m, in which case either Corollary 2 (A) or (B) holds. There is no guarantee that all steps of the procedure can be carried out successfully, let alone automatically. While step 1 and 4 usually are straightforward (only computational issues might arise for large networks), step 2,3, 5,6 and 7 might in particular require case specific approaches. However, there exist computationally feasible sufficient criteria that guarantee the conditions in each step can be checked efficiently. To illustrate several aspects of the algorithm we provide a detailed step-by-step analysis of a collection of examples. The computational complexity of some of the steps in the procedure are demanding. Some conditions can be checked using linear algebra and do not depend on parameter values, others depend on parameter values and require symbolic manipulations. In some situations, the calculation can be done for even large networks at the cost of time, while in other situations symbolic software (like Mathematica and Maple) have inherent limits to what it can process. We offer here a few remarks about computational strategies and time complexity. Dissipativity. There are efficient algorithms to check whether the network is conservative and strongly endotactic, using linear algebra or mixed-integer linear programming [21,47]. We are not aware of a systematic way to check if Proposition 1 is fulfilled or not. Finding the minimal siphons of a network requires in general exponential time and there might be exponentially many of these [59]. Different algorithms developed in Petri Net theory can be applied to find the minimal siphons; see for example [48,49,59] and references therein. The complexity of this computation can often be substantially reduced by removing so-called intermediates and catalysts from the network [50] (see §5. 1 in the S1 File for details). Finding all non-interacting and reactant-non-interacting sets requires in general exponential time. One strategy is the following. We first remove all species Si for which αij > 1 or βij > 1 for some reaction Rj (the latter constraint is omitted if we are looking for reactant-non-interacting sets only). Then we build non-interacting (reactant-non-interacting) sets by adding new species recursively until no more species can be added without having an interacting pair of species in the set. Calculation of the symbolic determinant of the matrix M (x), and hence also of a (x ^), often fails in our experience for networks with more than 20 variables on common laptops [60]. However, this clearly depends on the sparsity of the matrix M (x), that is, on the number and order of the reactions. Strategies to reduce the complexity of the computation by expanding the determinant along the non-symbolic rows (conservation relations) were inspected in [60]. Specialized software like Singular [61] and/or better hardware could probably push what is possible to something closer to 50 variables. At this size, however, what might best be called ‘cognitive limitations’ come into play: symbolic software typically has problems with collecting and simplifying terms ‘the right way’ if there are many variables and/or parameters. And if terms are not collected appropriately it might be difficult, if not impossible, to decide on the sign of the polynomial coefficients. Our approach is therefore best suited to systems of moderate size (say 20-30 variables). Furthermore, it is our experience that large non-linear models tend to be multistationary because of the many non-linear dependencies that typically are present [60]. Positive parameterizations: The worst case scenario involves checking ∑ i = 1 d (n i) different sets of variables, each with at most d variables. Finding the vertices of the Newton polytope can be done with existing symbolic software, for example Polymake [62] or Maple, as we demonstrate in the S1 File. We stress that it is always beneficial to guide the procedure/algorithm whenever possible in the sense that, if something is known for the network, there is no reason to go through many possibilities.
The main result of this paper, the procedure to identify parameter regions for unique and multiple equilibria, combines Brouwer degree theory and algebraic geometry. In particular, under the assumptions of Corollary 2, we show that there exist stoichiometric compatibility classes with at least two equilibria if, and only if, a certain multivariate polynomial can attain a specific sign. Discriminating regions of the parameter space where multistationarity occurs is a hard mathematical problem, theoretically addressable by computationally expensive means [25]. Our approach beautifully overcomes these difficulties by building on a simple idea, the computation of the Brouwer degree of a function related to a dissipative network. Additionally, not only closed-form expressions in the parameters are obtained, but, as illustrated in examples, these expressions are often interpretable in biochemical terms, providing an explanation of why multistationarity occurs. The procedure applies theoretically to any choice of algebraic reaction rate functions. However, in practice, the procedure works well with mass-action kinetics. For example, we have considered the two-site phosphorylation cycle depicted in the second row of Table 1, but now modelled with Michaelis-Menten kinetics instead of mass-action kinetics. This network is known to be multistationary [63], and the conditions to apply Corollary 1 and Corollary 2 are valid. However, a positive algebraic parameterization does not exist, and hence our approach cannot be used to find parameter conditions for multistationarity. However, Corollary 1 might be used with rational reaction rate functions for monostationary networks. This is the case for example for the one-site phosphorylation cycle S ⇌ S p with Michaelis-Menten kinetics [63]. This network has two species and rank one. The sign of det (M (x) ) is −1 for all parameter values and all x ∈ R > 0 2. By Corollary 1, the network admits exactly one positive equilibrium in every stoichiometric compatibility class P c with P c + ≠ ∅ for all parameter values. If a reaction network does not have any conservation relation, then the set of equilibria consists typically of a finite number of points. In this case an algebraic parameterization is an algebraic expression of the equilibria in terms of the parameters of the system. Since m = 0, then R m consists of a single point and it follows directly that there is a unique equilibrium. Such an expression rarely exists. Therefore the procedure applies mainly to reaction networks with conservation relations. In particular, this rules out reaction networks where each species is produced and degraded. Several natural questions remain outside the reach of our procedure. Firstly one would like to determine the particular stoichiometric compatibility classes for which there are multiple equilibria. As stated in Corollary 2, if sign (a (x ^) ) = (- 1) s + 1, then c: = W Φ (x ^) defines a stoichiometric compatibility class with multiple equilibria. However, this only establishes c indirectly through x ^. In some situations, it might be possible to find a positive parametrization that uses some of the conservation relations (ideally, all but one) and the stoichiometric compatibility classes with multiple/single equilibria would be determined up to a single parameter. Secondly, one could ask for parameter regions that differentiate between the precise number of equilibria (that is, 0,1, 2, …). This question should be seen in conjunction with the previous question: in typical examples, when there are two equilibria in a particular stoichiometric compatibility class, then there exists another class for which there are three. Hence the number of equilibria cannot be separated from the stoichiometric compatibility classes. A third question concerns the stability of the equilibria, which cannot be obtained from our procedure. It is, however, known that if the sign of the Jacobian evaluated at an equilibrium is (−1) s+1, then it is unstable [34]. This is in particular the case for an equilibrium fulfilling the condition in Corollary 2 (B). We have shown that for some reaction networks our procedure can be formulated as an algorithm. We consider therefore our research a step in the direction of providing ‘black box tools’ to analyse complex dynamical systems. Such tools would easily find their use in systems and synthetic biology, where it is commonplace to consider (many) competing models. A particular problem is to exclude models that cannot explain observed qualitative features, such as multistationarity.
We used Maple for the symbolic computations, such as finding det (M (x) ), the positive parameterizations, a (x ^) and the vertices of the Newton polytope.
Mathematical modelling has become an important tool in biology. As modelling requires separating the essential from the ordinary, there is never just one model but a collection thereof. To understand biology through modelling it is therefore crucial to be able to tell which of these models are capable of reproducing an observed behaviour and which are not. For example, to understand cellular decision making, models allowing multiple equilibria are studied and one asks which models allow for this behaviour. Here we describe a procedure that links the existence of a unique and of multiple equilibria to the sign of a single expression. We demonstrate the usefulness of the procedure by applying it to models of gene transcription and cellular signalling.
phosphorylation
chemical compounds
applied mathematics
reactants
organic compounds
dna transcription
simulation and modeling
algorithms
systems science
mathematics
chemical equilibrium
algebra
basic amino acids
amino acids
polynomials
physical chemistry
research and analysis methods
computer and information sciences
system stability
proteins
gene expression
chemistry
biochemistry
histidine
organic chemistry
post-translational modification
genetics
biology and life sciences
chemical reactions
physical sciences
Identifying parameter regions for multistationarity