Graph-theoretic analysis of multistationarity using degree theory
Biochemical mechanisms with mass action kinetics are often modeled by systems of polynomial differential equations (DE). Determining directly if the DE system has multiple equilibria (multistationarity) is difficult for realistic systems, since they are large, nonlinear and contain many unknown parameters. Mass action biochemical mechanisms can be represented by a directed bipartite graph with species and reaction nodes. Graph-theoretic methods can then be used to assess the potential of a given biochemical mechanism for multistationarity by identifying structures in the bipartite graph referred to as critical fragments. In this article we present a graph-theoretic method for conservative biochemical mechanisms characterized by bounded species concentrations, which makes the use of degree theory arguments possible. We illustrate the results with an example of a MAPK network.
Keywords. Biochemical mechanisms, mass-action kinetics, multistationarity, bipartite graph, MAPK network.
Biochemical mechanisms of chemical species and elementary reactions are often modeled by differential equations (DE) systems with the species concentrations as variables. Multistability, the existence of multiple stable positive equilibria (for some choice of parameter values) is ubiquitous in models of biochemical mechanisms, such as cell decision [15, 17]. And multistationarity, the existence of multiple positive equilibria is necessary for multistability or a biological switch, a term used in the biological literature.
The models in this work will be taken with mass action kinetics resulting in a polynomial right-hand side of the DE system. The DE system models of the biochemical mechanisms of interest are typically high-dimensional, nonlinear and contain many unknown parameters (rate constants and total concentrations). Thus determining parameter values such that multiple equilibria can be found by solving numerically large nonlinear polynomial systems is difficult, if not impossible. On the other hand solving a nonlinear polynomial system with unknown coefficients directly using methods from algebraic geometry has its limitations . Therefore other methods and approaches such as graph-theoretic are being developed to answer the question of the existence of multistationarity more easily.
A biochemical mechanism with mass action kinetics can be represented as a directed bipartite graph, which is a graph with two non-intersecting sets of nodes representing species and reactions, and directed edges starting at a species (reaction) node and ending at a reaction (species) node. Graph-theoretic methods can be used to identify structures referred to as critical fragments that are necessary for the existence of multistationarity .
Graph-theoretic methods have been used to determine the potential of various biochemical mechanisms for multistationarity [2, 7, 16, 18]. Many of these methods use the one-to-one correspondence between structures in the graph (fragments [16, 18] or cycle structures ) and the summands in the determinant of the Jacobian of the right-hand side of the DE system. However, many models have conservation relations with positive coefficients of the species concentrations. The existence of conservation relations results in a non-full rank Jacobian. This leads to considering a coefficient of the characteristic polynomial of the Jacobian different from the constant coefficient and its sign when studying multistationarity. Here we study conservative biochemical mechanisms where all species concentrations participate in at least one conservation relation. This means that all species concentrations are bounded from above and degree theory  can be used to study the number of equilibria of the DE model. So far degree theory has been used to study multistationarity in biochemical mechanism models, for example, in [5, 8, 9] and graph-theoretic methods have been developed in [2, 7, 16, 18]. Here we combine both approaches to develop a graph-theoretic method for multistationarity in a conservative biochemical mechanism DE model.
This article is organized as follows. Sec. 2 provides an introduction to conservative biochemical mechanisms with mass action kinetics and their properties. In Sec. 3 we discuss consequences of the well-known fact that solutions of the DE systems under study are confined to affine linear subspaces defined by these conservation relations. In Sec. 4 the Jacobian of the original DE system, its parametrization and the determinant of the Jacobian on the level sets is given. In Sec. 5 the degree of a nonlinear function and some of its properties related to a DE system’s right-hand side on a given level set is introduced. In Sec. 6 the bipartite graph of a biochemical mechanism with mass action kinetics is introduced. The main result in Sec. 7 (Theorem 4 and Corollary 5) gives a necessary condition for multistationarity for conservative biochemical mechanism models. An example of a MAPK network model studied for multistationarity in  is presented in the same section.
A (bio)chemical mechanism with species , , and elementary reactions is represented as
where , are the rate constants. The constants and are small integers called stoichiometric coefficients that account for the number of molecules of species participating in the elementary reaction in (1). An example of a chemical mechanism is given below
A true reaction is a reaction different from an inflow reaction or an outflow reaction . An autocatalytic reaction is a reaction of the form where . We assume that every species in (1) is consumed and produced in at least one true non-autocatalytic elementary reaction.
The chemical mechanism in (2) satisfies the above assumption.
We will denote by the vector of concentrations of species and by the vector of rate constants. If for , () for all we will write . Since each as a concentration, we have . Similarly since each rate constant .
If mass action kinetics is used for the mechanism (1), then the corresponding rate functions are
The vector of rate functions will be denoted as where .
The differential equation (DE) model of a mass-action mechanism such as (1) can be written as
are the entries of the stoichiometric matrix with dimension and is the vector of rate functions (3). In what follows we will at times interpret the right hand side of (4) as a function of the rate constants and the species concentrations and at other times as a function of the reaction rates and the concentrations , depending on the situation.
For the system (2) we obtain the stoichiometric matrix
The equations of the system (4) can be written componentwise as
The model equations of the mechanism (2) are given below
Initially , and we will denote a solution of (4) with initial condition as .
3 The dynamics on the level sets
Let the stoichiometric matrix have rank . Suppose that at least one solution of the system
exists. Then we have mass conservation laws
If for all in at least one solution , then by (8) it follows that all species concentrations are conserved, i.e., for all where . A biochemical mechanism (1) with mass-conserved species concentrations will be called conservative biochemical mechanism .
Here we study only biochemical mechanisms where the system
has a solution. As outlined above, in this situation all species concentrations are bounded.
The rank of the stoichiometry matrix of the system (6) equals . The left kernel of N is spanned by the vector . Hence the network is conservative and there exists one conservation relation,
We can rewrite (8) in a matrix form as
where is a full rank matrix whose columns span .
In what follows we interpret the entries of the dimensional vector as additional parameters and study the dynamics of the system (4) on the level sets
This is motivated by the observation that the sets are invariant under the dynamics of (4).
Lemma 1 ( convex, compact and forward invariant).
The set is convex, compact and forward invariant.
To study the dynamics of system (4) on invariant sets we let be the matrix of full column rank whose columns are an orthonormal basis of and we let the matrix be the matrix of full column rank whose columns are the orthonormal basis of . Then the linear transformation
sends to an element and to an element :
Note that and are unique for given , (as and are complementary subspaces). Since by assumption, and are orthonormal we recover
We further note that by construction and hence all elements are sent to the same element :
That is, is constant, reflecting the invariance of . We introduce the abbreviation:
where in complete analogy to above we interpret as a parameter vector. Then every identifies an -dimensional dynamical system
[ convex, compact and forward invariant] The set is convex, compact and forward invariant.
For a set we will denote its interior by and its boundary by .
From hereon we will assume level sets do not contain equilibria with zero coordinates, referred to as boundary equilibria.
Here we study only DE models of biochemical mechanisms such that the level set does not contain any boundary equilibria, that is,
4 The Jacobian parametrized at and its projection on space
The Jacobian matrix of (4) has entries
Recall (3), then the Jacobian can be written also as
Note that if the concentrations and the rate functions are evaluated at a positive equilibrium, they are positive and can be used as parameters in (19).
For example, the Jacobian matrix of the right-hand side of the system (6) is
For the remainder of this contribution we make the following assumption.
We assume that .
Since the rank of is , it follows under the above assumption that and that the characteristic polynomial of the Jacobian (19) is
where the coefficients , ,…, are computed as the sum of all principal minors of order of the negative Jacobian .
It is easy to verify, that the Jacobian with respect to of from (15) is given by
where we have suppressed the , dependence of .
The relation between the coefficient and the determinant of the negative Jacobian of the right-hand side of the reduced system is considered in the next lemma. A special case of this lemma with is available in .
The following equivalence is true
Let , …, be the eigenvalues of and hence the roots of the characteristic polynomial (21), where we assume that . By (21), there exist trivial eigenvalues that are identically zero for all values of and non-trivial eigenvalues that are nonzero for some values of . For simplicity we assume that , , …, are nontrivial and , are trivial.
First we show that is the product of the nontrivial eigenvalues
If we write the polynomial from the characteristic polynomial (21) in factored form and substitute we obtain
Next we show that , which will prove the claim . For this purpose we apply the orthonormal transformation to . The transformed matrix has a block form
Since is orthonormal, and have the same eigenvalues. Both matrices and have trivial eigenvalues. Thus
If is written in factored form and we let we obtain . Thus . ∎
5 The degree of
Let be an open and bounded set. The closure of will be denoted by and the boundary of by . Thus, is a compact set.
Let be a smooth function, where using the usual notation we write . We denote the Jacobian matrix of by
and its determinant by . A point is a regular point for if . A point is called a regular value if all such that are regular.
Next we define the (topological or Brouwer) degree of , denoted by . In the next definition we use the sign function .
[(topological) degree] If and is a regular value, the degree of F is defined by
The next lemma is similar to Lemma 5.4 available in the Supporting Information of .
We obtain the following corollaries on the degree of where . Similar corollaries for the special case of are available in .
First we need the following theorem on the homotopy invariance of the degree .
Let be a bounded and open set. Let , be a continuously varying set of functions such that does not have any zeroes on the boundary of for all . Then is constant for all .
Since is smooth on , therefore is
smooth on . Let be fixed but arbitrary so that .
We have , where is the interior of and is the boundary of . By Lemma 2, is bounded. We follow the proof of [9, Lemma 2]. Let be an arbitrary point and consider the function
By Definition 1 it follows that
Next we show that and are homotopic. We define the following homotopy
where . Therefore is continuous on , and . To apply Theorem 1 we need to show that for all and for all . The latter is true if since and if by Corollary 1. Suppose it is not true if , then there exists and such that
By the convexity of it follows that points strictly inwards at . Thus at points strictly outwards. This is a contradiction since is forward invariant by Lemma 2. Thus the claim in (28) follows by equation (29) since the degree is homotopy invariant by Theorem 1. ∎
Let and be given and note that . Assume that the boundary does not contain any equilibria of (16). If for all , then the equation
has a unique solution.
If all solutions of , are regular, then the number of solutions in is odd.
6 The bipartite digraph of a biochemical mechanism
For the convenience of the reader, in this section we present definitions regarding the bipartite digraph of a biochemical mechanism (1) [18, 16]. To illustrate the definitions we will continue to use as an example the mechanism (2).
A directed bipartite graph (bipartite digraph) has a node set that consists of two disjoint subsets, and , and each of its directed edges (arcs) has one end in and the other in .
The bipartite digraph of a biochemical reaction network (1) is defined as follows. The nodes are separated into two sets, one for the chemical species and one for the elementary reactions . We draw an arc from to if and only if species is a reactant in reaction , i.e., if the stoichiometric coefficient in (1). Similarly, we draw an arc from to if and only if is a product in reaction , i.e., if the stoichiometric coefficient in (1). Therefore the set of arcs consists of arcs such as and . Hence the bipartite digraph can be defined as where is the set of nodes and is the set of arcs. If an arc is not weighted explicitly, we assume that its weight equals . The corresponding bipartite digraph of the mechanism (2) is shown in Figure 1.
The element is an edge if , i.e., if species is a reactant in reaction . The weight of an edge is defined as
For example, the edge corresponding to the arc in Figure 1 has weight .
If , then the arcs and form a positive path that corresponds to the production of from in a reaction . The weight of the positive path is defined as . For example, the positive path in Figure 1 has weight .
If , then the arcs and form a negative path that corresponds to and interacting as reactants in reaction . The weight of the negative path is defined as . Note that the negative paths and are considered to be different since they start at a different species node. For example, both and in Figure 1 are negative paths with weight . We note that the direction of the arcs is followed in the positive paths but not in the negative paths.
A cycle of is a sequence of distinct paths with the last species node of each path being the same as the first species node of the next path , ,, , . A cycle will be denoted by , where the number of species nodes defines its order. The set of species nodes in a cycle is distinct, but there may be a repetition among the reaction nodes. This is because negative paths containing the same nodes are considered different depending on the starting species node. For example, in Figure 1 is a cycle formed by the two negative paths and .
A cycle is positive if it contains an even number of negative paths and negative if it contains an odd number of negative paths. The sign of a cycle can also be determined by the cycle weight which is a product of all corresponding weights of negative and positive paths of
For example, in Figure 2 is a positive cycle of order with weight .
A subgraph of consists of edges or cycles , , where each species is the beginning of only one edge, or one path participating in a cycle. The number of species nodes in a subgraph is defined as its order. The subgraph weight is defined as
Since more than one path can exist between species nodes via different reaction nodes in a bipartite digraph, the number of subgraphs through the same node sets may be greater than one. The set of all subgraphs of order with the same species nodes and reaction nodes sets is called a fragment of order and is denoted by . For a fragment we define the number
as the fragment weight. If , then is defined as a critical fragment.
For example, the fragment shown in Figure 4 together with its two subgraphs and . The first subgraph is a positive cycle and thus it has a negative weight. Therefore is a critical fragment since