On the multistationarity of chemical reaction networks

TL;DR

This paper proposes a new conjecture outlining a necessary condition for chemical reaction networks to exhibit multistationarity, supported by examples and extending known results for monotonic kinetics.

Contribution

It introduces a novel conjecture on multistationarity conditions in chemical networks, expanding understanding beyond strictly monotonic kinetics.

Findings

The conjecture is supported by multiple illustrative examples.

It extends the known theorem for strictly monotonic kinetics.

Provides a new perspective on conditions for multistationarity.

Abstract

We present a new conjecture about a necessary condition that a (bio)chemical network has to satisfy for it to exhibit multistationarity. According to a Theorem of Feliu and Wiuf [27, 12], the conjecture is known for strictly monotonic kinetics. We give several examples illustrating our conjecture.

Figures

B, for suitable parameter values, this slightly modified Brusselator becomes a bistable system without losing the capacity to exhibit oscillations (Fig. 7D). Excitability linked to multistationarity is also observed (Fig. 7C). In Fig. 7D one can see that $k_{4}$ favors bistability while $k_{3}$ favors sustained oscillations.

B, for suitable parameter values, this slightly modified Brusselator becomes a bistable system without losing the capacity to exhibit oscillations (Fig. 7D). Excitability linked to multistationarity is also observed (Fig. 7C). In Fig. 7D one can see that $k_{4}$ favors bistability while $k_{3}$ favors sustained oscillations.

Equations

$$r:y⟶y^{\prime }$$

$$K_j:{\mathbb{R}}_{+}^n⟶\mathbb{R}$$

$$F:{\mathbb{R}}_{+}^n⟶{\mathbb{R}}^n$$

$$F(c)=\sum _{j=1}^m{K}_j(c)(y_j^{\prime }-y_j).$$

$$\dot{c}=F(c).$$

$$F:{\mathbb{R}}_{+}^n⟶{\mathbb{R}}^n$$

$$\tilde{F}(c)=({\omega }^1\cdot c,\dots ,{\omega }^d\cdot c,F_{d+1}(c),\dots ,F_n(c))$$

$$\tilde{F}:{\mathbb{R}}_{+}^n⟶{\mathbb{R}}^n.$$

$$\det (J_c(\tilde{F}))=\sum _{\genfrac{}{}{0ex}{}{I\subset \{1,\dots ,n\}}{\#I=n-d}}\det (J_c(F{)}_{I,I}).$$

$$d(c,K)=\det (J_c(\tilde{F})).$$

$${\omega }^i\cdot c=k_i,1\le i\le d,$$

$$c_i=k_i+u_i(c_{d+1},\dots ,c_n),$$

$${\phi }_i(c_{d+1},\dots ,c_n)=F_i(k_1+u_1(c_{d+1},\dots ,c_n),\dots ,k_d+u_d(c_{d+1},\dots ,c_n),c_{d+1},\dots ,c_n).$$

$$\Phi (c)={\left(\frac{\partial {\phi }_i}{\partial c_j}(c)\right)}_{d+1\le i,j\le n}$$

$$d(c,K)=\det ({\varphi }^{\prime }(c)).$$

$$Z_j^{+}=\{i\mathrm{such}\mathrm{that}{Z}_{ji}=+1\},$$

$$Z_j^{-}=\{i\mathrm{such}\mathrm{that}{Z}_{ji}=-1\},$$

$$Z_j^0=\{i\mathrm{such}\mathrm{that}{Z}_{ji}=0\}.$$

$$\frac{dA}{dt}$$

$$\frac{dB}{dt}$$

$$\frac{dC}{dt}$$

$$J=\left[\begin{array}{ccc}-k_1{C}^2-{\gamma }_a& 0& -2k_{1}AC+3k_{-1}{C}^2\\ & & \\ 0& -k_{-2}-{\gamma }_b& k_2\\ & & \\ k_1{C}^2& k_{-2}& 2k_{1}AC-3k_{-1}{C}^2-k_2-{\gamma }_c\end{array}\right]$$

$$\det (J)=$$

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_2}{dt}$$

$$\frac{d{x}_3}{dt}$$

$$\frac{d{x}_4}{dt}$$

$$\frac{d{x}_5}{dt}$$

$$x_2+x_3+x_4=e_t$$

$$M=\left[\begin{array}{ccccc}-k_1-k_3{x}_2-k_5{x}_3& -k_3{x}_1& -k_5{x}_1& k_6& 0\\ -k_3{x}_2& -k_3{x}_1& k_4& 0& 0\\ k_3{x}_2-k_5{x}_3& k_3{x}_1& -k_4-k_5{x}_1& k_6& 0\\ 0& 1& 1& 1& 0\\ 0& 0& k_4& 0& -k_7\end{array}\right]$$

$$\det (M)=$$

$$+k_3{x}_2{k}_4{k}_6{k}_7-k_3{x}_1{k}_4{k}_5{x}_3{k}_7.$$

$$\frac{dX}{dt}$$

Full text

On the multistationarity of chemical reaction networks

Marcelle Kaufmana111Corresponding author  and Christophe Souléb

(a Université libre de Bruxelles (ULB), Service de Chimie physique et Biologie théorique, Campus Plaine ULB, CP 231, Boulevard du Triomphe, 1050 Bruxelles, Belgique, mkaufman@ulb.ac.be

b CNRS et IHES, Le Bois-Marie, 35 route de Chartres, 91440 Bures-sur-Yvette, France, soule@ihes.fr)

We present a new conjecture about a necessary condition that a (bio)chemical network has to satisfy for it to exhibit multistationarity. According to a Theorem of Feliu and Wiuf, the conjecture is known for strictly monotonic kinetics. We give several examples illustrating our conjecture.

Keywords: bistability, feedback circuit, reaction network, influence graph, Jacobian matrix.

1 Introduction

This paper is a contribution to the theory of multistationarity of (bio)chemical networks. The mathematical theory of chemical reaction networks is already well developed, thanks to the pioneering work of Feinberg [9, 10], and major contributions by Craciun [5, 6, 2], Banaji [1], Mincheva-Roussel [18], Wiuf and Feliu [27], and others.

A chemical reaction network is a finite set of reactions between complexes made of chemical species. These reactions depend upon a stoichiometric matrix, and to each of them is attached a rate function which controls the kinetics of the concentration value of the different species. For instance, in the classical case of mass-action kinetics, the rate functions are monomials in the concentration values.

One is interested in deciding whether such a system allows for multiple positive steady states, a phenomenon also known as multistationarity. More precisely, from the stoichiometry and the rate functions one can define a species formation rate function $F$ having $n$ variables and $n$ values, where $n$ is the number of species. The dynamic of the system consists of the finite set of differential equations

[TABLE]

where $c(t)$ is a vector of concentrations varying with time. We shall consider a theorem of Wiuf and Feliu on the injectivity properties of (A) when the rate functions are strictly monotonic ([27], Th. 10.2), a condition fulfilled not only by mass-action kinetics, but also by Michaelis-Menten and Hill kinetics, and by others.

Our work stems from results in the biological context, namely genetic networks theory. This theory models the interactions (activation and repression) between genes in a given cell: if $c(t)$ denotes the vector of concentrations of proteins in this cell at time $t$, one assumes that $c(t)$ obeys an equation like (A). Again, one problem is to give necessary conditions for (A) to allow for multistationarity. A well known conjecture of Thomas [24] (1981), states that such a necessary condition is the fact that the interaction graph of the jacobian of $F$ contains a positive circuit. Mathematical proofs of this fact were found by Plahte et al. [20], Snoussi [21], Gouzé [13], Cinquin and Demongeot [4], and Soulé [23].

In [16], another result on gene networks was proved. It says that multistationarity of (A) requires the existence of nuclei of given signs in the interaction graph of $F$, where a nucleus is a disjoint union of circuits which contains all vertices.

An important step occurred with the paper [22] by Soliman, which proves a result “à la Thomas” for chemical reaction networks, and uses in a non-trivial way the stoichiometric matrix. In fact, as Soliman explained to us, one can recover the system (A) for any function $F$ as the dynamic of a special chemical reaction network (with ad hoc stoichiometry, see 2.2 below), so that the original Thomas’ rule becomes a special case of [22].

In this paper we are looking for a common generalization of [16] and [22]. We state this as a conjecture about chemical reaction networks (Conjecture 1). Feliu and Wiuf proved this conjecture when the rate functions are strictly monotonic (Theorem 3).

Note that an important invariant of a system like (A) is the determinant of the Jacobian matrix of $F$. But, if the system (A) obeys conservation laws, this determinant is identically trivial. In [27], Wiuf and Feliu replace this invariant by the determinant of another matrix, and they give a formula for it as the sum of principal minors of the Jacobian (Section 2, Prop. 1). By solving the conservation laws we derive another formula (Prop. 2).

In Section 3, following [22], [6] and [24], we introduce the reaction labelled influence graph of a chemical reaction network, and we use it to formulate the main result of [22], with three important corollaries. Next, we state the main result of [16] (Theorem 2). In Section 4 we present our conjecture (Conjecture 1) and, after defining strictly monotonic kinetics, we state the theorem of Feliu and Wiuf (Theorem 3).

Section 5 is devoted to concrete examples which illustrate several aspects of the theory. In particular, we use Theorem 3 to modify a reaction network in order to produce bistability. In the last example we consider a case of non-monotone kinetics.

After that, Section 6 (resp. Section 7) gives a proof of Proposition 2 (resp. Theorem 3). Finally, in Section 8, we present some conclusions and remarks.

2 Conservation laws

2.1 Definitions

Let ${\mathbb{R}}_{+}$ be the set of positive real numbers, and ${\mathbb{N}}$ the set of natural integers. Let ${\mathcal{N}}=({\mathcal{S}},{\mathcal{C}},{\mathcal{R}})$ be a chemical network. As explained in [27], Def. 3.2, ${\mathcal{N}}$ consists in three finite sets:

1. $\bullet$

A finite set of species ${\mathcal{S}}=\{S_{1},\ldots,S_{n}\}$; 2. $\bullet$

A set of complexes ${\mathcal{C}}\subset{\mathbb{N}}^{n}$; 3. $\bullet$

A set of reactions ${\mathcal{R}}=\{r_{1},\ldots,r_{m}\}\subset{\mathcal{C}}\times{\mathcal{C}}$.

It is assumed that if $y\in{\mathcal{C}}$, then $(y,y)\notin{\mathcal{R}}$ and there exists $y^{\prime}\in{\mathcal{C}}$ such that $(y^{\prime},y)\in{\mathcal{R}}$ or $(y,y^{\prime})\in{\mathcal{R}}$. We write

[TABLE]

to mean that $r=(y,y^{\prime})\in{\mathcal{R}}$.

We consider differential kinetics $K=(K_{j})$ for the network ${\mathcal{N}}$ i.e. (see [27] Def. 3.3) to each reaction $r_{j}:y_{j}\to y^{\prime}_{j}$ in ${\mathcal{R}}$ we attach a rate function

[TABLE]

which is differentiable. Let $A$ be the stoichiometric matrix, i.e. the $n\times m$ matrix with $j$-th column $y^{\prime}_{j}-y_{j}$, where $(y_{j},y^{\prime}_{j})$ lies in ${\mathcal{R}}$. The species formation rate function is the map

[TABLE]

defined by the formula

[TABLE]

If $c(t)\in{\mathbb{R}}_{+}^{n}$ is a differentiable function of $t\in{\mathbb{R}}_{+}$, and $\dot{c}(t)$ its derivative, the dynamic of ${\mathcal{N}}$ consists of the finite set of $n$ ordinary differential equations

[TABLE]

2.2 A special case

Let

[TABLE]

be any differential map. One can define as follows a chemical network such that $F$ is the corresponding species formation rate function. One takes $m=n$, and one defines $y_{j}$ (resp. $y^{\prime}_{j}$) to be zero (resp. the complex $y^{\prime}_{ij}$ such that $y^{\prime}_{jj}=1$, and $y^{\prime}_{ij}=0$ when $i\neq j$). The rate function $K_{j}$ is defined to be the $j$-th component $F_{j}$ of the function $F$. The equation (1) is then clearly satisfied.

2.3 A formula of Wiuf and Feliu

We equip ${\mathbb{R}}^{n}$ with the standard scalar product. We denote by ${\rm Im}(A)$ the image of $A$ and by ${\rm Im}(A)^{\perp}$ the orthogonal complement of ${\rm Im}(A)$. Assume ${\rm Im}(A)^{\perp}$ has dimension $d$. Consider a reduced basis of ${\rm Im}(A)^{\perp}$ i.e. a basis $(\omega^{1},\ldots,\omega^{d})$ such that $\omega_{i}^{i}=1$ if $1\leq i\leq d$, and $\omega_{j}^{i}=0$ when $j=1,\ldots,\widehat{i},\ldots,d$ (such a basis exists after reordering the species, cf. [27], Def. 5.1). The formula

[TABLE]

defines a new map

[TABLE]

When $c\in{\mathbb{R}}_{+}^{n}$ we let $J_{c}(F)$ (resp. $J_{c}(\widetilde{F})$) be the Jacobian matrix of $F$ (resp. $\widetilde{F}$) at $c$.

If $I\subset\{1,\ldots,n\}$ we let $J_{c}(F)_{I,I}$ be the submatrix of $J_{c}(F)$ with entries having indices in $I\times I$. The following is proved in [27], Prop. 5.3:

Proposition 1:

[TABLE]

We simplify the notation by defining

[TABLE]

2.4 Another formula

Let $k_{1},\ldots,k_{d}$ be $d$ real positive constants. The conservation laws

[TABLE]

reduce by $d$ the number of variables. Namely, by solving (3), if $1\leq i\leq d$, we get

[TABLE]

where $u_{i}(c_{d+1},\ldots,c_{n})$ is a linear combination of the variables $c_{d+1},\ldots,c_{n}$. When $d+1\leq i\leq n$ we define

[TABLE]

to be the function obtained from $F_{i}(c)$ by replacing $c_{1},\ldots,c_{d}$ by their expression. If $c=(c_{d+1},\ldots,c_{n})$ we let

[TABLE]

be the corresponding Jacobian matrix, and $\Phi^{\prime}(c)$ the matrix obtained from $\Phi(c)$ by replacing $k_{i}$ by $\omega^{i}\cdot c$, $i=1,\ldots,d$.

Proposition 2:

[TABLE]

3 Previous results

3.1 The reaction labelled influence graph

We keep the definitions of 2.1 above and we fix $c\in{\mathbb{R}}_{+}^{n}$. The reaction labelled infuence graph (RLIG) $G(c)=G(c,K)$ of $K$ at $c$ is defined as follows (cf. [22]). The set of vertices of $G(c)$ is the set ${\mathcal{S}}$ of species of ${\mathcal{N}}$. There is a positive (resp. a negative) edge $i\to i^{\prime}$ in $G(c)$, labelled by the reaction $r_{j}:y_{j}\to y^{\prime}_{j}$, when $\frac{\partial F_{i^{\prime}}}{\partial c_{i}}(c)$ is positive (resp. negative) and $i$ (resp. $i^{\prime}$) is a component of $y^{\prime}_{j}$ (resp. $y_{j}$). In the special case considered in 2.2, the RLIG coincides with the interaction graph attached to $F$.

A circuit in $G(c)$ is a sequence of distinct vertices $i_{1},\ldots,i_{k}$ such that there is an edge from $i_{\alpha}$ to $i_{\alpha+1}$ if $1\leq\alpha\leq k-1$, and from $i_{k}$ to $i_{1}$. The sign of $C$ is the product of the signs of its edges. Several circuits are said disjoint when they do not share any vertex. A hooping is a union of one or several disjoint circuits. If $s\geq 1$ is an integer, an $s$-hooping is a hooping containing $s$ vertices. When $s=n$, an $s$-hooping is called a nucleus[16]. The sign of a hooping $H$ is $\varepsilon(H)=(-1)^{p+1}$, where $p$ is the number of positive circuits contained in $H$ (see [8]).

Given a hooping $H$, the restriction of ${\mathcal{N}}$ to $H$ is the sub-network ${\mathcal{N}}_{H}$ of ${\mathcal{N}}$ where reactions $r_{j}$ not appearing in $H$ are omitted. Note that, given a vertex $i\in{\mathcal{S}}$, there is at most one edge issued from $i$ in the RLIG of ${\mathcal{N}}_{H}$. Therefore the stoichiometric matrix of ${\mathcal{N}}_{H}$ is a square matrix. Let $\Lambda_{H}$ be its determinant. We introduce the following

Definition 1: A hooping $H$ of $G(c)$ is called admissible when the stoichiometric matrix of ${\mathcal{N}}_{H}$ is invertible (i.e. $\Lambda_{H}\neq 0$) .

3.2 A result of Soliman

A zero $c\in{\mathbb{R}}_{+}^{n}$ of the function $F$ is called non-degenerate when the Jacobian matrix $\left(\frac{\partial F_{i^{\prime}}}{\partial c_{i}}(c)\right)_{1\leq i,i^{\prime}\leq n}$ is invertible.

Theorem 1 [22]: Assume $F$ has several non-degenerate zeroes. Then there exists $c\in{\mathbb{R}}_{+}^{n}$ such that $G(c)$ has a positive circuit contained in an admissible hooping.

3.3

Soliman derived three corollaries from Theorem 1:

Corollary 1 [22]: A necessary condition for the multistationarity of a biochemical system is that there exists a positive circuit in its RLIG using at most once each reaction.

Indeed, if a circuit uses more than once some reaction, it is not admissible.

Corollary 2 [22]: A necessary condition for the multistationarity of a biochemical system is that there exists a positive circuit in its RLIG not using both forward and backward directions of any reversible reaction.

Indeed, if a circuit uses both forward and backward directions of some reversible reaction, it is not admissible.

And, by a similar argument:

Corollary 3 [22]: A necessary condition for the multistationarity of a biochemical system is that there exists a positive circuit in its RLIG not using all species involved in a conservation law.

3.4 A result of Kaufman, Soulé and Thomas

In this section we assume that we are in the special case considered in 2.2. Consider the finite set ${\mathcal{F}}$ of real functions on ${\mathbb{R}}_{+}^{n}$ of the form $\pm\underset{i\in I}{\prod}(\partial F_{i}/\partial\,x_{\tau(i)})$ where $I\subset\{1,\ldots,n\}$ is any subset and $\tau:I\to\{i,\ldots,n\}$ is an injective map. We consider the following condition:

(C)  Given two functions $f$ and $g$ in ${\mathcal{F}}$ such that $f$ is not identically zero and $g$ is strictly positive somewhere in ${\mathbb{R}}_{+}^{n}$, there exists $x\in{\mathbb{R}}_{+}^{n}$ such that $f(x)\neq 0$ and $g(x)>0$.

Note that (C) is very often fulfilled.

We say that the interaction graph $G$ has a variable nucleus when there exist two points $c$ and $d$ in ${\mathbb{R}}_{+}^{n}$ and a nucleus of opposite signs in $G(c)$ and $G(d)$ (see [16] 3.1 for a more precise definition).

Theorem 2 [16]: Let $F$ be such that (C) is satisfied. If $F$ has several non-degenerate zeroes, then:

1. i)

either there exists $c\in{\mathbb{R}}_{+}^{n}$ such that $G(c)$ has two nuclei of opposite signs; 2. ii)

or $G$ has a variable nucleus.

This theorem is presented as a conjecture in [25].

4 The main result

4.1 A conjecture

Assume $m\geq 1$ is any positive integer. In view of Theorems 1 and 2 it seems reasonable to make the following conjecture:

Conjecture 1: Assume that $F$ has several non-degenerate zeroes and that (C) is satisfied. If $F$ has $d$ conservation laws, and $s=n-d$,

1. i)

either there exists $c\in{\mathbb{R}}_{+}^{n}$ such that $G(c)$ has two admissible $s$-hoopings of opposite signs; 2. ii)

or $G$ has a variable admissible $s$-hooping.

4.2 Monotonic kinetics

The Conjecture 1 is known when the kinetics are strictly monotone. To describe this result, we introduce more definitions from [27].

Let $Z=(Z_{ji})$ be a matrix of type $m\times n$ such that each entry $Z_{ji}$ is equal to $+1,-1$ or [math]. We let

[TABLE]

[TABLE]

and

[TABLE]

We assume that, for every $c\in{\mathbb{R}}_{+}^{n}$ and every $j\in\{1,\ldots,m\}$, $K_{j}(c)>0$. We say that $K$ is strictly monotonic with respect to $Z$ ([27], Def. 9.3) if, for all $j\in\{1,\ldots,m\}$, the function $K_{j}$ is

1. i)

strictly increasing in $c_{i}$ when $i\in Z_{j}^{+}$, 2. ii)

strictly decreasing in $c_{i}$ when $i\in Z_{j}^{-}$, 3. iii)

constant in $c_{i}$ when $i\in Z_{j}^{0}$.

We denote by $K(Z)$ the set of kinetics $K$ which are strictly monotonic with respect to $Z$.

Let $A$ be the stoichiometric matrix. The matrix $A$ is called not $Z$-injective over $K(Z)$ when there exist $K\in K(Z)$, $a\in{\mathbb{R}}_{+}^{n}$ and $b\in{\mathbb{R}}_{+}^{n}$ such that $F(a)=F(b)$.

The network ${\mathcal{N}}$ is called non-degenerate when there exists $c\in{\mathbb{R}}_{+}^{n}$ and $K\in K(Z)$ such that $d(c,K)\neq 0$.

The following result relates the fact that $F$ is not injective to the existence of two admissible $s$-hoopings of opposite signs in the RLIG. It proves Conjecture 1 when $K$ is strictly monotone: there are no variable $s$-hoopings and, if F has several non degenerates zeroes, A is not Z-injective over K(Z) (we do not need to assume that hypothesis (C) is satisfied).

Theorem 3 [27, 12]

Assume that ${\mathcal{N}}$ is non-degenerate and that $A$ is not $Z$-injective over $K(Z)$. Let $s=n-d$ be the rank of $A$. Then, for every $c\in{\mathbb{R}}_{+}^{n}$ and $K\in K(Z)$, the RLIG $G(c)$ contains two admissible $s$-hoopings of opposite signs.

5 Examples

To illustrate Conjecture 1 we consider four simple examples. For the first three examples we simply use mass-action kinetics to describe the reaction rates. The last example involves a biochemical reaction with non-monotone kinetics.

Example 1. An autocatlytic reaction

Our first example consists of the simple autocatalytic reaction system shown in Fig. 1. It is described by the differential equations

[TABLE]

where we consider a constant source $S_{A}$ of $A$ and $S_{C}$ of $C$. The associated Jacobian matrix is

[TABLE]

which contains the contributions, if any, of each reaction for each of the variables. The corresponding reaction labelled influence graph (RLIG) is depicted in Fig. 2.

In this RLIG, there are two positive circuits using both the forward and backward directions of a reversible reaction. According to Definition 1 and corollary 2 of Theorem 1 they cannot form admissible hoopings because their corresponding stoichiometric matrices are singular. There remain two admissible positive nuclei formed by the positive loop $R_{1}$ on $C$, the negative loop $R_{a}$ on $A$ and, respectively, the negative loop $R_{b}$ or $R_{-2}$ on $B$. The RLIG also comprises several admissible negative nuclei composed of negative self-loops. The necessary conditions of Theorem 1 and Theorem 3 are satisfied and for suitable parameter values this system exhibits bistability (Fig. 3).

According to (20) below, the determinant expansion of the Jacobian matrix (5) is a sum of terms indexed over all the admissible nuclei

[TABLE]

where the terms associated with the admissible positive and negative nuclei are of opposite sign.

If $\gamma_{a}=0$, Theorem 1 does not rule out the possibility of several steady states as there still exists a positive circuit $R_{1}$ included in two admissible 2-hoopings $R_{1}R_{b}$ and $R_{1}R_{-2}$. However, in this case, Theorem 3 is no longer satisfied. This can be seen on the RLIG in Fig. 2 and on the determinant expansion (6) which, for $\gamma_{a}=0$, contain only negative admissible nuclei composed of negative self-loops. The system cannot have multiple steady states and its unique steady state can be determined analytically.

Example 2. Conservation laws

To illustrate Theorem 3 in the presence of a conservation law, we consider a model of reversible substrate inhibition studied by Mincheva and Roussel [18]. This biochemical reaction system and corresponding RLIG are shown in Figs. 41 and 42.

As can be seen on Fig. 42, the 3-circuit between $S_{1}$, $ES_{1}$ and $E$, involving the reactions $R_{5}$, $R_{4}$ and $R_{3}$, is the only positive circuit that satisfies the necessary conditions of the three corollaries of Theorem 1. It forms an admissible positive 4-hooping with the negative loop $R_{7}$ on $P$.

The system of differential equations is

[TABLE]

where $x_{i}$, $i=1,\ldots 5$, are the concentrations of species $S_{1}$, $E$, $ES_{1}$, $S_{1}$ $ES_{1}$ and $P$, respectively. The enzymatic species are linked by the conservation relation

[TABLE]

where $e_{t}$ is the total enzyme concentration.

Due to the conservation relation, the Jacobian matrix of system (7) is singular. Therefore, following the theory developed in the sections 2 and 6 (see also [11], 3.1), we construct a $5\times 5$ matrix $M$ by replacing the fourth row of the Jacobian matrix of eqs (7) by the vector $w^{1}=(0,1,1,1,0)$ defined by the conservation law (3) with $c=x$ and $k_{1}=e_{t}$

[TABLE]

and

[TABLE]

The determinant of $M$ contains a sum of terms associated with the positive admissible $s$-hooping (with $s=n-d=4$) and with the four admissible negative $s$-hoopings which are all formed by a combination of self-loops. The necessary condition for multistationarity is thus satisfied and for suitable parameter values this system is bistable [18].

Example 3. Synthetic approach

Let us consider the well-known Brusselator [17, 19, 26] shown in Fig. 5(a). The initial and final product concentrations $A,B,D,E$ are maintained time-independent.

This model system is described by the differential equations

[TABLE]

with the associated Jacobian matrix

[TABLE]

Equations (11) admit a unique steady state solution

[TABLE]

that becomes unstable through a Hopf bifurcation leading to sustained oscillations (Fig. 7A).

From the corresponding RLIG (Fig. 5(b)) or Jacobian matrix (12) it can be seen that the Brusselator, in its original form, does not satisfy the conditions of Theorem 3 and Theorem 1: the positive loop $R_{2}$ of $X$ on itself and the positive circuit between $X$ and $Y$, involving reactions $R_{1}$ and $R_{2}$, are not admissible. The determinant of Jacobian (12) contains only one, negative nucleus

[TABLE]

However, inspection of the RLIG in Fig. 5(b) suggests that a slight modification of this influence graph consisting in the addition of a linear decay reaction for species $Y$ (Figs. 6(a) and (b)) allows to satisfy the conditions of both theorems. In addition to three terms corresponding to negative nuclei, the determinant expansion of the modified Jacobian now contains a term associated with a positive nucleus formed by the positive loop $R_{2}$ on $X$ and the negative loop $R_{4}$ on $Y$

[TABLE]

where $k_{4}$ is the kinetic constant of reaction $R_{4}$.

The steady state solutions are now given by a cubic equation and as shown in Fig. 7B, for suitable parameter values, this slightly modified Brusselator becomes a bistable system without losing the capacity to exhibit oscillations (Fig. 7D). Excitability linked to multistationarity is also observed (Fig. 7C). In Fig. 7D one can see that $k_{4}$ favors bistability while $k_{3}$ favors sustained oscillations.

Example 4. A variable nucleus

To illustrate point ii) of conjecture 1, we consider a substrate cycle model derived from the system studied in detail by Hervagault and Canu [14, 3]. In this system two metabolites $S_{1}$ and $S_{2}$ are interconverted by two antagonist enzymes. Metabolite $S_{1}$ is converted into metabolite $S_{2}$ by enzyme $E_{1}$ and $S_{2}$ is converted in turn into $S_{1}$ by enzyme $E_{2}$. In addition, enzyme $E_{1}$ is inhibited by excess of its substrate $S_{1}$. This cyclic pathway is depicted in Fig. 8 together with the corresponding RLIG. It is governed by the differential equations

[TABLE]

where we consider a constant source $v_{s}$ of $S_{1}$ and linear decays of $S_{1}$ and $S_{2}$. In this example, the reaction rate of the first reaction is a non-monotone function of its argument (see insert in Fig. 9) and does not satisfy the premises of Theorem 3. Note that this function is a variant of Example 2.

The RLIG contains a sign-variable self-loop $R_{1}$ on $S_{1}$ which together with the negative sef-loop $R_{4}$ on $S_{2}$, forms a variable nucleus.

From the Jacobian matrix of (16)

[TABLE]

one observes that the variable self-loop has a positive effect on $S_{1}$ when $S_{1}>\sqrt{(K_{1}/k)}$ and it provides the possibility of a change of sign of the determinant expansion of $J$

[TABLE]

As shown in Fig. 9, for appropriate values of the parameters, this model system has several steady states.

6 Proof of Proposition 2

To prove Proposition 2 we first notice that, since $(\omega^{1},\ldots,\omega^{d})$ is reduced, when $k=1,\ldots,d$, we have

[TABLE]

Therefore, if $j=d+1,\ldots,n$ and $i=d+1,\ldots,n$, we have

[TABLE]

Consider the $n$ by $n$ matrix $J_{c}(\tilde{F})$. For every $k=1,\ldots,d$, substract from the $j$-th column of this matrix, $j=d+1,\ldots,n$, the product of the $k$-th column by $\omega_{j}^{k}$. The determinant does not change by this operation, and the matrix $M$ obtained this way, is such that the first $d$ lines of $M$ are equal to $\delta_{ij}$, $i=1,\ldots,d$, $j=1,\ldots,n$.

From (19) we deduce that

[TABLE]

Since $\det(M)=\det(J_{c}(\widetilde{F}))$ Proposition 2 follows.

7 Proof of Theorem 3

7.1

Theorem 3 is due to Feliu and Wiuf, see [12] (based on [27], Section 11). The proof can be described as follows. According to [27], Theorem 10.2, and the line after it, since $A$ is not injective over $K$, property i) in [27], loc. cit., does not hold, i.e. there exist $c_{1}\in{\mathbb{R}}_{+}^{n}$ and $K_{1}\in K(Z)$ such that $d(c_{1},K^{1})=0$.

On the other hand, since ${\mathcal{N}}$ is non degenerate, there exists $c_{2}\in{\mathbb{R}}_{+}^{n}$ and $K^{2}\in K_{d}(Z)$ such that $d(c_{2},K^{2})\neq 0$, say $d(c_{2},K^{2})>0$. By (2) this implies that there exists $I\subset\{1,\ldots,n\}$ with $\#\,I=s$ such that $\det(J_{c_{2}}(F^{2})_{I,I})>0$, where $F^{2}$ is the species formation rate function defined from $K^{2}$. Let $H_{s}(c_{2},K^{2})$ be the set of s-hoopings in $G(c_{2},K^{2})$ and ${\mathcal{H}}_{2}\subset{\mathcal{H}}_{s}(c_{2},K^{2})$ (resp. ${\mathcal{H}}_{1}\subset{\mathcal{H}}_{s}(c_{1},K^{1})$) be the collection of hoopings in $G(c_{2},K^{2})$ (resp. $G(c_{1},K^{1})$) having $I$ as set of vertices. We also write ${\mathcal{H}}={\mathcal{H}}_{s}(c,K)$.

7.2

We shall need some notions and results from [15] and [22]. If $s_{i}\in{\mathcal{S}}$, $r_{j}\in{\mathcal{R}}$, $c\in{\mathbb{R}}_{+}^{n}$ and $K\in K_{d}(Z)$ we let

[TABLE]

and

[TABLE]

For a path $P$ in $G(c,K)$ we let $\lambda_{SR}(P)$ (resp. $\lambda_{RS}(P)$) be the product of all the numbers $\lambda(s_{i},r_{j})$ (resp. $\lambda(r_{j},s_{i})$) where $S_{i}\xrightarrow{\ R_{j}\ }$ (resp. $\xrightarrow{\ R_{j}\ }S_{i}$) is an edge in $P$. We define also

[TABLE]

If $H,H^{\prime}\in{\mathcal{H}}={\mathcal{H}}_{s}(c,K)$ we write that $H\sim H^{\prime}$ when $H$ and $H^{\prime}$ have the same species to reactions arcs $S_{i}\xrightarrow{\ R_{j}\ }$. Note that $\lambda_{SR}(H)$ and $\Lambda(H)$ depend only on the equivalence class $[H]\in{\mathcal{H}}/\sim$. According to [15] Thm. 6.7 (see also [22], Thm. 2), we have

[TABLE]

and

[TABLE]

where $\sigma(H^{\prime})=\pm\,1$. As a consequence ([15] Lemma 6.2 and [22], (1)) we get

[TABLE]

7.3

Let us come back to the situation of §7.1. Let $J=J_{c_{2}}(F^{2})_{I,I}$. Since $\det(J)>0$, it follows from (20), (21) and (22) that there exists $H_{2}\in{\mathcal{H}}_{2}$ such that $\Lambda(H_{2})\neq 0$ and $\sigma(H_{2})\,\lambda(H_{2})>0$.

We claim that there exists also $H_{3}\in{\mathcal{H}}_{2}$ such that $\Lambda(H_{3})\neq 0$ and $\sigma(H_{3})$ $\lambda(H_{3})<0$. Suppose, by contradiction, that, for any $H\in{\mathcal{H}}_{2}$ such that $\Lambda(H)\neq 0$, we have $\sigma(H)\,\lambda(H)\geq 0$. Since the kinetics are strictly monotonic, for any $H\in{\mathcal{H}}_{s}(c,K)$ the sign of $\sigma(H)\,\lambda(H)$ does not depend on $c$ and $K$, but only on the arcs contained in $H$. Therefore, for any $s$-hooping $H^{\prime}$ in ${\mathcal{H}}_{1}$ such that $\Lambda(H^{\prime})\neq 0$, one has $\sigma(H^{\prime})\,\lambda(H^{\prime})\geq 0$, and there exists $H_{1}\in{\mathcal{H}}_{1}$ such that $\sigma(H_{1})\,\lambda(H_{1})>0$. By (22) this implies that $\det(J_{c_{1}}(F^{1})_{I,I})>0$. But $d(c_{1},K_{1})=0$, so we conclude by (2) that there exists $I^{\prime}$ of cardinality $s$ such that $\det(J_{c_{1}}(F^{1})_{I^{\prime},I^{\prime}})<0$. Therefore, when $H\in{\mathcal{H}}$ and $\Lambda(H)\neq 0$, the number $\sigma(H)\,\lambda(H)$ takes both positive and negative values.

To conclude the proof of Theorem 3, it remains to note that the sign $\varepsilon(H)$ is $(-1)^{s+1}$ times the sign of $\sigma(H)\,\lambda(H)$ (cf. [8], Appendix, Lemma 2). Therefore $G(c)$ contains two admissible $s$-hoopings of opposite signs.

8 Conclusion

In this paper, we conjecture that the multistationarity of a (bio)chemical network requires the existence, in the reaction labelled influence graph, of two admissible hoopings of maximal rank and opposite signs or a variable admissible hooping of maximal rank. This conjecture has been proved [12] under the hypothesis that the kinetic is strictly monotonic in the sense of [27] . This includes mass-action kinetics, Michaelis-Menten and Hill kinetics, and many more. The proof that we present makes an essential use of the connection between the reaction labelled influence graph and the determinant of the Jacobian (or its substitute when there exists conservation laws). We hope that this theorem can be extended to include arbitrary differentiable metrics (as suggested by Example 4).

This result can be used to rule out the possibility of multistationarity in situations where other criteria do not (Example 1). Given a network which does not fulfill our condition for multistationarity, following Theorem 3, inspection of the reaction labelled influence graph or Jacobian matrix also suggests ways of modifying the network so that it fulfills our condition (Example 3).

This theorem is thus an illustration of R.Thomas’ insight that, in biology, one should often look for necessary rather than sufficient conditions.

Acknowledgements

We dedicate this work to the memory of our friend René Thomas. C.S. expresses his admiration for the scientific trajectory of R.Thomas, from fundamental biology to pure mathematics.

We are grateful to E. Feliu, D. Gonze and S. Soliman for their help. The bifurcation diagrams were generated with AUTO [7]. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.