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:
$$ r_{j}\colon \sum_{i=1}^{n} \alpha_{ij} X_{i} \rightarrow \sum_{i=1}^{n}\beta_{ij} X_{i}, \qquad j=1,\dots,m $$
((3))
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:
$$a_{ij} = \beta_{ij}\alpha_{ij}, $$
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
$$ v_{j}(x)=\kappa_{j} \,x_{1}^{\alpha_{1j}}\cdot \dots \cdot x_{n}^{\alpha_{nj}},\quad x\in\Omega_{v}, $$
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:
$$ \dot{x} = Av(x),\qquad x\in \Omega_{v}. $$
((4))
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:
$$\begin{array}{*{20}l} \dot{x}_{1} & = \kappa_{1} x_{1} + \kappa_{2} x_{2} \\ \dot{x}_{2}& = \kappa_{1} x_{1}  \kappa_{2} x_{2}  \frac{x_{2} (x_{2}+x_{3})^{4} }{K^{4}+(x_{2}+x_{3})^{4}} + \kappa_{4} x_{3} \\ \dot{x}_{3} & = \frac{x_{2} (x_{2}+x_{3})^{4} }{K^{4}+(x_{2} +x_{3})^{4}}  \kappa_{4} x_{3}, \end{array} $$
((5))
where κ
_{1},…,κ
_{4},K>0 are parameters. It takes the form (4) with
$$ A=\left(\begin{array}{rrrr} 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{array} \right), $$
((6))
$$ v(x)=\left(\kappa_{1} x_{1},\kappa_{2} x_{2},\frac{x_{2} (x_{2}+x_{3})^{4} }{K^{4}+(x_{2}+x_{3})^{4}},\kappa_{4} x_{3}\right), $$
((7))
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:
$$ v(x)=\left(\kappa_{1} x_{1},\kappa_{2} x_{2}, \frac{{x_{2}^{5}} }{K^{4}+{x_{2}^{4}}},\kappa_{4} x_{3}\right). $$
((8))
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])
$$ \mathcal C_{0}= \left(x_{0}+ \text{im}(A)\right) \cap\mathbb{R}_{\geq 0}, $$
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
$$ z_{ij} = \begin{cases} \gamma_{ij} & \textrm{if }v_{j}(x)\textrm{increases in }x_{i}\\ \gamma_{ij} & \textrm{if }v_{j}(x)\textrm{decreases in }x_{i} \\ 0 & \textrm{if }v_{j}(x)\textrm{is constant in }x_{i}, \end{cases} $$
((9))
where γ
_{
ij
} are symbolic variables.
The influence matrices associated with the two reaction rate vectors in (7) and (8) are given by
$$ Z_{1}=\left(\begin{array}{cccc} \gamma_{1,1} & 0 & 0 & 0 \\ 0 & \gamma_{2,2} & \gamma_{2,3} & 0 \\ 0 & 0 & \gamma_{3,3} & \gamma_{3,4} \end{array} \right), $$
((10))
and
$$ Z_{2}=\left(\begin{array}{cccc} \gamma_{1,1} & 0 & 0 & 0 \\ 0 & \gamma_{2,2} & \gamma_{2,3} & 0 \\ 0 & 0 & 0& \gamma_{3,4} \end{array} \right), $$
((11))
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
$$ \sigma(D)\ell(D) =(1)^{a_{2}} \prod_{i=1}^{a} \ell(C_{i}). $$
((12))
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
$$\text{sign}(\sigma(D)\ell(D)) = (1)^{a_{2}+a} = (1)^{a_{1}+2a_{2}} = (1)^{a_{1}}. $$
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
$$ \text{sign}(\sigma(D)\ell(D))=(1)^{s} $$
((13))
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
$$p_{A,Z}=\det(\widetilde{M}), $$
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
$$\begin{array}{*{20}l} p_{A,Z_{1}} & = \gamma_{2,2}\gamma_{3,3} \gamma_{1,1}\gamma_{3,3} \\ & \gamma_{2,2} \gamma_{3,4} + \gamma_{1,1}\gamma_{2,3} + \gamma_{1,1}\gamma_{3,4}. \end{array} $$
((14))
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,
$$\begin{array}{*{20}l} p_{A,Z_{2}} & = \gamma_{2,2} \gamma_{3,4} + \gamma_{1,1}\gamma_{2,3} + \gamma_{1,1}\gamma_{3,4}. \end{array} $$
((15))
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
$$ p_{A,Z} = \sum_{I,J\subseteq \{1,\dots,n\}} \sum_{D \in D_{s}(I,J)} \sigma(D)\ell(D), $$
((16))
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.