 Research Article
 Open Access
 Published:
Finding the positive feedback loops underlying multistationarity
BMC Systems Biologyvolume 9, Article number: 22 (2015)
Abstract
Background
Bistability is ubiquitous in biological systems. For example, bistability is found in many reaction networks that involve the control and execution of important biological functions, such as signaling processes. Positive feedback loops, composed of species and reactions, are necessary for bistability, and generally for multistationarity, to occur. These loops are therefore often used to illustrate and pinpoint the parts of a multistationary network that are relevant (‘responsible’) for the observed multistationarity. However positive feedback loops are generally abundant in reaction networks but not all of them are important for understanding the network’s dynamics.
Results
We present an automated procedure to determine the relevant positive feedback loops of a multistationary reaction network. The procedure only reports the loops that are relevant for multistationarity (that is, when broken multistationarity disappears) and not all positive feedback loops of the network. We show that the relevant positive feedback loops must be understood in the context of the network (one loop might be relevant for one network, but cannot create multistationarity in another). Finally, we demonstrate the procedure by applying it to several examples of signaling processes, including a ubiquitination and an apoptosis network, and to models extracted from the Biomodels database. The procedure is implemented in Maple.
Conclusions
We have developed and implemented an automated procedure to find relevant positive feedback loops in reaction networks. The results of the procedure are useful for interpretation and summary of the network’s dynamics.
Background
Bistability, and multistationarity in general, is ubiquitous in biological systems ranging from biochemical networks to epidemiological and ecosystems [14]. It is considered an important biological mechanism for controlling cellular and bacterial behaviour and developmental processes in organisms, and it is closely linked to the idea of the cell as a decision making unit, where a continuous input is converted to an on/off response corresponding to two distinct states of the cell [5,6].
The question of bistability therefore arises naturally in many contexts. Many studies aim to demonstrate that in a given biochemical system, bistability can or cannot occur [2,3,710]. There are several methods that can be used to address whether a network is multistationarity or not, see for example [1119] and the references therein. Some of these methods are implemented in the CRNT toolbox [20] or in CoNtRol [21]. More general there has been some interest in formal methods that connect the network structure to the dynamic behaviour of the system, see e.g. [15,2227].
One qualitative network feature has in particular been linked to multistationarity, namely the existence of a positive feedback loop. A positive feedback loop consists of a sequence of species such that each species affects the production of another species, either positively or negatively, and such that the number of negative influences is even. The idea of associating positive feedback loops with bistability goes back to Jacob and Monod who introduced it in the context of gene regulatory networks [28]. It was later formalised by Thomas in the form of a conjecture [29], which was finally proved by Soulé [30], see also [31,32].
Soulé considers dynamical systems of the form
where x=x(t), x=(x _{1},…,x _{ n }) is the vector of species concentrations, $\dot {x}=dx/dt$ is the derivative of x with respect to time t, and f is the socalled speciesformation rate function, which specifies the instantaneous change in the concentrations.
The work of Soulé is based on the socalled interaction graph [30]. This graph encodes how the variation of one species concentration depends on the concentration of the other species. It is built from the Jacobian matrix J _{ f }(x ^{∗}) of f evaluated at a point x ^{∗}, such that the nonzero entries of J _{ f }(x ^{∗}) correspond to directed edges of the graph and the signs of the entries are edge labels. Soulé proved that the existence of a positive feedback loop in the interaction graph is a necessary condition for f(x) to have multiple zeros. In other words, it is a necessary condition for multistationarity to exist in the ODE system (1).
Soulé’s approach is often not useful for many reaction networks, such as enzymatic signaling networks, because the edge labels are not constant, but depend on the concentrations of the species, that is, for some concentrations a label might be positive, for others it might be negative. A refinement of Soulé’s work is based on the socalled directed speciesreaction graph (DSRgraph) [13,3335]. If f in (1) is obtained from a reaction network, then it decomposes in the form
where A is the stoichiometric matrix of the network and v(x) the vector of reaction rates. The DSRgraph uses this particular structure.
The DSRgraph is a bipartite graph with nodes labeled by the species and the reactions of the reaction network. Labeled directed edges from species nodes to reaction nodes and from reaction nodes to species nodes are derived from the vector of reaction rates v(x) and the stoichiometric matrix A, respectively. Compared to the interaction graph, the DSRgraph makes use of the explicit decomposition (2) of f.
It has been shown that the existence of positive feedback loops in the DSRgraph is a necessary condition for the system (2) to admit multistationarity [34].
Based on these results it has become standard to highlight positive feedback loops in multistationary reaction networks, eg. [1,2]. The loops are typically found using intuitive reasoning that might overlook the existence of other relevant positive feedback loops or might select positive feedback loops that are not related to the existence of multistationarity. Here we provide a method, based on theoretical considerations, to classify all positive feedback loops of a multistationary network into those that are related to the observed multistationarity and those that are not. In other words, we determine the positive feedback loops that when they all are broken, multistationarity disappears.
The question needs to be understood in the context of the whole network and not in isolation: a particular positive feedback loop that is responsible for multistationarity in one network might appear in another network that cannot have multiple steady states.
We present an automated procedure to determine the positive feedback loops that contribute to multistationarity. The procedure is based on various ideas from previous work by us and others. In particular, it builds on the injectivity property applied to an ODE system of the form (2), as described in [13]. The procedure can be applied to any network, provided that the components of the vector v(x) are strictly monotone with respect to all variables. This is a mild assumption that is fulfilled for typical kinetics such as massaction kinetics and MichaelisMenten kinetics.
In Methods, we introduce the necessary background material. As part of this we explain why positive feedback loops are necessary for multistationarity and how this relates to the DSRgraph. In Results, we present the automated procedure and how it selects the relevant positive feedback loops out of all positive feedback loops. We further apply the procedure to examples of multistationary reaction networks involved in cell signaling. We also consider the networks in the Biomodels database [36] and apply the procedure to all noninjective networks (injective networks cannot be multistationary, see below). This provides an overview of the landscape of relevant positive feedback loops occurring in documented reaction networks.
Methods
In this section we introduce the different ideas we need to construct the automated procedure. We use the formalism of Chemical Reaction Network Theory (CRNT) [37]. An ODE system is built from a set of reactions and reaction rates, which we explain in the section below.
Reaction networks
A reaction network, or simply a network, consists of a set of species {X _{1},…,X _{ n }} and a set of reactions of the form:
where α _{ ij },β _{ ij } are nonnegative integers, called the stoichiometric coefficients of the reactants and the products, respectively. As a running example we use the network in Figure 1. It has three species, X _{cyt}, X _{nuc}, X $^{*}_{\text {nuc}}$ , which are different forms of the Cdk1cyclin B1 complex, and four reactions [1].
We denote the concentration of the species X _{1},…,X _{ n } by lowercase letters x _{1},…,x _{ n }. The evolution of the species concentrations with respect to time is modelled as an ODE system in the following way. We let A=(a _{ ij }) be the stoichiometric matrix of the network:
that is, the (i,j)th entry encodes the net stoichiometric coefficient of species X _{ i } in reaction r _{ j }. The vector (a _{1j },…,a _{ nj }) is called the reaction vector of reaction r _{ j }.
The rate of reaction r _{ j } is a function $v_{j}\colon \Omega _{v} \rightarrow \mathbb {R}_{\geq 0}$ , where $\mathbb {R}_{>0}^{n}\subseteq \Omega _{v} \subseteq \mathbb {R}^{n}_{\geq 0}\) and Ω _{ v } is the set of possible species concentrations. A typical choice of v=(v _{1},…,v _{ m }) is massaction kinetics. In this case
where κ _{ j } is a positive reaction rate constant and 0^{0}=1 by convention. Putting the pieces together provides a model for the evolution of the species concentrations over time:
Returning to Figure 1, we let x _{1},x _{2},x _{3} be the concentrations of X _{cyt}, X _{nuc}, X $^{*}_{\text {nuc}}$ , respectively. Following [1], one model of the network is:
where κ _{1},…,κ _{4},K>0 are parameters. It takes the form (4) with
and $\Omega _{v}=\mathbb {R}^{n}_{\ge 0}$ . Observe that the phosphorylation reaction $\mathrm {X}_{\text {nuc}} \rightarrow \mathrm {X}^{*}_{\text {nuc}}$ has a reaction rate that depends on both the concentration of the reactant X _{nuc} and the concentration of the product X $^{*}_{\text {nuc}}$ . We also consider an alternative model in which the rate of X _{nuc} phosphorylation depends on x _{2} only:
This alternative model is also consistent with the set of reactions in Figure 1, but the third reaction is now independent of the amount of X $^{*}_{\text {nuc}}$ .
Multistationarity
In general the trajectory of the ODE system (4) determined by an initial positive condition is confined to a particular subset of $\mathbb {R}_{\geq 0}^{n}$ , called a stoichiometric compatibility class [37]. For instance, in the running example, the quantity T(x)=x _{1}+x _{2}+x _{3} is conserved through time and determined by its value at time 0. This equation (called a conservation law), and the value it takes, characterises the stoichiometric compatibility class. Two trajectories with different initial conditions but with the same value of T(x) are confined to the same stoichiometric compatibility class.
The stoichiometric compatibility classes are defined as (see [37])
where x _{0}=x(0) in $\mathbb {R}_{>0}\) is the initial condition and im(A) denotes the image of A. That is, the trajectories are restricted to the space spanned by the reaction vectors. Any trajectory that starts in the interior of $\mathcal C_{0}$ , stays there, but might be attracted towards the boundary.
A reaction network is said to be multistationary if there exist two distinct positive steady states in a stoichiometric compatibility class (but not necessarily in all classes) [37]. Equivalently, if there exist distinct positive $x,y\in \mathbb {R}^{n}_{>0}\) such that A v(x)=A v(y)=0 and x−y∈im(A). A network with one positive steady state and one steady state at the boundary is therefore not multistationary in this terminology.
The reaction network in Figure 1 is multistationary for some choice of parameters with the rate vector in (7) [1], but not with the rate vector in (8) for all choice of parameter values (which will be shown later).
Influence matrix
The concept of a positive feedback loop is associated with structural network properties and qualitative features of the reaction rates. Therefore, we assume some regularity on the reaction rates which we will encode into an abstract symbolic matrix, called the influence matrix [13] (see also [35]). A feedback loop does not depend on specific parameters or the specific functional form of the reaction rates.
To proceed, we assume that the function v _{ j }(x) is strictly monotone in each variable x _{ i } and define the influence matrix Z=(z _{ ij }) as
where γ _{ ij } are symbolic variables.
The influence matrices associated with the two reaction rate vectors in (7) and (8) are given by
and
respectively. In (10) and (11), all influences are zero or positive.
In the following sections we will develop the graphical framework that we use to find the relevant positive feedback loops. It builds on the DSRgraph [13] (see also [34,35] for alternative definitions of the DSRgraph). Subsequently, we define circuits and nuclei based on Soulé’s work [30].
DSRgraph
We define the DSRgraph as a labelled bipartite directed graph with node set {X _{1},…,X _{ n },r _{1},…,r _{ m }} and such that:

(a)
There is an edge from X _{ i } to r _{ j } with label z _{ ij } if z _{ ij }≠0.

(b)
There is an edge from r _{ j } to X _{ i } with label a _{ ij } if a _{ ij }≠0.
We refer to the signed DSRgraph as the graph identical to the DSRgraph given by (a)(b), but with the labels replaced by their signs. The (signed) DSRgraph of the running example with A as in (6) and Z as in (10) is shown in Figure 2. The (signed) DSRgraph with Z as in (11) is identical to that in Figure 2, with the edge from X $^{*}_{\text {nuc}}$ to r _{3} removed.
Circuits and nuclei
The positive feedback loops of the DSRgraph are instances of circuits, cycles of directed edges in a graph. Formally, a circuit in a graph G is a sequence of distinct nodes i _{1},…,i _{ q } such that there is a directed edge from i _{ k } to i _{ k+1} for all k≤q−1 and one from i _{ q } to i _{1}. A circuit must involve at least one edge. The label of a circuit C, denoted ℓ(C), is the product of the labels of the edges in the circuit. Two circuits are disjoint if they do not have any common nodes. Three circuits are highlighted in Figure 2(A) (two blue and one red), each involving two nodes.
A circuit with positive label is a positive feedback loop, and similarly, a circuit with negative label is a negative feedback loop. The three positive feedback loops of the running example are shaded in Figure 2B. They correspond to shuttling of the complex between the nucleus and the cytoplasm, activation and deactivation of X _{nuc}, and selfactivation of X _{nuc} (the rate of reaction r _{3} increases with x _{3}, that is, the production of X $^{*}_{\text {nuc}}$ increases with the amount of X $^{*}_{\text {nuc}}$ ).
A knucleus is a collection of disjoint circuits which involves k nodes [30]. The label ℓ(D) of a knucleus D is the product of the labels of the edges in the nucleus. Let a _{1},a _{2} be the number of circuits in the nucleus that have odd (resp. even) number of species nodes and let a=a _{1}+a _{2}. The sign of a knucleus is defined as σ(D)=(−1)^{a} _{2}. That is, if D=C _{1}∪⋯∪C _{ a } is a disjoint union of circuits, then
In the DSRgraph, any circuit involves an equal number of species and reaction nodes and, hence, an even number of edges. We will consider nuclei with k=2s edges, where s is the rank of the matrix A. The reason for considering k=2s will become clear in the Section ‘The polynomial and circuits’. Let D=C _{1}∪⋯∪C _{ a } be a 2snucleus of the DSRgraph. If none of the circuits are positive feedback loops, then the sign of σ(D)ℓ(D) is (−1)^{s}. Indeed, if all circuits have negative labels, that is, ℓ(C _{ i }) has negative sign for all i, then
Because D is a 2snucleus, it contains s species nodes. Let n _{ i } be the number of species nodes in circuit C _{ i }, such that s=n _{1}+⋯+n _{ a }. We have that n _{ i } is odd for a _{1} of the circuits and even for a _{2} of the circuits. Therefore, $(1)^{s} = (1)^{n_{1}+\dots +n_{a}} = (1)^{a_{1}},$ and
if there are no positive feedback loops in D. This result is also in [35], where it is stated using a different terminology.
We next we introduce the concept of injectivity and then link it to circuits and positive feedback loops in a separate section.
Injectivity
In this section we study the injectivity of the function x↦A v(x), $x\in \mathcal C_{0}$ (4). The injectivity property has been discussed in different contexts, see for example [12,13,15,35,38,39]. We use here the approach given in [13].
In the next section we link this injectivity property for all positive stoichiometric compatibility classes to the nonexistence of positive feedback loops in the DSRgraph. With other words, if all feedback loops are negative then the function is injective on all positive stoichiometric compatibility classes. In particular, there cannot exist two distinct points $x,y\in \mathbb {R}^{n}_{>0}\) in the same stoichiometric compatibility class such that A v(x)=A v(y)=0, that is, the network cannot be multistationary.
To decide whether the function A v(x) is injective on all positive stoichiometric compatibility classes for any v(x) with given influence matrix Z, we rely on a method previously developed by us [12,13]. We will now explain this method.
Given a stoichiometric matrix A and an influence matrix Z, we define a polynomial p _{ A,Z } of degree s=rank(A), in as many variables as there are nonzero entries of Z. For this, let M=A Z ^{t} and let {ω ^{1},…,ω ^{d}} be a basis of im(A)^{⊥}, which we assume to be Gaussreduced. Further, let i _{1},…,i _{ d } be the indices of the first nonzero entries of ω ^{1},…,ω ^{d}, respectively. We define a symbolic n×n matrix, $\widetilde {M}$ , by replacing the i _{ j }th row of M with ω ^{j} (cf. [13], Section 5). The polynomial p _{ A,Z } is defined as
which can be written as a sum of terms depending on the variables γ _{ ij }, by expanding the determinant. Each nonzero term is of the form $c\prod _{k=1}^{s} \gamma _{i_{k}j_{k}}$ where c is a coefficient (positive or negative) and all i _{ k }, respectively j _{ k }, are distinct.
It is a result of [12,13] that if p _{ A,Z } is not identically zero and all nonzero coefficients of p _{ A,Z } have the same sign, then the function A v(x) is injective on each positive stoichiometric compatibility class and, hence, the network cannot be multistationary. As a consequence, p _{ A,Z } has coefficients of opposite sign whenever the network is multistationary. If the coefficients do not have the same sign, then the network might be multistationary, but it cannot be concluded from the test.
Consider the matrix A given in (6) and Z _{1} in (10). We choose {(1,1,1)} as a basis of im(A)^{⊥} and obtain
There are both positive and negative terms, hence multistationarity cannot be excluded. For Z _{2} in (11), the polynomial $p_{A,Z_{2}}$ is obtained from (14) by setting γ _{3,3}=0,
In this case all terms have the same sign and thus, the network cannot be multistationary. This holds for any choice of rate functions with influence matrix Z _{2}.
The polynomial and circuits
Finally, we link injectivity and the polynomial p _{ A,Z } to positive feedback loops. It is shown in [13] that each term of the polynomial p _{ A,Z } can be identified with a collection of disjoint unions of circuits in the DSRgraph G. Specifically, given subsets I,J⊆{1,…,n} of cardinality s, let D _{ s }(I,J) be the set of 2snuclei of G with node set {X _{ i } i∈I}∪{r _{ j } j∈J}. Then
where the sets I,J in the sum have cardinality s (cf. [13], Section 11).
Let D∈D _{ s }(I,J) be a 2snucleus. If none of the circuits of the nucleus are positive feedback loops then the sign of σ(D)ℓ(D) is (−1)^{s}, see Section ‘Circuits and nuclei’ and (13). Hence if the DSRgraph contains a positive feedback loop in a 2snucleus then the sign of σ(D)ℓ(D) must be (−1)^{s+1}. This observation will be crucial for the automated procedure.
Results
We begin by summarising the key ingredients described in Methods that lead to the proposed automated procedure. The (computable) polynomial p _{ A,Z } might have positive and negative terms. On one hand, each term in the polynomial corresponds to a collection of nuclei, as described in the Section ‘The polynomial and circuits’. Each of these nuclei consists of circuits, some of which might be positive feedback loops. A term of sign (−1)^{s+1} in the polynomial necessarily corresponds to a collection of nuclei involving at least one positive feedback loop, see Section ‘The polynomial and circuits’.
On the other hand, the network can only be multistationary if the polynomial p _{ A,Z } has terms of opposite sign. Consequently, if the network is multistationary, then some terms of p _{ A,Z } have the sign (−1)^{s+1}. Therefore we conclude that (i) the network must contain positive feedback loops in order to be multistationary and (ii) the positive feedback loops underlying multistationarity can be found by considering the terms with sign (−1)^{s+1} in the polynomial. From these terms we identify the associated nuclei and extract the positive feedback loops in these nuclei.
It is an empirical observation that in most realistic or real reaction networks, the predominant sign of the coefficients of p _{ A,Z } is (−1)^{s}. Therefore, the number of terms to inspect, that is, the number of terms with sign (−1)^{s+1}, is usually low. For this reason, we call the sign (−1)^{s+1} the ‘wrong sign’.
Based on these considerations, we develop a procedure to extract the positive feedback loops that correspond to terms with the wrong sign in p _{ A,Z }. The procedure is based on the following steps. For a nonzero term with the wrong sign, say
(c is positive) consider the following edges from species to reactions
The 2snuclei corresponding to the term (17) must contain these edges. Therefore, we add to these edges all possible choices of s edges from reactions $\{r_{j_{1}},\dots,r_{j_{s}}\}$ to species $\{X_{i_{1}},\dots,X_{i_{s}}\}$ such that the resulting graph is a 2snucleus. We keep only the nuclei D for which the sign of σ(D)ℓ(D) is (−1)^{s+1}. The positive feedback loops in these nuclei are those that do contribute to the existence of multiple steady states. Indeed, if all these loops are broken, then the network cannot be multistationary because the polynomial will then only have terms of sign (−1)^{s}. We find these loops in the signed DSRgraph.
For example, consider the polynomial $p_{A,Z_{1}}$ in (14) and the DSRgraph shown in Figure 2. In this case, the rank of A is s=2, and hence we focus on the negative terms since (−1)^{s+1}=−1. These are γ _{2,2} γ _{3,3} and γ _{1,1} γ _{3,3}. The corresponding 4nuclei are depicted in Figure 2(A): there are two 4nuclei obtained as the union of the red circuit and one of the two blue circuits. The only positive feedback loop that appears is therefore the selfactivation positive feedback loop, and this is the only loop that is related to the observed multistationarity. The other two positive feedback loops (termed the spatial and the activation loop, respectively in Figure 2(B)) are therefore not relevant for the observed multistationarity.
We next state the automated procedure.
Automated procedure
The procedure to select positive feedback loops that contribute to multistationarity, applies to any reaction network defined as in (4), which fulfils the monotonicity criteria (9). This criterion states that the reaction rates v(x) must be strictly monotone in each variable x _{ j }, otherwise the influence matrix in not uniquely defined.
The procedure consists of the following steps.

1.
For a network with stoichiometric matrix A of rank s and influence matrix Z, compute p _{ A,Z } and select the terms with sign (−1)^{s+1}.

2.
Construct the DSRgraph. For each selected term of p _{ A,Z } with the wrong sign, determine the corresponding 2snuclei of the DSRgraph that have the wrong sign.

3.
For each of the nuclei, select the connected components that form positive feedback loops.
These steps have been implemented in Maple. The Maple script is provided as Additional file 1. The procedure might fail for practical reasons (such as lack of computational memory) if the number of species and reactions is too big. In our experience, this number depends heavily on the sparsity of the influence matrix [12].
It is worth emphasising that the procedure is meaningful only if we know that the input network is multistationary. We might apply the procedure to a nonmultistationary network and obtain a list of positive feedback loops, but these loops will lack proper interpretation.
Examples
We have applied the procedure to find positive feedback loops that are responsible for multistationarity in several reaction networks. These examples are also shown in the Additional file 1, together with some other systems such as the threesite phosphorylation system.
Ring1B/Bmi1 ubiquitination system
We consider an ODE model of histone H2A ubiquitination that involves the E3 ligases Ring1B and Bmi1 [3]. When degradation of species is not taken into account, the model has 10 species and 15 reactions. We let B and B $^{d}_{\textit {ub}}$ denote the protein Bmi1 in isolation and targeted for degradation by ubiquitination, respectively. The protein Ring1B is denoted by R, and R _{ ub }, R $_{\textit {ub}}^{a}$ , R $^{d}_{\textit {ub}}$ denote three different forms of selfubiquitinated R, with R $^{d}_{\textit {ub}}$ being the form targeted for degradation. Bmi1 and Ring1B form a complex Z, that also might be ubiquitinated, Z _{ ub }. Finally, Ring1B (either alone or in the complex Z) is responsible for the ubiquitination of the histone H2A, with states H, H _{ ub }.
The reactions describing the mechanism are [3]:
We let x _{1},…,x _{10} denote the concentrations of B, B $_{\textit {ub}}^{d}$ , R, R $_{\textit {ub}}^{d}$ , R _{ ub }, R $_{\textit {ub}}^{a}$ , Z, Z _{ ub }, H, H _{ ub }, respectively. The associated reaction rates are [3]:
where κ _{ i }>0 are constants.
Selfubiquitination of B is taken into account in the rate functions v _{7} and v _{11} for reactions r _{7} and r _{11}, respectively. These incorporate a positive influence from the reaction products. With these specific rate functions, the system is multistationary for some values of the reaction rate constants [3]. We apply the automated procedure to find positive feedback loops that are responsible for multistationarity and obtain the circuits depicted in Figure 3(A). In [3], it is postulated that selfubiquitination of Z and R are crucial steps for the emergence of multiple steady states, and we confirm the statement here.
Phosphorylation systems
We have analysed different networks of signal transmission based on phosphorylation. We consider models for sequentially distributed phosphorylation and dephosphorylation cycles and some modifications of these, see e.g. [10,40].
We consider first a twosite sequential phosphorylation cycle for a substrate S, where phosphorylation of the two sites is catalysed distributively by a kinase E, and dephosphorylation of the two sites uses different phosphatases F_{1}, F_{2}. Assuming a MichaelisMenten mechanism, the reaction network consists of the following reactions:
In S_{0}, S_{1}, S_{2} the subindex denotes the number of phosphorylated sites. With massaction kinetics, this system is multistationary for some choice of reaction rate constants [2,10]. However, the positive feedback loops that can account for the observed multistationarity are not trivial.
We apply the automated procedure and obtain the positive feedback loops given in Figure 3(B). The first of the two loops has two edges with negative labels. It implies that S_{0} and S_{1} enhance their own production. Indeed, because S_{0} and S_{1} both compete for the same kinase, an increase in S_{0} increases the rate of reaction r _{1}, which in turn lowers the amount of available enzyme E. This implies that reaction r _{7} slows down and hence S_{1} is consumed at a slower rate. The circuit closes through the reactions r _{4} and r _{6}, which shows that an increase in S_{1} implies an increase in S_{0}.
These type of loops are recurrent in phosphorylation motifs. It is worth emphasising that the loops do not have independent meaning outside the network. Another network with the same positive feedback loop might not be multistationary. For example, the second loop of Figure 3(B) also occurs in the MichaelisMenten mechanism $\mathrm {S}_{0}+\mathrm {E}\leftrightharpoons \text {ES}_{0}\rightarrow \mathrm {S}_{1}+\mathrm {E}$ , but these reactions alone are not multistationary.
Signalling cascades
We consider a 2layer cascade with an explicit positive feedback. The first layer is a phosphorylation cycle with kinase E, phosphatase F_{1}, and phosphorylated and unphosphorylated substrate S_{0}, S_{1}. The second layer has kinase S_{1}, phosphatase F_{2}, and phosphorylated and unphosphorylated substrate P_{0}, P_{1}. Assuming a MichaelisMenten mechanism, the reaction network consists of the following reactions:
This network is multistationary for some choice of reaction rate constants. The automated procedure finds three positive feedback loops, as shown in Figure 3(C). The first loop is expected and appears because the product of the second layer, P_{1}, activates the kinase of the first layer, E. The other two loops are similar to those in Figure 3(C).
Apoptosis
We finally consider a basic model of caspase activation for apoptosis [9]:
With massaction kinetics, this network is multistationary for some choice of reaction rate constants [9] and has two relevant positive feedback loops, Figure 3(D). The second loop consists of two species, each with positive influence on a reaction rate, while at the same time, decreasing the amount of the other.
Analysis of the Biomodels database
We investigated the models in the database PoCaB [41], which consists of 365 models from the publicly available database Biomodels [36] (see also the page http://www.ebi.ac.uk/biomodelsmain/). The database PoCaB contains precomputed stoichiometric matrices, massaction exponent matrices, and kinetic data from the selected models.
In a previous paper [12] we extracted influence matrices of the reported kinetics. Of the 365 networks, 323 have strictly monotone kinetics such that the influence matrix is well defined. On these 323 networks we ran the injectivity method and ended up with a total of 64 noninjective networks, excluding 27 very large networks where the injectivity method failed (as described in [12]). Noninjectivity is a prerequisite for being multistationary.
We applied the automated procedure on the 64 networks to determine the positive feedback loops. We obtained a total of 341 different positive feedback loops with size distribution shown in Figure 4(C). Most loops involve only 2 or 3 species (112 and 108 loops out of 341, respectively). In Figure 4(AB), we show the positive feedback loops involving one or two species and conclude that all possible types occur in the database. However, for two species, their frequencies vary from 9 to 35 (out of 108), indicating that the motifs are not equally represented in the database (Pearson’s chisquare test, p<0.0005). For one species, there appears to be no difference (p=0.28). We show in Table 1 the most frequent positive feedback loops for different number of species nodes.
Discussion
We have presented an automated procedure to find the positive feedback loops in a multistationary network that are contributing to the multistationarity. The procedure relies on structural and qualitative features of the network together with a kinetics, and it is insensitive to the specific form of the rate functions. Only positive feedback loops that are contributing to the multistationarity of the network are reported, excluding those positive feedback loops that are not relevant.
Whether a loop is relevant or not, depends on the entire DSRgraph of the network (that is, the reactions and the influence matrix) and not just on the loop itself. In this sense, being a positive feedback loop related to an observed multistationarity, is a context or network dependent property. This fact has also been observed in [3335]. In these papers, it is shown that the existence of a positive feedback loop fulfilling a certain extra condition is a requirement for multistationarity to occur. Specifically, the loop must either intersect another positive feedback loop in a specific way (called an StoR intersection) or fulfil a certain condition on the labels (called an scycle). The first possibility is network dependent. It is worth mentioning that there can be positive feedback loops that are scycles or that intersect another positive feedback loop in an StoR intersection without being relevant for the observed multistationarity. This is the case for most of the examples presented here. For example, the apoptosis network has another positive feedback loop, Y →r _{14}→IAP→r _{3}→Y (all edges are positive), which intersects the loop on the right side in Figure 3(D) in an StoR intersection.
The property of network dependence is further illustrated in Figure 2(A)(B), where the procedure is applied to the reaction network in Figure 1 with influence matrix given by (10). The network models translocation and phosphorylation of a cyclin dependent kinase X in the onset of mitosis. Only one of the three positive feedback loops shown in Figure 2(B) can be associated with multistationarity in the model, namely the selfactivating loop that stimulates the creation of phosphorylated X $^{*}_{\text {nuc}}$ in the nucleus. In [1], it is argued by different means than in this paper, that the spatial redistribution of the cyclin dependent kinase is important for creating the observed bistability. In contrast, our results suggest that the observed bistability is due to the selfactivation loop of the phosphorylated complex in the nucleus.
The presented procedure cannot establish whether a reaction network is multistationary or not. Other means will here be required. If the procedure is applied to a reaction network that might or might not be multistationary, then the absence of positive feedback loops implies that the network cannot be multistationary, whereas the presence of positive feedback loops leaves the question open.
Whether a positive feedback loop found by the procedure is important in a biological or experimental context, is not addressed in this paper, but must be established in other ways, for example by experimental verification. A reaction network might only be multistationary for very specific choices of reaction rates, which might not be relevant in a particular experimental setting. As the procedure is parameter independent, any such verification and subsequent interpretation is beyond the scope of the procedure.
Conclusions
It is well known that multistationarity requires the existence of a positive feedback loop. However, positive feedback loops are abundant in most biochemical reaction networks and it remained unclear how to select the positive feedback loops that could be accounted for underlying the observed multistationarity. In this work we have proposed an automatised method to determine the relevant positive feedback loops. The method is based on the mathematical arguments underlying the proof that positive feedback loops are necessary for multistationarity to occur. We have tested and illustrated the method with several reaction networks. An implementation of the method in Maple is provided as supplementary material.
References
 1
Santos SD, Wollman R, Meyer T, Ferrel JE Jr.Spatial positive feedback at the onset of mitosis. Cell. 2012; 149(7):1500–13.
 2
Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J. Cell Biol. 2004; 164:353–9.
 3
Nguyen LK, MuñozGarcía J, Maccario H, Ciechanover A, Kolch W, Kholodenko BN. Switches, excitable responses and oscillations in the ring1b/bmi1 ubiquitination system. PLoS Comput Biol. 2011; 7(12):1002317.
 4
Zhdanov VP. Bistability in gene transcription: interplay of messenger RNA, protein, and nonprotein coding RNA. Biosystems. 2009; 95(1):75–81.
 5
Legewie S, Bluthgen N, Schäfer R, Herzel H. Ultrasensitization: switchlike regulation of cellular signaling by transcriptional induction. PLoS Comput Biol. 2005; 1(5):54.
 6
Palani S, Sarkar CA. Synthetic conversion of a graded receptor signal into a tunable, reversible switch. Mol Sys Biol. 2011; 7:480.
 7
Liu D, Chang X, Liu Z, Chen L, Wang R. Bistability and oscillations in gene regulation mediated by small noncoding rnas. PLoS ONE. 2011; 6(3):17029.
 8
Harrington HA, Feliu E, Wiuf C, Stumpf MPH. Cellular compartments cause multistability in biochemical reaction networks and allow cells to process more information. Biophys J. 2013; 104:1824–31.
 9
Eissing T, Conzelmann H, Gilles ED, Allgower F, Bullinger E, Scheurich P. Bistability analyses of a caspase activation model for receptorinduced apoptosis. J Biol Chem. 2004; 279(35):36892–97.
 10
Feliu E, Wiuf C. Enzymesharing as a cause of multistationarity in signalling systems. J. R. S. Interface. 2012; 9(71):1224–32.
 11
Conradi C, Flockerzi D, Raisch J, Stelling J. Subnetwork analysis reveals dynamic features of complex (bio)chemical networks. Proc Nat Acad Sci. 2007; 104(49):19175–80.
 12
Feliu E, Wiuf C. A computational method to preclude multistationarity in networks of interacting species. Bioinformatics. 2013; 29:2327–34.
 13
Wiuf C, Feliu E. Powerlaw kinetics and determinant criteria for the preclusion of multistationarity in networks of interacting species. SIAM J Appl Dyn Syst. 2013; 12:1685–721.
 14
Conradi C, Flockerzi D. Switching in mass action networks based on linear inequalities. SIAM J Appl Dyn Syst. 2012; 11(1):110–34.
 15
Craciun G, Feinberg M. Multiple equilibria in complex chemical reaction networks. I, The injectivity property. SIAM J Appl Math. 2005; 65(5):1526–46.
 16
Pérez Millán M, Dickenstein A, Shiu A, Conradi C. Chemical reaction systems with toric steady states. Bull Math Biol. 2012; 74:1027–65.
 17
Feinberg M. Chemical reaction network structure and the stability of complex isothermal reactors I. The deficiency zero and deficiency one theorems. Chem Eng Sci. 1987; 42(10):2229–68.
 18
OteroMuras I, Banga JR, Alonso AA. Characterizing multistationarity regimes in biochemical reaction networks. PLoS One. 2012; 7(7):39194.
 19
Domijan M, Kirkilionis M. Bistability and oscillations in chemical reaction networks. J Math Biol. 2009; 59:467–501.
 20
Ellison P, Feinberg M, Ji H, Knight D. Chemical Reaction Network Toolbox. Version 2.2. 2012, Available online at http://www.crnt.osu.edu/CRNTWin.
 21
Donnell P, Banaji M, Marginean A, Pantea C. CoNtRol: an open source framework for the analysis of chemical reaction networks. Bioinformatics. 2014; 30(11):1633–4.
 22
Tyson JJ, Chen KC, Novak B. Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003; 15:221–31.
 23
Rand DA, Shulgin BV, Salazar D, Millar AJ. Design principles underlying circadian clocks. J R Soc Interface. 2004; 119–130:2874–80.
 24
Joshi B, Shiu A. Atoms of multistationarity in chemical reaction networks. J Math Chem. 2013; 51(1):153–78.
 25
Feinberg M, Horn FJM. Chemical mechanism structure and the coincidence of the stoichiometric and kinetic subspaces. Arch Rational Mech Anal. 1977; 66(1):83–97.
 26
Angeli D, De Leenheer P, Sontag E. Graphtheoretic characterizations of monotonicity of chemical networks in reaction coordinates. J Math Biol. 2010; 61:581–616.
 27
Horn FJM, Jackson R. General mass action kinetics. Arch Rational Mech Anal. 1972; 47:81–116.
 28
Jacob F, Monod J. Genetic regulatory mechanisms in the synthesis of proteins. J Mol Biol. 1961; 3:318–56.
 29
Thomas R. On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations. Springer Syne. 1981; 9:180–3.
 30
Soulé C. Graphic requirements for multistationarity. ComPlexUs. 2003; 1:123–33.
 31
Gouze JL. Positive and negative circuits in dynamical systems. J Biol Syst. 1998; 6:11–15.
 32
Kaufman M, Soulé C, Thomas R. A new necessary condition on interaction graphs for multistationarity. J Theor Biol. 2007; 248:675–85.
 33
Craciun G, Feinberg M. Multiple equilibria in complex chemical reaction networks. II, The speciesreaction graph. SIAM J Appl Math. 2006; 66(4):1321–38.
 34
Banaji M, Craciun G. Graphtheoretic criteria for injectivity and unique equilibria in general chemical reaction systems. Adv Appl Math. 2010; 44:168–84.
 35
Banaji M, Craciun G. Graphtheoretic approaches to injectivity and multiple equilibria in systems of interacting elements. Commun Math Sci. 2009; 7(4):867–900.
 36
Li C, Donizelli M, Rodriguez N, Dharuri H, Endler L, Chelliah V, et al.BioModels Database: an enhanced, curated and annotated resource for published quantitative kinetic models. BMC Syst Biol. 2010; 4:92.
 37
Feinberg M. Lectures on chemical reaction networks. 1980. Available online at http://www.crnt.osu.edu/LecturesOnReactionNetworks.
 38
Feliu E, Wiuf C. Preclusion of switch behavior in reaction networks with massaction kinetics. Appl Math Comput. 2012; 219:1449–67.
 39
Müller S, Feliu E, Regensburger G, Conradi C, Shiu A, Dickenstein A. Sign conditions for the injectivity of polynomial maps in chemical kinetics and real algebraic geometry. Foundations Comput Math. 2015. http://link.springer.com/article/10.1007/s1020801492393. Published online Jan 5.
 40
Gunawardena J. Distributivity and processivity in multisite phosphorylation can be distinguished through steadystate invariants. Biophys J. 2007; 93:3828–34.
 41
Samal SS, Errami H, Weber A. Pocab: A software infrastructure to explore algebraic methods for biochemical reaction networks. Lect Notes Comput Sci. 2012; 7442:294–307.
Acknowledgements
EF was supported by project MTM201238122C0301/FEDER from the Ministerio de Economía y Competitividad, Spain, the Danish Research Council and the Lundbeck Foundation (Denmark). CW was supported by the Danish Research Council, the Lundbeck Foundation (Denmark) and the Carlsberg Foundation (Denmark). This paper was completed while the authors were visiting the Isaac Newton Institute in Cambridge, UK.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
Both authors contributed equally to all aspects of this work. Both authors read and approved the final manuscript.
Additional file
Additional file 1
Maple Script contains the automated procedure as well as examples.
Rights and permissions
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Bistability
 DSRgraph
 Feedback loop
 Injectivity
 Interaction graph
 Reaction network