Graph-theoretic analysis of multistationarity using degree theory

Graph-theoretic analysis of multistationarity using degree theory

Abstract

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.

1 Introduction

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 [14]. 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 [16].

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 [7]) 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 [10] 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 [3] is presented in the same section.

2 Preliminaries

A (bio)chemical mechanism with species , , and elementary reactions is represented as

 n∑i=1αijAikj−→n∑i=1βijAi,j=1…m, (1)

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

 A2+A3k1−→2A1,A3k2−→A1,A1k3−→A3,A2k4−→A1A1k5−→A2. (2)
Assumption 1.

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

 vj(k,x)=kjxα1j1…xαnjn,j=1,…,m. (3)

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

 ˙x(t)=Nv(k,x)=f(k,x)=f(v,x) (4)

where

 Nij=βij−αij

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

 N=⎡⎢⎣−2−1−1−1−1−1−0−0−1−1−1−1−1−0−0⎤⎥⎦ and the % rate functions v(k,x)=(k1x2x3,k2x3,k3x1,k4x2,k5x1)T .

The equations of the system (4) can be written componentwise as

 ˙xi(t)=m∑j=1Nijvj(k,x),i=1,…,n. (5)

The model equations of the mechanism (2) are given below

 ˙x1 =2k1x2x3+k2x3−k3x1+k4x2−k5x1=2v1+v2−v3+v4−v5, (6) ˙x2 =−k1x2x3−k4x2+k5x1=−v1−v4+v5, ˙x3 =−k1x2x3−k2x3+k3x1=−v1−v2+v3.

Initially , and we will denote a solution of (4) with initial condition as .

3 The dynamics on the level sets ωc0

Let the stoichiometric matrix have rank . Suppose that at least one solution of the system

 n∑i=1λiNij=0,j=1,…,m (7)

exists. Then we have mass conservation laws

 n∑i=1λixi=n∑i=1λixi(0). (8)

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 [8].

Assumption 2.

Here we study only biochemical mechanisms where the system

 λTN=0,λ>0

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

 WTx(t,x0)≡WTx0=c0 (9)

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

 ωc0={x≥0|WTx=c0}. (10)

This is motivated by the observation that the sets are invariant under the dynamics of (4).

Lemma 1 (ωc0 convex, compact and forward invariant).

The set is convex, compact and forward invariant.

The proof of Lemma 1 is available in [5].

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

 x→(STx,ZTx) (11)

sends to an element and to an element :

 ξ:=STx and η:=ZTx. (12)

Note that and are unique for given , (as and are complementary subspaces). Since by assumption, and are orthonormal we recover

 x≡x(ξ,η)=Sξ+Zη. (13)

We further note that by construction and hence all elements are sent to the same element :

 x1,x2∈ωc0⇒ZTx1=ZTx2=:η0,∀x1,x2∈ωc0.

We now apply the linear transformation (11) to the system (4) to obtain:

 ˙ξ =ST˙x=STNv(k,x(ξ,η)) (14a) ˙η =ZT˙x=ZTNv(k,x(ξ,η))≡0. (14b)

That is, is constant, reflecting the invariance of . We introduce the abbreviation:

 gη(k,ξ):=STNv(k,x(ξ,η)), (15)

where in complete analogy to above we interpret as a parameter vector. Then every identifies an -dimensional dynamical system

 ˙ξ=gη0(k,ξ). (16)

Now studying the system (4) restricted to a level set is equivalent to studying the system (16) with for some .

Solutions of (16) give rise to solutions of (4)

 x(t,x0)=Sξ(t,ξ0)+Zη0.

Since a solution of the system (4), for all , it follows that the corresponding solution of (16) remains in the set

 Ωη0={ξ∈IRr|Sξ≥−Zη0}. (17)

The set has similar properties as the set . We have the following lemma for which will be used in Corollary 2. The proof is available in [5].

Lemma 2.

[ convex, compact and forward invariant] The set is convex, compact and forward invariant.

The following lemma compares the number and type of equilibria of (4) in the set to the number and type of equilibria of (16) in the set .

For a set we will denote its interior by and its boundary by .

Lemma 3.

[Equilibrium points.]

• A positive point is an equilibrium of (4) in with , if and only if is an equilibrium of (16) for (where positivity of entails ).

• The number of equilibria in of (4) equals the number of equilibria of (16) in .

• The boundary of the set contains an equilibrium point of (4) if and only if the boundary of contains an equilibrium point of (16).

Proof.
• This follows since , where corresponds to such that .

• This follows from the fact that and are full rank matrices, and thus the correspondence between and is one-to-one.

• We prove the contrapositive. An equilibrium of (4) is not on if and only if an equilibrium of (16) is such that .

From hereon we will assume level sets do not contain equilibria with zero coordinates, referred to as boundary equilibria.

Assumption 3.

Here we study only DE models of biochemical mechanisms such that the level set does not contain any boundary equilibria, that is,

 if x∈∂ωc0⇒Nv(k,x)≠0.

The following corollary follows by Lemma 3 and Assumption 3.

Corollary 1.

Under Assumption 3 the boundary of does not contain an equilibrium of (16).
The number of equilibria of (4) in equals the number of equilibria of (16) in .

Remark 1.

In essence, in Sec. 5 we will study the number of equilibria of the reduced system (16) in . By Corollary 1 we will obtain the corresponding result on the number of equilibria of the system (4) in .

4 The Jacobian J parametrized at (v,x) and its projection on (ξ,η) space

The Jacobian matrix of (4) has entries

 Jil(k,x)=m∑j=1Nijαljkjxα1j1…xαlj−1l…xαnjn. (18)

Recall (3), then the Jacobian can be written also as

 Jil(k,x)=Jil(v,x)=m∑j=1Nijαljvjxl. (19)

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

 J(v,x)=⎛⎜ ⎜ ⎜ ⎜⎝−v3x1−v5x12v1x2+v4x22v1x3+v2x3v5x1−v1x2−v4x2−v1x3v3x1−v1x2−v1x3−v2x3⎞⎟ ⎟ ⎟ ⎟⎠. (20)

For the remainder of this contribution we make the following assumption.

Assumption 4.

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

 det(λI−J(v,x))=λn−r(λr+a1λr−1+…+ar−1λ+ar)=λn−rq(λ), (21)

where the coefficients , ,…, are computed as the sum of all principal minors of order of the negative Jacobian [11].

The coefficients of (21) are rational functions in and by (19). For example, the non-zero coefficients of the characteristic polynomial of the Jacobian (20) are

 a1(v,x) =v3+v5x1+v1+v4x2+v1+v2x3 (22) a2(v,x) =v1v3−v1v5+v3v4x1x2+−v1v3+v1v5+v2v5x1x3+v1v2+v1v4+v2v4x2x3. (23)

It is easy to verify, that the Jacobian with respect to of from (15) is given by

 Gη(v,ξ)=STJ(v,x)S, (24)

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 [5].

Lemma 4.

The following equivalence is true

 det(−Gη(v,ξ))≡det(−STJ(v,x)S)≡ar(v,x) .
Proof.

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

 ar=r∏i=1λi.

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

 ϕTJϕ=[STJSSTJZ00].

Since is orthonormal, and have the same eigenvalues. Both matrices and have trivial eigenvalues. Thus

 det(λI−ϕTJ(v,x)ϕ)=λn−rdet(λI−STJ(v,x)S)=λn−r~q(λ).

If is written in factored form and we let we obtain . Thus . ∎

5 The degree of gη(k,ξ)

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

 ~J(x)=[∂Fi∂xj] (25)

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 [10], denoted by . In the next definition we use the sign function .

Definition 1.

[(topological) degree] If and is a regular value, the degree of F is defined by

 deg(F)=deg(F,U,y)=∑F(x)=ysign(det(−~J(x))). (26)
Remark 2.

Note that we use in (26) in place of to avoid the case of the degree depending on the dimension of similarly to [13].

The sum in (26) is over all solutions of such that . If does not have solutions , then we set . Since we are interested in the equilibrium solutions of that satisfy , we will let in (26).

Next we study the degree of the function defined in (16). Recall that and that by Lemma 4, where is the Jacobian of the function given in (24).

The next lemma is similar to Lemma 5.4 available in the Supporting Information of [5].

Remark 3.

Note that in the lemma and corollaries below the assumption that does not contain any boundary equilibrium of the system (16) is automatically satisfied by Assumption 3 and Corollary 1.

Lemma 5.

Let be as in (15). Fix and assume that the boundary does not contain any equilibria of (16). If all equilibria are regular, then

 deg(gη,int(Ωη),0)=∑{ξ∈int(Ωη)|gη(ξ)=0}sign(ar(v,x(η,ξ))). (27)
Remark 4.

By Corollary 1 the equilibria of (16) are in . Therefore the degree of , given by (27) is well defined.

We obtain the following corollaries on the degree of where . Similar corollaries for the special case of are available in [5].

First we need the following theorem on the homotopy invariance of the degree [8].

Theorem 1.

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 .

Corollary 2.

Let be as in (15) and assume that the boundary does not contain any equilibria of (16).Then the the following holds true:

 deg(gη,int(Ωη),0)=1. (28)
Proof.

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

 G(ξ)=¯ξ−ξ.

By Definition 1 it follows that

 deg(G,int(Ωη),0)=1. (29)

Next we show that and are homotopic. We define the following homotopy

 H(ξ,s)=sgη(ξ)+(1−s)G(ξ)

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

 gη(~ξ)=−1−ssG(~ξ).

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. ∎

Corollary 3.

Let and be given and note that . Assume that the boundary does not contain any equilibria of (16). If for all , then the equation

 gη(v,ξ)=0,ξ∈Ωη

has a unique solution.

If all solutions of , are regular, then the number of solutions in is odd.

Proof.

Recall that by Lemma 4, . Suppose that for all . Then by Corollary 2, . Thus by (27), it follows that the number of equilibria in equals one.

If all equilibria are regular, then is either positive or negative at . Since by Corollary 2, it follows by (27) that the number of equilibria has to be odd. ∎

Now we turn to the system (4). In the next theorem we show that the system (4) has an interior equilibrium solution in any set where is fixed.

Theorem 2.

Let be as in (4) and recall Assumption 3. If

 ar(v,x)>0, for all v>0 and x>0,

then

 f(k,x)=0,x∈ωc0

has a unique solution for all .

Proof.

Pick and choose any . Compute , . By assumption

 ar(v,x(η0,ξ))>0,∀v>0 and ∀ξ∈Ωη0.

Thus, by Corollary 3, , has a unique solution . And by Lemma 3 (a)-(b), an equilibrium of , is unique. ∎

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 [12].

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

 KE=−α2kj. (30)

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

 KC=∏¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯[Ak,Bj,Ai]∈C(−αkjαij)∏[Ak,Bj,Ai]∈Cαkjβij. (31)

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

 Kg=(−1)c∏C∈gKC∏E∈g(−KE), (32)

where is the number of cycles in , is the cycle weight (31) and is the edges weights (30) of the cycles and edges in . For example, the subgraph with weight is shown in Figure 3 .

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

 KSk=∑g∈SkKg (33)

as the fragment weight. If