Families of Toric Chemical Reaction Networks

TL;DR

This paper investigates families of chemical reaction networks with toric steady states, providing algebraic tools to analyze their multistationarity, model distinguishability, and implications for experimental design, exemplified by phosphorylation networks.

Contribution

It introduces an algorithmic construction of reaction network families with toric steady states and uses matroids to analyze their steady state behavior and distinguishability.

Findings

Multistationarity propagates through network families.

Matroids encode steady state dependencies and model distinguishability.

Necessary conditions for network family distinguishability are established.

Abstract

We study families of chemical reaction networks whose positive steady states are toric, and therefore can be parameterized by monomials. Families are constructed algorithmically from a core network; we show that if a family member is multistationary, then so are all subsequent networks in the family. Further, we address the questions of model selection and experimental design for families by investigating the algebraic dependencies of the chemical concentrations using matroids. Given a family with toric steady states and a constant number of conservation relations, we construct a matroid that encodes important information regarding the steady state behaviour of the entire family. Among other things, this gives necessary conditions for the distinguishability of families of reaction networks with respect to a data set of measured chemical concentrations. We illustrate our results using multi-site phosphorylation networks.

Equations

$$\begin{array}{rl}& E+S_0\rightleftharpoons E{S}_0\to E+S_1,\\ & F+S_1\rightleftharpoons F{S}_1\to F+S_0.\end{array}$$

$$\displaystyle\begin{split}E+S_{0}\rightleftharpoons ES_{0}\rightarrow E+S_{1}&{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\rightleftharpoons ES_{1}\rightarrow E+S_{2}},\\ &{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}F+S_{2}\rightleftharpoons FS_{2}\rightarrow}F+S_{1}\rightleftharpoons FS_{1}\rightarrow F+S_{0},\end{split}$$

$$\displaystyle\begin{split}E+S_{0}\rightleftharpoons ES_{0}&{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\rightleftharpoons ES_{1}\rightarrow}E+S_{2},\\ F+S_{2}&{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\rightleftharpoons FS_{2}}\rightarrow FS_{1}\rightarrow F+S_{0}.\end{split}$$

$$E+S\rightleftharpoons ES\to E+P$$

$$\sum _l{\alpha }_{il}{X}_l\stackrel{\kappa }{\to }\sum _l{\alpha }_{jl}{X}_l;$$

$$x^{{\alpha }_i}=x_1^{{\alpha }_{i1}}\cdots x_n^{{\alpha }_{in}}.$$

$$x_1{x}_3\underset{{\kappa }_2}{\overset{{\kappa }_1}{\rightleftharpoons }}{x}_5\stackrel{{\kappa }_3}{\to }{x}_1{x}_4,$$

$$x_2{x}_4\underset{{\kappa }_5}{\overset{{\kappa }_4}{\rightleftharpoons }}{x}_6\stackrel{{\kappa }_6}{\to }{x}_2{x}_3.$$

$$\frac{dx}{dt}={\alpha }^T{A}_{\kappa }{x}^{\alpha }$$

$$V(f_1,\dots ,f_r)=\{x\in {\mathbb{C}}^n|f_1(x)=\cdots =f_r(x)=0\}.$$

$$\frac{dx}{dt}={\alpha }^T{A}_{\kappa }{x}^{\alpha }$$

$$X_A=\overline{\{(t^{a_1},\dots ,t^{a_n}):t\in ({\mathbb{C}}^{\ast }{)}^d\}}\subset {\mathbb{C}}^n.$$

$$I=⟨x^{c^{+}}-x^{c^{-}}|c\in \ker (A)⟩,$$

$$I=⟨c_i{x}^{b_i^{+}}-k_i{x}^{b_i^{-}}|i=\{1,\dots ,\nu \},c_i,k_i\in \mathbb{R}(\kappa )⟩\subseteq \mathbb{R}(\kappa )[x_1,\dots ,x_n].$$

$${\psi }_A:({\mathbb{C}}^{\ast }{)}^d\to {\mathbb{C}}^n\mathrm{where}t\mapsto (x_1^{\ast }{t}^{a_1},\dots ,x_n^{\ast }{t}^{a_n})$$

$${\varphi }_A:\mathbb{C}[x_1,\dots ,x_n]$$

$${\varphi }_A(x_i)=x_i^{\ast }{t}^{a_i}$$

$$I=⟨$$

$$-{\kappa }_1{x}_1{x}_3+{\kappa }_2{x}_5+{\kappa }_6{x}_3,{\kappa }_3{x}_5-{\kappa }_4{x}_2{x}_4+{\kappa }_5{x}_6,$$

$${\kappa }_1{x}_1{x}_3-({\kappa }_2+{\kappa }_3)x_5,{\kappa }_4{x}_2{x}_4-({\kappa }_5+{\kappa }_6)x_6⟩,$$

$$=⟨$$

$$A=\left(\begin{array}{cccccc}1& 0& -1& 0& 0& 0\\ 0& 1& 0& -1& 0& 0\\ 0& 0& -1& -1& -1& -1\end{array}\right).$$

$$E+S_0$$

$$F+S_1$$

$$E+S_0$$

$$\displaystyle{\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}F+S_{2}\rightleftharpoons FS_{2}\rightarrow}F+S_{1}$$

$$E+S_0$$

$$F+P$$

$$E+S_0\rightleftharpoons E{S}_0$$

$$F+P$$

$$X_i+X_{i+1}\to 2X_{i+1},$$

$$X_i\rightleftharpoons \varnothing \text{for }i=1,\dots ,N\text{ and setting }{X}_{N+1}=X_1.$$

$$\pi :\mathbb{R}[{\kappa }_1,\dots ,{\kappa }_{m+m^{\prime }}][x_1,\dots ,x_n]$$

$$\pi ({\kappa }_i)$$

$$\mu :\mathcal{C}\times \mathcal{C}$$

$$\mu \left((C_i,C_j)\right)$$

Full text

1 Introduction

Many of the fundamental processes in biological cells can be described by a set of interlinked chemical reactions. Prominent examples of cellular processes regulated via biochemical interactions include immune response [1], cell signalling [2], cell death [3, 4], and toxin formation [5]. For this reason the study of chemical reaction networks forms a central part of algebraic systems biology [6, 7, 8, 9, 10]. One approach focuses on the long term behaviour of networks by investigating their steady states and the relation of the number and stability of steady states to the network structure [10, 11, 13]. In this paper we investigate the positive steady states for algorithmically constructed reaction networks, which we call families, for which the positive steady states may be parameterized by monomials (i.e. polynomials with a single term). Further, we use the algebraic dependencies of the variables representing the chemical concentrations to investigate experimental design and model identification for entire families.

Families of networks are formally defined in Definition 3.2. To obtain an intuition for what could be described as a family, we introduce phosphorylation networks [13, 19]. Phosphorylation is a vital signalling process in biochemistry and it is one of the most widely studied protein modifications. During phosphorylation a phosphoryl group ($PO_{3}^{-}$) is added to an organic molecule which acts as a substrate. The chemical reaction is catalysed by enzymes and often many phosphoryl groups can be added to the same substrate. This process is known as multisite phosphorylation. There are two limiting mechanisms which have been investigated in the literature and which will serve as the running examples in this paper, namely, distributive phosphorylation and processive phosphorylation. In the distributive system the enzyme unbinds from the substrate after every time a $PO_{3}^{-}$ is added and in the processive mechanism the enzyme stays bound until the substrate is fully phosphorylated. The desphosphorylation mechanisms work analogously. The distributive and processive mechanisms for one site substrates follow the same reaction scheme

[TABLE]

However, on a two-site substrate the mechanisms differ, namely, the distributive mechanism is given by

[TABLE]

and the processive system is

[TABLE]

It is clear that the constructions for both the distributive and the processive system can be extended to an $N$-site substrate and that two different procedures are needed to do so. Hence, we can consider these two mechanisms to be in two distinct families of networks. In this paper we study graph theoretic constructions which allows us to rigorously identify families of networks and also construct them from a base graph. Most importantly, though, we investigate which steady state properties are conserved throughout a family.

One of the central goals of this paper is to develop criteria to establish which families have members with multiple positive steady states, so-called multistationary networks. Confirming that a network is multistationary and finding the associated parameter regions is highly nontrivial; a range of different approaches have been applied previously (see [20] for a survey). In particular, monomial positive steady state parameterizations have proved fruitful due to their relations to toric varieties which are well understood in algebraic geometry [18].

Previous work relating to the concept of families in this paper considers so-called atoms of multistationarity, which are the smallest multistationary subnetworks which can induce multistationarity in their parent networks [25]. Network properties resulting from the gluing of networks are investigated in [26]. Other network modifications which preserve or destroy multistationarity are studied in [27]. Recent results extend the techniques for identifying multistationarity to highly structured networks [23] and networks with intermediate species [24, 42]. In this paper we develop a concept of families of networks which unifies the notions of highly structured networks, embedded subnetworks as defined in [25] (i.e. distributive systems), and networks with intermediate species [42] (i.e. processive systems). We show in Theorem 5.2 that if a member of a family is multistationary then so are all larger networks obtained by the recursive construction of Definition 3.2.

Going beyond multistationarity, we also investigate the necessary conditions for model rejection among members of a family if only limited steady state data is available. In particular, in Section 4, we encode algebraic dependencies between the variables using a combinatorial object called an algebraic matroid [29, 30]. From the algebraic matroid we find binomial relations which have to be satisfied by the chemical concentrations at any positive steady state, so-called steady state invariants. The results in Section 4.1 extend the previous research of [33, 34, 32] and the application of matroids for experimental design presented in [6, 34]. Using these previous results, we give novel necessary conditions for the distinguishability of two members of a family of reaction networks with respect to a data set of measured chemical concentrations. Consequently, we can also prescribe which species to measure to be able to reject a family member.

In summary, the biochemical questions we would like to address are:

What conditions are needed to define families of chemical reaction networks and what is the relation between their steady states? (Section 3) 2. 2.

Can we use the family construction in model selection or parameter estimation for the entire family? (Section 4) 3. 3.

Can we find conditions such that multistationarity of one family member implies multistationarity for all subsequent members? (Section 5)

These motivating questions from chemistry translate into the following algebraic questions which we answer using techniques from toric geometry and matroid theory:

What are the relations between the toric varieties defined by recursively constructed reaction graphs? 2. 2.

What is the connection between the circuit polynomials of matroids associated to different family members? 3. 3.

What is the relation between the positive parts of the steady state varieties of subsequent family members?

This paper is organised as follows. Section 2 introduces chemical reaction networks and relevant definitions from toric geometry and matroid theory. In Section 3 we give a rigorous definition of a family of reaction network graphs and some preliminary results. In Section 4 we focus on biochemical and algebraic question 2 using matroid theory. We also introduce some new terminology which simplifies the proofs in the remainder of the paper. In Section 5 we prove the main result on multistationarity using the matroidal language developed in the previous section. We discuss our results and suggest further directions in Section 6.

2 Background

In this section we briefly introduce aspects of chemical reaction network theory and discuss essential notions of algebraic geometry and matroid theory.

2.1 Chemical Reaction Network Theory

Informally a chemical reaction network (CRN) can be described by a multiset $\mathfrak{N}=\{\mathcal{S},\mathcal{C},\mathcal{R}\}$, where $\mathcal{S}$ is the set of species, $\mathcal{C}$ is the set of linear combinations of species (complexes) and $\mathcal{R}$ is the set of reactions.

Example 2.1** (Michaelis-Menten).**

The set $\mathcal{S}$ of chemical species present in the network

[TABLE]

is defined by $\mathcal{S}=\{E,S,ES,P\}$. The complexes, which are linear combinations of species, are $\mathcal{C}=\{E+S,ES,E+P\}$. The reaction set is $\mathcal{R}=\{E+S\to ES,ES\to E+S,ES\to E+P\}$.

The multiset $\mathfrak{N}$ defines a directed graph (digraph) $\mathcal{G}$ whose vertex set is $\mathcal{C}$ and whose edge set is defined by the reaction set $\mathcal{R}$. The reaction from complex $C_{i}$ to $C_{j}$ is an element of the reaction set $\mathcal{R}$ if and only if there is a directed edge $C_{i}\to C_{j}$ in $\mathcal{G}$. Let $X_{l}\in\mathcal{S}$ and $\{\alpha_{il}\}\in\mathbb{Z}_{\geq 0}$. A reaction from complex $C_{i}=\sum_{l}\alpha_{il}X_{l}$ to $C_{j}=\sum_{l}\alpha_{jl}X_{n}$, with rate constant $\kappa$ is written as

[TABLE]

the constants $\alpha_{il}$ are called the stoichiometric coefficients of the complex $C_{i}$. Let the reaction vector for the $\ell^{\text{th}}$ reaction $C_{i}\to C_{j}$ be $r_{\ell}=\alpha_{j}-\alpha_{i}$ where $\alpha_{i}$, $\alpha_{j}$ are the column vectors of the stoichiometric coefficients of the complexes $C_{i}$ and $C_{j}$. The $n\times m$ matrix of all reaction vectors $\Gamma=(r_{1},\dots,r_{m})$ is called the stoichiometric matrix. The rate constant $\kappa$ assigns a weight to each edge of the digraph $\mathcal{G}$, making $\mathcal{G}$ a weighted digraph. Definition 2.2 gives the description of a CRN which we are going to adopt for the remainder of this paper. While this is not the standard definition of a CRN, but equivalent to it, it allows us to focus more on the graphical structure of CRNs.

Definition 2.2**.**

A chemical reaction network is a weighted directed graph $\mathcal{G}=(\mathcal{C},\mathcal{R})$ with vertex set $\mathcal{C}$, edge set $\mathcal{R}$ and edge weights $\kappa=(\kappa_{1},\dots,\kappa_{m})^{T}$.

To connect the graphical structure of a CRN to its dynamical properties a description of reaction kinetics is needed. Previous work introduced a number of reaction laws such as mass action, rational function kinetics, Michaelis-Menten kinetics or Hill function kinetics [14]. In this paper we use mass action kinetics which assigns a monomial to each complex in the network. Let the chemical concentration of chemical species $X_{n}$ be $x_{n}$, then the monomial for complex $C_{i}$ is obtained by

[TABLE]

Hence, using the representation of complexes as monomials we represent $\mathcal{C}$ as a vector of monomials $x^{\alpha}=(x^{\alpha_{1}},\dots,x^{\alpha_{m}})^{T}$.

Example 2.3**.**

Revisiting the one site phosphorylation network (1), we map each species to its concentration $E\to x_{1},\>F\to x_{2},\>S_{0}\to x_{3},\>S_{1}\to x_{4},\>ES_{0}\to x_{5},\>FS_{1}\to x_{6}$ and introduce a vector of rate constants $\kappa=(\kappa_{1},\dots,\kappa_{6})^{T}$. Then, the reaction network is represented by the weighted digraph

[TABLE]

The dynamics of the network can be expressed in terms of the network structure and the stoichiometric coefficients as a set of ordinary differential equations

[TABLE]

where $A_{\kappa}$ is the negative weighted graph Laplacian of $\mathcal{G}$, $\alpha^{T}$ is the matrix of stoichiometric coefficients and $x^{\alpha}$ is a vector of monomials. An alternative representation of equation (6) assigns a monomial $\kappa_{\ell}x^{\alpha_{\ell}}$ to the $\ell^{\text{th}}$ reaction in the network to build a vector $R(x)=(\kappa_{1}x^{\alpha_{1}},\dots,\kappa_{m}x^{\alpha_{m}})^{T}$. Then, the dynamical system (6) is given by $dx/dt=\Gamma R(x)$, where $\Gamma$ is the stoichiometric matrix as above.

The left kernel of the stoichiometric matrix, $\Gamma$, is of biochemical importance as it describes conservation relations. Conservation relations induce linear relations between the variables. Informally we say that a set of species is conserved if their total concentration, $c$, is constant. In particular, suppose $z\in\ker\left(\Gamma^{T}\right)$, then $z^{T}(dx/dt)=0$ which implies $z^{T}x=c\in\mathbb{R}_{>0}$; this latter equation is then a conservation relation.

Remark 2.4**.**

Given a CRN, $\mathfrak{N}$, with a maximum of $d$ independent conservation relations, the linear subspace defined by the conservation relations, often referred to as compatibility class, may be compactly written as $Zx-c=0$, where $c\in\mathbb{R}_{\geq 0}^{d}$. The rows of $Zx-c$ define the subspace and the matrix $Z=(z_{1},\dots,z_{d})^{T}$ represents the conservation relations.

Example 2.5**.**

In the one site phosphorylation system (1) (also Example 2.3), the enzyme $E$ is conserved, but can exist in two states, the free state $E$ and the bound state $ES_{1}$. Thus, the total concentration $c_{1}=x_{1}+x_{3}$ is conserved. In total there are three independent conservation relations as the total mass of $E$, $F$, and the substrate is conserved respectively.

We now proceed to defining the steady states of a chemical reaction network [19].

Definition 2.6**.**

A vector $x^{*}\in\mathbb{C}^{n}$ is called a steady state of a CRN if $\alpha^{T}A_{\kappa}\left(x^{*}\right)^{\alpha}=0$. A positive steady state is a steady state such that $x^{*}\in\mathbb{R}_{>0}^{n}$ and $\alpha^{T}A_{\kappa}\left(x^{*}\right)^{\alpha}=0$.

Often one is interested in whether a chemical reaction network can have multiple positive steady states for a given set of rate constants $\kappa$ and total concentrations $c$.

Definition 2.7**.**

A chemical reaction network is multistationary if there exists a set of parameters $\{\kappa_{1},\dots,\kappa_{m}\}$ such that $\alpha^{T}A_{\kappa}\left(x^{*}\right)^{\alpha}=\alpha^{T}A_{\kappa}\left(y^{*}\right)^{\alpha}=0$ and $x^{*}-y^{*}\in\ker(Z)$ for two distinct positive steady states $x^{*}$ and $y^{*}$.

2.2 Algebraic Geometry and Toric Steady States

In this subsection we briefly introduce the notion of toric varieties from algebraic geometry, and their connections to chemical reaction networks. For an introduction to toric varieties we refer the reader to [18, 38].

Given polynomial equations $f_{1}(x_{1},\dots,x_{n}),\dots,f_{r}(x_{1},\dots,x_{n})$ in the polynomial ring $\mathbb{C}[x_{1},\dots,x_{n}]$ we define the algebraic variety

[TABLE]

Note that, by definition, $g(x)=0$ for any $x\in V(f_{1},\dots,f_{r})$ and any polynomial $g$ in the ideal $\langle f_{1},\dots,f_{r}\rangle=\{\sum_{i}h_{i}f_{i}\;|\>h_{i}\in\mathbb{C}[x_{1},\dots,x_{n}]\}\subset\mathbb{C}[x_{1},\dots,x_{n}]$; hence we may also associate a variety $V(I)=V(f_{1},\dots,f_{r})$ to an ideal $I=\langle f_{1},\dots,f_{r}\rangle$.

As in §2.1 we will consider a chemical reaction network $\mathcal{G}=(\mathcal{C},\mathcal{R})$ defining a ODE model

[TABLE]

as in (6). The set of all steady states of this system of ODEs (in an affine space $\mathbb{C}^{n}$) defines an algebraic variety generated by the steady state ideal $I=\langle\alpha^{T}A_{\kappa}x^{\alpha}\rangle\subseteq R=\mathbb{C}[x_{1},\dots,x_{n}]$. The ring $R$ is generated by the chemical concentrations and defined over the field $\mathbb{C}$.

A toric variety is an algebraic variety which can be described as the image of a monomial map. Work over the field $\mathbb{C}$ and fix an integer $d\times n$-matrix $A$, with columns $a_{1},a_{2},\ldots,a_{n}$, and rank $d$. A given column vector $a_{i}$ defines a (Laurent) monomial $t^{a_{i}}=t_{1}^{a_{1i}}t_{2}^{a_{2i}}\cdots t_{d}^{a_{di}}$ where $t\in(\mathbb{C}^{*})^{d}$, and $\mathbb{C}^{*}=\mathbb{C}\backslash\{0\}$ denotes the non-zero elements of the field (this is sometimes called the complex torus). The affine toric variety over $\mathbb{C}$ defined by $A$ is

[TABLE]

In words, the variety $X_{A}\,$ is the (Zariski) closure in $\mathbb{C}^{n}$ of the image of the monomial map defined by $A$. The defining implicit equations for $X_{A}$ can always be chosen to be binomials, that is $X_{A}=V(I)$ where I is an ideal defined by binomial equations. Specifically by [12, Corrollary 4.3] we have that

[TABLE]

where $c^{+}_{i}$ is equal to $c_{i}$ if $c_{i}>0$ and [math] otherwise, and where $c^{-}_{i}$ is equal to $|c_{i}|$ if $c_{i}<0$ and [math] otherwise. The ideal $I$ above is, additionally, a prime ideal; we say an ideal $I\subset\mathbb{C}[x_{1},\dots,x_{n}]$ is prime if whenever $f\cdot g\in I$ then either $f\in I$ or $g\in I$ (or both are in $I$). Geometrically this means the associated variety is irreducible (i.e. it cannot be written as a non-trivial union of other varieties) and reduced (meaning a generic point has multiplicity one, e.g. the roots of $x^{2}=0$ have multiplicity two).

Conversely, any prime polynomial ideal $\langle f_{1},\dots,f_{m}\rangle$ in $\mathbb{C}[x_{1},\dots,x_{n}]$, where each $f_{i}$ is a binomial, defines a toric variety ${X}_{A}$ in $\mathbb{C}^{n}$. Chemical reaction networks whose steady state ideal is generated by such prime ideals have been studied in the literature e.g. [19]. To simplify notation we make the following definition.

Definition 2.8**.**

A chemical reaction network whose steady state ideal has an associated prime that is binomial is a toric chemical reaction network.

The fact that toric varieties can be parameterized by monomials over $\mathbb{C}^{*}$ also yields a parametric representation of the positive part of the toric variety, i.e. of the positive steady states of the associated reaction network. For this reason we restrict our analysis to toric networks. We will also wish to explicitly incorporate the reaction rates in our analysis, hence our starting point will be a binomial ideal with coefficients in $\mathbb{R}(\kappa):=\mathbb{R}(\kappa_{1},\dots,\kappa_{m})$, the field of rational functions of the reaction rates with real coefficients. More precisely we consider a prime binomial steady state ideal defined by $\nu$ equations,

[TABLE]

The vectors $b^{+}_{i}$ and $b^{-}_{i}$ are positive column vectors with disjoint support and let $B=(b_{1},\dots,b_{\nu})$ be the matrix with columns $b_{i}=b^{+}_{i}-b^{-}_{i}$, called the exponent matrix of $I$. Then there exists a matrix $A$ such that $b_{i}\in\ker(A),$ for $i=1,\dots,\nu$. The matrix $A$ is a $d\times n$ integer matrix where $d$ is the dimension of the toric variety and $n$ is the dimension of the ambient affine space. Note that, since the dimension of a variety cannot exceed the dimension of its ambient space, $\text{rank}(A)\leq d\leq n$.

Remark 2.9**.**

Note that there is some freedom in the choice of the matrix A, as integer row operations on rows of A will not change its kernel and, hence, the variety it parameterizes is left unaltered. Therefore, we informally refer to any matrix $A$ as the $A$-matrix (associated to a toric variety $X_{A}=V(I)$).

To explicitly incorporate the reaction rates in the monomial parameterization we will now define the scaled toric variety of a matrix $A$ and a vector $x^{*}\in\mathbb{C}^{n}$. The value $x^{*}$ can be thought of a particular steady state which in turn depends on the reaction rates, see Remark 2.13.

Definition 2.10** (Scaled Toric Variety).**

Given an $x^{*}=(x_{1}^{*},\dots,x_{n}^{*})\in\mathbb{C}^{n}$ define a (Laurent) monomial map $\psi_{A}^{(x^{*})}:=\psi_{A}$ with

[TABLE]

for $a_{i}$ a column of $A$. The (Zariski) closure of the image of this monomial map defines a scaled toric variety, $X_{A,x^{*}}=\overline{\psi_{A}((\mathbb{C}^{*})^{d})}$.

The monomial map $\psi_{A}$ also induces a parameterization map.

Definition 2.11**.**

The parameterization map defined by the $A$-matrix is the $\mathbb{C}$-algebra homomorphism

[TABLE]

Note, as above, the map $\phi_{A}$ depends on $x^{*}\in\mathbb{C}^{n}$.

Example 2.12**.**

The steady state ideal of (1) is

[TABLE]

which is a (prime) binomial ideal in the ring $\mathbb{R}(\kappa_{1},\dots,\kappa_{6})[x_{1},\dots,x_{6}]$. Hence, we can find a parameterization $x_{1}=x_{1}^{*}t_{1},\>x_{2}=x_{2}^{*}t_{2},\>x_{3}=x_{3}^{*}t_{1}^{-1}t_{3}^{-1},\>x_{4}=x_{4}^{*}t_{2}^{-1}t_{3}^{-1},\>x_{5}=x_{5}^{*}t_{3}^{-1},\>x_{6}=x_{6}^{*}t_{3}^{-1}$ corresponding to the A-matrix

[TABLE]

Remark 2.13**.**

As shown in [19] the entries of $x^{*}$ are algebraic functions in $\kappa$. In practice we wish to restrict our choices of $\kappa$ to those such that when we evaluate $x^{*}$ at some fixed $(\kappa_{1},\dots,\kappa_{m})\in\mathbb{R}_{>0}^{m}$ we have that $x^{*}\in\mathbb{R}_{>0}^{n}$.

The proposition below connects the number of conserved quantities to the dimension of the toric variety and is straightforward, however, we include a proof as we were unable to locate one in the literature.

Proposition 2.14**.**

Consider a chemical reaction network with steady state ideal $I=\langle\alpha^{T}A_{k}x^{\alpha}\rangle$ in the polynomial ring $\mathbb{R}[x_{1},\dots,x_{n}]$. Let $V=\overline{V(I)\cap(\mathbb{R}^{*})^{n}}$ be the real steady state variety and denote the dimension of $V$ by $\dim(V)$. Let the linear equations $\ell_{1},\dots,\ell_{d}$ be the corresponding conservation relations. Suppose that $\dim({V\cap V(\ell_{1},\dots\ell_{d})})=d^{\prime}$. Then $\dim(V)=d+d^{\prime}$. Further if $V$ is a toric variety we have that $\dim(V\cap(\mathbb{R}_{>0})^{n})=d+d^{\prime}$.

Proof.

Each linear form $\ell_{j}$ contains a constant term $c_{j}\in\mathbb{R}$ which can be chosen freely. For a sufficiently general choice of $c_{1},\dots,c_{d}$ the intersection $V\cap V(\ell_{1},\dots\ell_{d})$ is transverse, it follows that $d^{\prime}=\dim(V\cap V(\ell_{1},\dots\ell_{d}))=\dim(V)-\dim(V(\ell_{1},\dots\ell_{d}))=\dim(V)-d$. Because $V$ is toric it is parameterized by monomials, and hence its dimension in the positive orthant is the same as its dimension over $\mathbb{R}^{n}$.

∎

Remark 2.15**.**

In the notation of Proposition 2.14 any $d^{\prime}>0$ implies that there is an infinite number of steady states for each general choice of $c_{1},\dots,c_{d}$. Hence, the case $d^{\prime}>0$ renders the discussion of multistationarity irrelevant. Therefore for the remainder of this paper we assume that $d^{\prime}=0$ meaning that the dimension of the toric variety is equal to the number of conservation relations.

Example 2.16** (Example 2.12 cont.).**

The network of Example 2.12 it is known to have a finite number of positive steady states for each set of initial conditions ($d^{\prime}=0$) and there are three independent conservation relations $x_{1}+x_{5}=c_{1}$, $x_{2}+x_{6}=c_{2}$ and $x_{3}+x_{4}+x_{5}+x_{6}=c_{3}$. This implies that the dimension of the (toric) steady state variety, $V(I)$, is $\dim(V(I))=3$ which is indeed the case.

2.3 Matroids

One of the focuses of this paper is to find and study independent subsets of chemical species. Independent subsets can give valuable information about which species concentrations have to be measured and which concentrations can be determined from measurements [6]. A simple example of linear independence is the set of three vectors $v_{1}=(1,0)^{T},\>v_{2}=(0,1)^{T},\>v_{3}=(2,1)^{T}$. The set of vectors $\{v_{1},v_{2}\}$ is linearly independent as there exists no $\lambda\in\mathbb{R}$ such that $\lambda v_{1}=v_{2}$, whereas $v_{3}$ can be obtained by $v_{1}+v_{2}$ and $\{v_{1},v_{2},v_{3}\}$ is therefore dependent. The vectors $v_{1}$ and $v_{2}$ are said to from a basis of the set $\{v_{1},\>v_{2},\>v_{3}\}$ which is a familiar notion of linear algebra. The set $\{v_{1},\>v_{2},\>v_{3}\}$ is minimally dependent as it contains only a single element other than a basis and is a circuit. Matroid theory extends the notion of independence to polynomials rings [29]. First, we define a matroid.

Definition 2.17** (Matroid).**

A matroid $\mathcal{M}(E,\mathcal{I})$ is a pair of two finite sets $(E,\mathcal{I})$ where $E$ is the ground set and $\mathcal{I}$ is a set of subsets of $E$, called independent sets, satisfying the following conditions.

The empty set, $\emptyset$, is independent such that $\emptyset\in\mathcal{I}$ and, hence, $\mathcal{I}\neq\emptyset$. 2. 2.

If $i\in\mathcal{I}$ and $i^{\prime}\subseteq i$ then $i^{\prime}\in\mathcal{I}$. This is called the hereditary property. 3. 3.

If $i_{1},\>i_{2}\in\mathcal{I}$ and $|i_{1}|<|i_{2}|$ then there exists an element $x\in i_{2}-i_{1}$ such that $i_{1}\cup x\in\mathcal{I}$. This is the exchange property.

The notions of matroid basis, rank and circuit will also be useful.

Definition 2.18**.**

Matroid bases, rank and circuits are defined as follows:

• •

A basis $S$ of a matroid $\mathcal{M}$ is a maximally independent subset, i.e. a subset $S\in\mathcal{I}$ such that $S\cup k$ is a dependent set for all $k\in E-S$. Define the set of all bases to be $\mathcal{B}$.

• •

The rank $\rho$ is a function on a subset of $E$ which takes a set $e\subseteq E$, and returns the cardinality of the largest subset $i\subseteq e$ which also satisfies $i\in\mathcal{I}$.

• •

The circuits $C$ are minimally dependent sets, i.e. subsets of $E$ which satisfy $C\not\in\mathcal{I}$ and $(C-z)\in\mathcal{I}$ for any $z\in C$.

In particular for a matroid $\mathcal{M}$ it holds that all bases of $\mathcal{M}$ are equicardinal. An important notion in matroid theory is also the rank of a matroid.

Definition 2.19**.**

The rank of a matroid $\mathcal{M}$, $\rho(\mathcal{M})$, is the cardinality of a basis.

Due to the general definition of a matroid many mathematical objects have a matroid structure such as sets of vectors in a vector space or graphs. One class of matroids relevant for chemical reaction networks are algebraic matroids. Algebraic matroids encode the algebraic dependencies between the variables of a polynomial ring $\mathbb{C}[x_{1},\dots,x_{n}]$ in a prime ideal $P$. A more detailed discussion of these ideas in given in Section 4. For information on how to compute algebraic matroids we refer the reader to [29, 39]. Since we are interested in the relation between the matroids of families of reaction networks we introduce the notion of a submatroid.

Definition 2.20** (Submatroid).**

Given a subset $E^{\prime}\subset E$, a matroid $\mathcal{M}^{\prime}(E^{\prime},\mathcal{I}^{\prime})$ is a submatroid of $\mathcal{M}(E,\mathcal{I})$ if $\mathcal{I}^{\prime}=\mathcal{I}\cap\mathcal{P}(E^{\prime})$, where $\mathcal{P}(E^{\prime})$ is the power set of $E^{\prime}$. The rank function of $\mathcal{M}^{\prime}$ is the restriction $\rho_{\mathcal{M}^{\prime}}(S)=\rho_{\mathcal{M}}(S)$ for $S\subseteq E^{\prime}$.

Submatroids of algebraic matroids contain a subset of the polynomial relations between the variables of circuits of the original matroid. The relations between the variables of circuits are referred to as circuit polynomials. Below we will consider algebraic matroids of toric networks and the relations of matroids within a family of networks.

Example 2.21**.**

The algebraic matroid of 2.12 can be calculated using the computer algebra system Macaulay2 [31]. It is a matroid on the ground set $E=\{x_{1},\dots,x_{6}\}$ of rank $3$ with $12$ bases. Hence, as we will show in later sections, measuring only four species will suffice to reject the one site model, any three species in a basis and one extra species. Further, the matroid of the two site distributive system is a rank $3$ matroid with $61$ bases and the matroid of the two site processive network is a rank $3$ matroid with $19$ bases. Both two-site model matroids contain the one site matroid as a submatroid.

3 Families of Reaction Networks

In this paper we obtain results which hold for a range of networks rather than a single network. When the properties studied (i.e. multistationarity) are present in a “small network”, these properties lift to all larger networks constructed by the procedure outlined below. We call the collection of all such networks a family $\mathcal{N}$ with members $\mathcal{N}_{i}$. In this section we rigorously define families of networks and give examples of chemical systems fulfilling the family property. Furthermore, we prove some results on families which have toric steady states.

3.1 Family Definition and Properties

Families of networks can be found in a wide range of biochemical settings such as multi-site phosphorylation [13] (e.g. cellular signalling, DNA transcription, cell death), kinetic proofreading [1] (immune response) or compartmentalised diffusion [40] (spatial models).

We identify a family by properties of its reaction graphs. If we can construct the graph of another network from a given network by a fixed set of procedures, the networks are in the same family; an example of such a construction is given in Figure 1.

Remark 3.1** (Informal Construction of Families of Graphs).**

Fix a labelled digraph $\mathcal{G}_{n}$ and an unlabelled digraph $M$; from these build new family members step-by-step. At each step assign a new set of labels to $M$, rendering $M$ a labelled digraph. The vertices of $\mathcal{G}_{n+1}$ are the union of those of $\mathcal{G}_{n}$ and of $M$. The edges of the graph $\mathcal{G}_{n+1}$ are obtained from $\mathcal{G}_{n}$ as follows:

Add a set $\mathcal{E}$ of edges adjacent to both some vertices of $M$ and some vertices of $\mathcal{G}_{n}$; the resulting graph $\mathcal{G}_{n}^{\rm int}$ is referred to as the intermediate graph (see also Definition 3.8). 2. 2.

Delete some subset of the edges of $\mathcal{G}_{n}$ from the graph $\mathcal{G}_{n}^{\rm int}$ to form $\mathcal{G}_{n+1}$.

We formalise this construction using a definition from graph theory [41]. Graphs constructed as in Remark 3.1 form a family if they satisfy the conditions of Definition 3.2 below.

Definition 3.2** (See also Section 2 of [41]).**

For a reaction graph $\mathcal{G}=(\mathcal{C},\mathcal{R})$ and for $U\subseteq\mathcal{C}$ let $\mathcal{G}[U]$ be the subgraph induced by $U$ and $N_{\mathcal{G}}(U)$ the set of vertices adjacent to some vertex in $U$. Consider a set of graphs $\{\mathcal{G}_{i}\}_{i\geq 0}$ where $\mathcal{G}_{i}=(\mathcal{C}_{i},\mathcal{R}_{i})$. Set $\mathcal{W}_{0}=\mathcal{C}_{0}$ and $\mathcal{E}_{0}=\mathcal{R}_{0}$; additionally let $\mathcal{W}_{i+1}=\mathcal{C}_{i+1}-\mathcal{C}_{i},$ and $\mathcal{E}_{i+1}=\mathcal{R}_{i+1}-\mathcal{R}_{i}$ for $i\geq 0$. Further, there exists a positive integer $n$, a labelled graph $M$ and a set of reactions $Y$, such that the set $\{\mathcal{G}_{i}\}_{i\geq 0}$ is called a family of graphs if the following properties hold:

$N_{\mathcal{G}_{n}}(\mathcal{W}_{n})\subseteq\mathcal{W}_{0}\cup\left(\bigcup_{i=0}^{r}\mathcal{W}_{n-i}\right)$* for $n>r$,* 2. 2.

$\mathcal{R}_{n}=\left(\mathcal{R}_{n-1}-Y\right)\cup\mathcal{E}_{n}$, where $Y\subseteq\bigcup_{i=1}^{r}\mathcal{E}_{n-i}$, 3. 3.

the graph $\mathcal{G}_{n}\left[\mathcal{W}_{0}\cup\left(\bigcup_{i=0}^{r}\mathcal{W}_{n-i}\right)\right]$ is equal to $M$ for $n>r$. In addition, $\mathcal{G}_{n}[\mathcal{W}_{n}]$ is always the same graph.

Put simply the last condition says that the graph “added” to the previous graph must be the same for each step throughout the family; this is illustrated in Figure 2 and further elucidated in [41].

In the context of chemical reaction networks we add another condition, namely, that in every subsequent family member there appears at least one new chemical species:

1. 4.

$\mathcal{S}_{M}\subsetneq\mathcal{S}_{N}$ for $M<N$.

We assume that in a network with a given set of chemical species every complex which is chemically possible and that can be formed with this set of species is already present in the network; i. e. if in a network with species $H_{2}$ and $O_{2}$ chemistry dictates that the only possible reaction is $2\,H_{2}+O_{2}\to 2\,H_{2}O$ then only the complexes $2\,H_{2}+O_{2}$ and $2\,H_{2}O$ exist for this species set. Therefore, to give rise to a larger network new species have to be added.

Remark 3.3**.**

Note that not every infinite sequence of networks is a family as defined in Definition 3.2. For a non-example see Example 3.7.

Remark 3.4**.**

For ease of notation, the unlabelled graph $M$ is drawn in a “dot and arrow” representation; e.g. the graph $M=\{\bm{\cdot},\bm{\cdot}\to\bm{\cdot}\}$ corresponds to a graph with two connected components, an isolated vertex and a component with two unlabelled vertices and a directed edge between them.

Example 3.5** (Distributive Phosphorylation).**

Distributive phosphorylation with $N=1$ phosphorylation sites follows the reaction scheme

[TABLE]

The reaction scheme is a digraph $\mathcal{G}_{1}=(\mathcal{W}_{1},\mathcal{E}_{1})$ with $\mathcal{W}_{1}=\{E+S_{0},ES_{0},E+S_{1},F+S_{1},FS_{1},F+S_{0}\}$ and $\mathcal{E}_{1}=\{(E+S_{0},ES_{0}),(ES_{0},E+S_{0}),(ES_{0},E+S_{1}),(F+S_{1},FS_{1}),(FS_{1},F+S_{1}),(FS_{1},F+S_{0})\}$. To construct $\mathcal{G}_{2}$ we use \color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\mathcal{W}_{2}=\{ES_{1},E+S_{2},F+S_{2},FS_{2}\} and \color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\mathcal{E}_{2}=\{(E+S_{1},ES_{1}),(ES_{1},E+S_{1}),(ES_{1},E+S_{2}),(F+S_{2},FS_{2}),(FS_{2},F+S_{2}),(FS_{2},F+S_{1})\}. Defining $\mathcal{G}_{2}=(\mathcal{W}_{1}\cup\mathcal{W}_{2},\mathcal{E}_{1}\cup\mathcal{E}_{2})$ gives

[TABLE]

Hence (in the notation of Definition 3.2) , $M=\{\bm{\cdot}\to\bm{\cdot},\>\bm{\cdot}\rightleftharpoons\bm{\cdot}\}$ and $r=1$. The digraphs $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ define the family members $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$. Inductively, this procedure can be continued to member $\mathcal{N}_{N}$. By inspection of the green subgraph we can see that condition 3 of Definition 3.2 is fulfilled for any $N>0$ and, therefore, the distributive phosphorylation networks form a family. Further, $\mathcal{N}_{1}$ is a subnetwork of $\mathcal{N}_{2}$ as defined in [25, Definition 2.2].

Example 3.6** (Processive Phosphorylation).**

Processive phosphorylation of a substrate with $N=1$ phosphorylation sites follows the reaction scheme

[TABLE]

For $N=1$ this is the same reaction scheme as for distributive phosphorylation, except for the relabelling of the fully phosphorylated substrate as product $P$. However, to construct the next family member from the digraph $\mathcal{G}_{1}$ we use \color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\mathcal{W}_{2}=\{ES_{1},FS_{2}\} and \color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}\pgfsys@color@cmyk@stroke{0.64}{0}{0.95}{0.40}\pgfsys@color@cmyk@fill{0.64}{0}{0.95}{0.40}\mathcal{E}_{2}=\{(ES_{0},ES_{1}),(ES_{1},E+P),(F+P,FS_{2}),(FS_{2},F+P),(FS_{2},FS_{1})\}. Next, delete edges $Y_{1}=\{(ES_{0},E+P),(F+P,FS_{1}),(FS_{1},F+P)\}$ which results in the graph $\mathcal{G}_{2}=\left(\mathcal{W}_{1}\cup\mathcal{W}_{2},\left(\mathcal{E}_{1}\cup\mathcal{E}_{2}\right)-Y_{1}\right)$ to give

[TABLE]

Where the vertices $M=\{\bm{\cdot},\bm{\cdot}\}$ were added and $r=1$ as in Definition 3.2. Again, the digraphs $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ define the family members $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$. Condition 3 of Definition 3.2 is fulfilled for any $N>1$ and, therefore, the processive phosphorylation networks form a family.

Example 3.7** (Non-example).**

As mentioned in Remark 3.3 not every infinite sequence of graphs forms a family. Consider the autocatalytic networks

[TABLE]

In the reaction above $\emptyset$ represents production and degradation of a molecule by a mechanism not further studied. Without loss of generality consider $\mathcal{N}_{7}$ and $\mathcal{N}_{8}$; we see that $\mathcal{C}_{8}$ is not a union of $\mathcal{C}_{7}$ and another graph as $X_{7}+X_{1}\in\mathcal{C}_{7}\not\subseteq\mathcal{C}_{8}$.

We conclude this subsection by considering the polynomial equations arising from the reaction graphs of successive members of a family. First, we see from Section 2.1 that every reaction has a unique rate constant, $\kappa_{i}$, which is the edge weight of the reaction digraph. Further, in the dynamical system given by the equations (6) every monomial, which represents a reactant complex, is multiplied by a constant $\kappa_{i}$. Note that isolated vertices in the reaction graph do not contribute to the dynamics of the network. Hence, for families of networks which satisfy Definition 3.2 with $Y=\emptyset$, we can derive the equations of the $N^{\text{th}}$ member of a family from the $(N+1)^{\text{th}}$ member by setting all the edge weights of the edges added to the reaction graph to zero. Formally, we define the evaluation map

[TABLE]

For networks which require the deletion of a set of edges, $Y$, we need the concept of an intermediate network.

Definition 3.8** (Intermediate Network).**

The intermediate network $\mathcal{G}_{n}^{\rm int}$ is the network constructed by only adding reactions to join the newly labelled graph $M$ to $\mathcal{G}_{n}$ and before any edges are deleted from the reaction graph. See step 1 of Remark 3.1.

Remark 3.9**.**

The next member of a family is a subnetwork of the intermediate network. Hence, for families in which the next family member coincides with the intermediate network (i.e. when $Y=\emptyset$ in Definition 3.2) this terminology is not required.

The dynamical equations of the $N^{\text{th}}$ and $(N+1)^{\text{th}}$ member can be constructed from their intermediate network by defining the appropriate evaluation maps. As can be seen from condition 2 of Definition 3.2 the reaction (edge) set of the $(N+1)^{\text{th}}$ member of a family is given by $\mathcal{R}_{N+1}=(\mathcal{R}_{N}-Y)\cup\mathcal{E}_{N}$. The edge sets for the $N^{\text{th}}$ and intermediate network are given by $\mathcal{R}_{N}$ and $\mathcal{R}_{N}\cup\mathcal{E}_{N}$ respectively. Define a function $\mu$ which associates a unique edge weight, $\kappa_{i,j}$, to every edge of the reaction graph,

[TABLE]

This gives the rate constants $\mu\left(\mathcal{R}_{N}\right)=\{\kappa_{1},\dots,\kappa_{m}\}$, $\mu\left(\mathcal{E}_{N}\right)=\{\kappa_{m+1},\dots,\kappa_{m+m^{\prime}}\}$ and $\mu\left(Y\right)=\{\kappa_{1},\dots,\kappa_{m^{\prime\prime}}\}$ (where we relabel rate constants produced by (10) as required). By construction, an edge is not present in the reaction graph if it has an edge weight of [math]. Hence, two evaluation maps can be defined to map the edge sets of the intermediate network to the $(N+1)^{\text{th}}$ and $N^{\text{th}}$ network respectively,

[TABLE]

This results in the edge weights $\pi^{+}\left(\mu\left(\mathcal{R}_{N}\cup\mathcal{E}_{N}\right)\right)=\{\kappa_{m^{\prime\prime}+1},\dots,\kappa_{m+m^{\prime}}\}=\mu\left((\mathcal{R}_{N}-Y)\cup\mathcal{E}_{N}\right)$ and $\pi^{-}\left(\mu\left(\mathcal{R}_{N}\cup\mathcal{E}_{N}\right)\right)=\{\kappa_{1},\dots,\kappa_{m}\}=\mu\left(\mathcal{R}_{N}\right)$, which, after taking the inverse of $\mu$, can be associated to the reaction sets of the $(N+1)^{\text{th}}$ and $N^{\text{th}}$ network respectively. In particular, if $\mathcal{G}_{N}^{\rm int}$ is the intermediate graph associated to $\mathcal{G}_{N}$ then applying $\pi^{-}$ to (the edge set of) $\mathcal{G}_{N}^{\rm int}$ will yield $\mathcal{G}_{N}$ and applying $\pi^{+}$ to (the edge set of) $\mathcal{G}_{N}^{\rm int}$ will yield $\mathcal{G}_{N+1}$.

Example 3.10** (Processive Phosphorylation).**

[TABLE]

The reaction scheme for the first intermediate network of the processive phosphorylation family. If $\kappa_{1}=\kappa_{2}=\kappa_{3}=0$ then the $N=2$ site network is obtained. Equivalently, if we let $\kappa_{7}=\kappa_{8}=\kappa_{9}=\kappa_{10}=\kappa_{11}=0$, then we obtain the $N=1$ site network. The corresponding steady state ideals can be obtained in this manner from the ideal from the intermediate network.

3.2 Toric Families

As described in Definition 2.8, a reaction network is called toric if it has toric steady states. Showing whether a chemical reaction network is toric is a non-trivial task; previous results have used deficiency theory **[11]**, network structure based methods **[19]** or methods based on detecting binomiality **[44]**. Formally, one needs to find the associated primes of a steady state ideal and compute a Gröbner basis of each. If the reduced Gröbner basis is binomial, then the ideal is a prime binomial ideal **[43]**. However, as every toric variety has non-zero components, removing any components of the steady state variety contained in the coordinate axes is often a first step towards identifying the toric components. Algebraically this removal is accomplished by computing the saturation $I^{\infty}=I:(x_{1}\cdots x_{n})^{\infty}$, **[45]**. In this subsection we prove a series of small results regarding the steady state ideals of families of toric networks. First, we define a toric family.

Definition 3.11** (Toric family).**

If every member of a family of networks is a toric chemical reaction network, the family is called a toric family.

Remark 3.12**.**

Prominent examples of toric families are multisite phosphorylation networks or compartmentalised diffusion networks.

Lemma 3.13 states that saturation and the evaluation map, $\pi$, commute. Therefore, the network operation of deleting edges in the reaction graph can be carried out before, or after, finding the prime binomial ideal in the steady state ideal.

Lemma 3.13**.**

Let $\pi^{-}$ be the evaluation map (12) and let $I_{N+1},\;I_{N}\subseteq\mathbb{R}[\kappa_{1},\dots,\kappa_{m+m^{\prime}},x_{1},\dots,x_{n+n^{\prime}}]$ be the steady state ideals of the $(N+1)^{\text{th}}$ and $N^{\text{th}}$ family member respectively. Then it holds that

1. (i)

$\pi^{-}(I_{N+1}):(x_{1},\dots,x_{n+n^{\prime}})^{\infty}=\pi^{-}(I_{N+1}:(x_{1},\dots,x_{n+n^{\prime}})^{\infty})$, 2. (ii)

$\pi^{-}(I_{N+1}:(x_{1},\dots,x_{n+n^{\prime}})^{\infty})=I_{N}:(x_{1},\dots,x_{n})^{\infty}$**

The proof of this lemma uses the the geometric interpretation of the statements of the lemma and can be found in Appendix A. Suppose that $\mathcal{N}$ is a toric family of reaction networks with the $(N+1)^{th}$ member of the family having toric steady states specified by the binomial ideal $I_{N+1}$. Lemma 3.13 establishes a useful connection between species (corresponding to variables $x$), and reactions (corresponding to $\kappa$), or, equivalently, between vertices and edges of a reaction graph. By the extra Condition 4 of Definition 3.2 every new edge must either originate or end on a complex containing a new species. Hence, by applying the evaluation map $\pi^{-}$ an ideal of an intermediate network $I^{\text{int}}_{N}$ ($=I_{N+1}$* for nested subnetworks) is automatically mapped from $\mathbb{R}[x_{1},\dots,x_{n+n^{\prime}}]$ (corresponding to a network with $n+n^{\prime}$ chemical species) to an ideal in the ring $\mathbb{R}[x_{1},\dots,x_{n}]$ (corresponding to a network with $n$ species).*

Example 3.14** (Processive Phosphorylation (cont.)).**

The steady state ideal corresponding the network of Example 3.10 is

[TABLE]

As can be seen using the evaluation map $\pi^{-}$, which sets $\kappa_{7}=\cdots=\kappa_{11}=0$, the variables $x_{7}$ and $x_{8}$ disappear from the ideal. Further, computing the saturation $I:(x_{1}\cdots x_{8})^{\infty}$ and the appropriate evaluation maps we can show that the relations of Lemma 3.13 hold using Macaulay2 [31].

Proposition 3.15 below shows that the image of a binomial ideal under the evaluation map $\pi$ is either the constant ideal (corresponding to an empty variety) or another binomial ideal. The proof can be found in Appendix A.

Proposition 3.15**.**

Let $R=\mathbb{R}(\kappa_{1},\dots,\kappa_{m+m^{\prime}})[x_{1},\dots,x_{n+n^{\prime}}]$ and define the partial evaluation map $\pi:R\to\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x_{1},\dots,x_{n+n^{\prime}}]$. Suppose that $I_{N+1}:(x_{1}\cdots x_{n+n^{\prime}})^{\infty}\subset R$ is a binomial ideal. Then $I_{N}=\pi(I_{N+1}):(x_{1}\cdots x_{n})^{\infty}\subset\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x_{1},\dots,x_{n}]$ is either a binomial ideal or the ideal $(1)$. Further, for a fixed choice of rate constants $I_{N+1}:(x_{1}\cdots x_{n+n^{\prime}})^{\infty}\subset\mathbb{R}[x_{1},\dots,x_{n+n^{\prime}}]$. Then if $V(I_{N+1}:(x_{1}\cdots x_{n+n^{\prime}})^{\infty})\subset\mathbb{R}^{n}$ is a toric variety and if $V(I_{N}:(x_{1}\cdots x_{n})^{\infty})\cap\mathbb{R}_{>0}^{n}\neq\emptyset$ then $\overline{V(I_{N})\cap\mathbb{R}_{>0}^{n}}$ is a toric variety.

Hence, by Proposition 3.15 if the intermediate network of the $(N+1)^{\text{th}}$ and the $N^{\text{th}}$ member of a family is toric then both, the $N^{\text{th}}$ and $(N+1)^{\text{th}}$ members are either toric or have empty positive steady state varieties.

Example 3.16** (Example 3.14 cont.).**

Explicit computations is Macaulay2 show that the ideal $I:(x_{1}\cdots x_{8})^{\infty}$ is toric. We also know that the processive phosphorylation network has positive steady states. Hence, the positive steady states of the $N=1$ and $N=2$ networks are cut out by toric varieties.

A relationship between the $A$-matrices of successive members of toric families is now found by considering the exponent matrices of the binomial ideals and their corresponding $A$-matrices.

Theorem 3.17**.**

Let $I_{N+1}:(x_{1}\cdots x_{n+n^{\prime}})^{\infty}\subset\mathbb{R}(\kappa_{1},\dots,\kappa_{m+m^{\prime}})[x_{1},\dots,x_{n+n^{\prime}}]$ be a binomial ideal defining the positive steady states of the $(N+1)^{\text{th}}$ member of a family of reaction networks. Also let $\pi$ be the evaluation map which sends $\kappa_{m+1}=\cdots=\kappa_{m+m^{\prime}}=0$ and assume that $I_{N}=\pi(I_{N+1})\subset\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x_{1},\dots,x_{n}]$. Then $I_{N}:(x_{1}\cdots x_{n})^{\infty}$ is binomial. Further, for a choice of generators, let $B_{N+1}$, $B_{N}$ be the matrices associated to the exponents of the binomial ideals $I_{N+1}$ and $I_{N}$. Similarly let $A_{N+1}$ and $A_{N}$ be the $A$-matrices corresponding to $B_{N+1}$ and $B_{N}$. Then:

1. (i)

$B_{N}$* is a submatrix of $B_{N+1}$ and,* 2. (ii)

$A_{N}$* is a submatrix of $A_{N+1}$, if the submatrix of $B_{N+1}$ formed by the columns which have at least one non-zero entry in the last $n^{\prime}$ rows has rank at most $n^{\prime}$.*

The proof can be found in Appendix A. An interesting Corollary of the above Theorem occurs if the number of conservation relations is constant.

Corollary 3.18**.**

Suppose the assumption of Remark 2.15 holds and that the number of conservation relations remains constant. Then, the matrix $A_{N}$ is obtained from $A_{N+1}$ by deleting the last $n^{\prime}$ of its columns.

Proof.

When the dimensions of the varieties of the networks $\mathcal{N}_{N+1}$ and $\mathcal{N}_{N}$ are the same, then so do the dimensions of the left kernels of their $B$-matrices. Hence, by Theorem 3.17 this is the case if and only if $\text{rank}(\tilde{B})=n^{\prime}$ and the result follows by the Rouché-Capelli theorem. ∎

This section is concluded by relating the mathematical insights of the above theorems to families of toric chemical reaction networks.

Theorem 3.19**.**

The $A$-matrices of members of a family of chemical reaction networks are submatrices of each other when the assumption of Remark 2.15 is valid, all networks considered have the same number of conservation relations and one of the following conditions hold:

1. (i)

The family consists of a sequence of subnetworks with toric steady states. 2. (ii)

The family is toric with toric intermediates.

The proof can be found in Appendix A.

Example 3.20** (Processive Phosphorylation (cont.)).**

The processive network is toric with toric intermediates, every family member has a finite number of positive steady states for each set of initial conditions and the number of conservation relations is constant. The $A$-matrix of the intermediate network of Example 3.10 is given by

[TABLE]

which is identical to the $A$-matrix of the $N=2$ network. The $A$-matrix of the $N=1$ member is obtained by deleting the last two columns (marked in red) of the $A$-matrix associated to the intermediate network.

4 Matroid Theory for Toric Chemical Reaction Networks

In this section we study biochemical and algebraic question 2 regarding parameter estimation and model rejection. We start by showing the equivalence of two matroids, an algebraic matroid and a much simpler linear matroid. We then decorate the linear matroid with binomials and finally apply these results to model rejection. First, we prove the that following is a matroid.

Proposition 4.1**.**

Let $E=\{x_{1},\dots,x_{n}\}$, $S=\{x_{i_{1}},\dots,x_{i_{j}}\}\subseteq E$ and define $\phi_{A}(S)$ by $\phi_{A}(S)=\{\phi_{A}(x_{i_{1}}),\dots,\phi_{A}(x_{i_{j}})\}$, where $\phi_{A}$ is the parameterization map of Definition 2.11. The ground set $E$ with the set

[TABLE]

is a matroid $\mathcal{M}(E,\mathcal{I})$. Further, this matroid is isomorphic to the matroid defined by the column vectors of the $A$-matrix.

Proof.

Consider the image of $E$ under $\phi_{A}$; this is a set of monomials. A set $S$ of monomials is algebraically independent if and only if there exists no polynomial $p\in k[t^{\pm}_{1},\dots,t^{\pm}_{d}]$ such that $p(S)=0$. This is exactly the condition to give an independent set $i\in\mathcal{I}$. Further, by [48, Lemma 4.2.10], a set of monomials is independent if and only if their exponent vectors are linearly independent. Hence, algebraic independence of $\phi_{A}(S)$ is equivalent to linear independence of the columns of the matrix $A$ corresponding to $S$. ∎

Definition 4.2**.**

The matroid $\mathcal{M}(A)$ defined by the column vectors, $a_{i}$, of a full rank integer matrix $A$ encodes the algebraic dependencies of the concentrations at a non-zero (usually positive) steady state of a chemical reaction network and is therefore is called the positive steady state (PSS) matroid.

The next theorem shows that the PSS matroid and the algebraic matroid defined by the binomial ideal are identical. Therefore, the algebraic matroid can be studied directly using linear algebra operations on the columns of the $A$-matrix defining the PSS matroid. This way bases, circuits and even circuit polynomials can be inferred trivially.

Theorem 4.3**.**

Let $I_{b}\subseteq R=\mathbb{C}[x_{1},\dots,x_{n}]$ be a prime binomial ideal defining a toric variety $X_{A,x^{*}}=V(I_{b})$ where $A\in\mathbb{Z}^{d\times n}$ is the exponent matrix of the corresponding monomial parameterization. Denote the algebraic matroid defined by $I_{b}$ as $\mathcal{M}(I_{b})$ and the matroid defined by the linear independence of the columns of $A$ by $\mathcal{M}(A)$. Then $\mathcal{M}(I_{b})=\mathcal{M}(A)$.

The proof can be found in Appendix A. We now proceed to find decorations of the PSS matroid’s circuits which reveal the full power of the PPS matroid formulation.

Definition 4.4** (Laurent Binomial Associated to a PSS matroid circuit).**

Let $X_{A,x^{*}}$ be the toric variety defined by a full rank matrix $A$ and a positive vector $x^{*}$ as in Definition 2.10. Let $\mathcal{M}(A)$ be the associated PSS matroid (Definition 4.2). Consider a circuit $C=\{a_{i_{1}},\dots,a_{i_{j-1}}\}\cup\{a_{i_{j}}\}\subseteq E$ of the matroid $\mathcal{M}(A)$, where $\{a_{i_{1}},\dots,a_{i_{j-1}}\}\in\mathcal{I}$. We define the Laurent polynomial

[TABLE]

with $\lambda_{i_{l}}\in\mathbb{Z}$ chosen such that $\sum_{l=1}^{j-1}\lambda_{i_{l}}a_{i_{l}}={\lambda_{i_{j}}}a_{i_{j}}$ and ${\lambda_{i_{j}}}$ is positive (this is possible since $C$ is a circuit). The expression $\Phi(C)$ is called the Laurent binomial associated to $C$.

Definition 4.5** (Binomial Associated to a PSS matroid circuit).**

Let $\Phi(C)$ be a Laurent binomial of a circuit $C$ of a PSS matroid as in Definition 4.4. The binomial associated to $C$ is $\overline{\Phi(C)}$ which is $\Phi(C)$ with the denominator cleared, i.e.

[TABLE]

where $\lambda^{+}_{j}=\lambda_{j}$ if $\lambda_{j}>0$ and zero otherwise and where $\lambda^{-}_{j}=|\lambda_{j}|$ if $\lambda_{j}<0$ and zero otherwise.

Example 4.6**.**

For the $N=1$ site processive phosphorylation network the columns $C=\{a_{1},a_{2},a_{3}\}\cup\{a_{4}\}$ form a circuit with a linear relation $a_{1}-a_{2}+a_{3}=a_{4}$. The Laurent binomial associated to this circuit is

[TABLE]

and clearing the denominator gives the binomial

[TABLE]

By the following lemma the (Laurent) binomials associated to a PSS matroid are zero on the steady state variety. Hence, they form a model invariant.

Lemma 4.7**.**

Let $X_{A,x^{*}}$ be the toric variety defined by a full rank matrix $A$ and a positive vector $x^{*}$ as in Definition 2.10 and let $\mathcal{M}(A)$ be the associated PSS matroid. If $C$ is a circuit in $\mathcal{M}(A)$ then $\Phi(C)(x)=0$ if $x\in(\mathbb{C}^{*})^{n}\cap X_{A,x^{*}}$.

Proof.

Given a linearly dependent set of vectors it holds that, without loss of generality, $\sum_{l=1}^{j-1}\lambda_{i_{l}}a_{i_{l}}={\lambda_{i_{j}}}a_{i_{j}}$ for some integers $\lambda_{i_{l}}\in\mathbb{Z}$ and ${\lambda_{i_{j}}}\in\mathbb{Z}_{>0}$. It follows that

[TABLE]

Taking the preimage under $\phi_{A}$ this can be rewritten as

[TABLE]

∎

The following theorem states that a set of suitably chosen binomials contains only all positive steady states as our original ideal.

Theorem 4.8**.**

Let $\mathcal{M}(A)$ be the matroid associated to a toric chemical reaction network $\mathfrak{N}$ and choose a basis $S$ and $n-d$ circuits $C_{i}$ such that $\bigcap_{i}C_{i}=S$. Then, the following holds

[TABLE]

Hence, proving that the variety $V(\langle\overline{\Phi(C_{1})},\dots,\overline{\Phi(C_{n-d})}\rangle)$, when intersected with the subspace spanned by the conservation relations, contains at least two positive points is necessary for proving multistationarity of the original toric network given by $I_{\mathfrak{N}}$.

Proof.

Let $F_{k}(x_{1},\dots,x_{n})=\overline{\Phi(C_{k})}$ to give $W=V(F_{1},\dots,F_{n-d})$ and also let

[TABLE]

The containment $X_{A,x^{*}}\cap\mathbb{R}_{>0}^{n}\subseteq W\cap\mathbb{R}_{>0}^{n}$ follows from Lemma 4.7. Now prove the other containment. For a positive real point $w\in W\cap\mathbb{R}_{>0}^{n}$ we must have for each $k=1,\dots,n-d$ that $F_{k}(w)=0$. Hence, by Definition 4.5,

[TABLE]

For each circuit define the vector $\tilde{\lambda}\in\mathbb{Z}^{n}$ as the vector with non-zero entries $i_{\ell}$ corresponding to the values of the coefficients of $\lambda_{i_{j}}a_{i_{j}}-\sum_{y=1}^{j-1}\lambda_{i_{y}}a_{i_{y}}=0$. and let $\Lambda$ be the matrix with rows $\tilde{\lambda}_{k}$; note that the rows of $\Lambda$ form a basis of $\ker(A)$. Set $\tilde{w}=(w_{i_{1}},\dots,w_{i_{j}})$ and set $\tilde{x^{*}}=(x^{*}_{i_{1}},\dots,x^{*}_{i_{j}})$; then

[TABLE]

Further, $\log(\tilde{w}/\tilde{x^{*}})\in\ker(\Lambda)$ and, hence, $\tilde{w}/\tilde{x^{*}}$ is in the image of $x^{*}t^{A}$. ∎

As a result of Theorem 4.8 we can restrict the study of positive steady states to the study of the PSS matroid. The next proposition establishes a connection between the PSS matroids of all members of a family.

Proposition 4.9**.**

Fix a family of toric reaction networks $\mathcal{N}$ and let $\mathcal{N}_{M},\,\mathcal{N}_{N}\in\mathcal{N}$ with $M<N$. If both members of the family have the same number of conservation relations then $\mathcal{M}(A_{M})$ is a submatroid of $\mathcal{M}(A_{N})$.

Proof.

If the dimensions of the varieties are the same, then by Theorem 4.3 the $A$-matrix of $\mathcal{N}_{M}$ is a submatrix of the $A$-matrix of $\mathcal{N}_{N}$. Therefore, $E_{M}\subset E_{N}$ and, trivially, any set of linearly independent columns of any matrix is also linearly independent set of columns in any submatrix containing these. Hence, $\mathcal{I}_{M}=\mathcal{I}_{N}\cap\mathcal{P}(E_{M})$. ∎

4.1 Experimental Design and Compatibility

We now study how steady state matroids and submatroids can be employed in experimental design and model rejection. For the remainder of this subsection we assume that the number of conservation relations within a family is fixed. Hence, by Proposition 4.9 the matroids of smaller family members are submatroids of the matroids of larger family members.

Previous related work includes the study of “complex-linear steady state invariants” **[32]** and data coplanarity **[33]**. A study using the language of algebraic matroids explicitly can be found in **[6]**. In this section we obtain results similar to **[33, 34]** using the PSS matroid and the binomials of Definition 4.5 and show how they can be used in model rejection and experimental design for entire toric families. In particular, we focus on two extreme cases, the one of all parameters being known and the other of perfect, noise-free measurements. Statistical version of our results could form part of a future work.

First, we determine which species need to be measured to be able to construct the steady state locus for an entire family (if the rate constants are known).

Proposition 4.10**.**

Fix a family of chemical reaction networks $\mathcal{N}$ for which the rate constants $\kappa$ associated to every family member are known. There exists a subset of species of the smallest family member which is sufficient to measure at steady state in order to construct a corresponding positive steady state for every subsequent family member.

Proof.

Recall that $x^{*}$ is a positive vector for known constants. Choose a basis $S$ of the PSS matroid of the smallest family member $\mathcal{N}_{1}$. Hence, by Definition 4.5 binomial relations can be constructed to determine the steady state concentrations of the chemical species not in the basis. By Proposition 4.9 any basis of $\mathcal{N}_{1}$ is a basis of the subsequent family members and, hence, Definition 4.5 applies. ∎

Hence, by Proposition 4.10, measuring the concentrations of all species in a basis of the smallest network in a family is sufficient to determine the expected values of the chemical concentrations at steady state of the entire family.

Example 4.11**.**

A detailed study of positive steady state parameterizations for the processive system has been carried out in [13]. Suppose we know all $11$ rate constants of the one and two site processive phosphorylation networks of the running example. In this case the functions $x^{*}_{i}(\kappa)$ can be found and evaluated explicitly (e.g. see [13]). Further, the columns $a_{1},a_{2},a_{3}$ form a basis of $\mathbb{R}^{3}$ and hence the column space of the $A$ matrix. Therefore, given measurements of $x_{1},x_{2},x_{3}$, the expected steady state concentrations of any species $x_{j}$ with $j>3$ can be found as a zero of the Laurent binomial

[TABLE]

Next, we use the PSS matroids for model selection or model rejection for perfect, that is noise-free, data. However, the result of Lemma 4.12 could also be applied to noisy data by following the construction in **[34]**. To determine whether a model is compatible with observed data it is necessary to determine whether there exists a set of parameters $\{\kappa_{1},\dots,\kappa_{m}\}>0$ such that a measured data point $\{\xi_{i_{1}},\dots,\xi_{i_{j}}\}$ is an element of the steady state variety. We proceed by formulating a parameter-free condition for model compatibility.

Lemma 4.12**.**

Let $\mathfrak{N}$ be a reaction network with PSS matroid $\mathcal{M}(A)$. Assume that the rate constants are fixed but unknown, but the initial conditions are varied. Fix a circuit $C$ corresponding to the linear relations among the columns of $A$ via $\sum_{l=1}^{j-1}\lambda_{i_{l}}a_{i_{l}}=\lambda_{i_{j}}a_{i_{j}}$. Given a pair of steady state measurements $\{\xi_{i_{1}},\dots,\xi_{i_{j}}\}$ and $\{\zeta_{i_{1}},\dots,\zeta_{i_{j}}\}$ of the concentrations of $C$, the corresponding model is compatible only if

[TABLE]

Proof.

As in Remark 2.13 we fix rate constants $\kappa=(\kappa_{1},\dots,\kappa_{m})^{T}\in\mathbb{R}^{m}_{>0}$ such that $x^{*}\in\mathbb{R}^{n}_{>0}$. Rearrange $\Phi(C)$ of Definition 4.4 to give

[TABLE]

For the measurements to be compatible we must have that when we evaluate the expression above at $x_{i_{l}}=\xi_{i_{l}}$ and at $x_{i_{l}}=\zeta_{i_{l}}$ we obtain the same value, $\theta$. The conclusion follows. ∎

Example 4.13**.**

Suppose we wanted to test whether some data we obtained could originate from a processive phosphorylation system. We set up the experiment using two different initial conditions and measured the concentrations of the elements of the circuit $\{x_{1},x_{2},x_{3},x_{4}\}$ and based on the results of the previous section we know that $(x_{1}x_{3})/(x_{2}x_{4})=\,const.$ Let the two ideal data points be $\xi=\{5/2,3/10,10/8,875/96\}$ and $\zeta=\{1/2,3/5,1/8,35/384\}$. A priori these data seem unrelated, however, by using the circuit information we find that

[TABLE]

Therefore, the data are compatible with the processive model. Note that the relation holds equally for the one-site and two-site model and we can conclude that the entire family is compatible.

Due to the additional structure provided by the PSS matroid the condition of Lemma 4.12 is much simpler than the linear algebra condition of **[33]**. By Proposition 4.9 it is easy to see that if Lemma 4.12 holds for a given PSS matroid it holds for all of its submatroids. Hence, measuring only a subset of species which belongs to a small member of the family tells us that the measurements of this subset of species is compatible for all subsequent family members. This allows us to determine that a given data set is compatible with a family of networks, but we cannot specify which network in the family is ‘most’ compatible with the data.

Identifying model parameters for perfect data has been studied extensively in previous work **[49, 50, 51]** and, therefore, the discussion in this chapter is restricted to the novel methods obtained by using the PSS matroid of a family of networks. In particular, the PSS matroid allows for straightforward conversion of a variety in the space of species to a variety in the space of parameters. We first show that, in the case of toric steady states the biochemically viable parameter sets, $\kappa=(\kappa_{1},\dots,\kappa_{m})^{T}\in\mathbb{R}^{m}_{>0}$, are the positive part of an algebraic variety and then generalise this result to the entire family.

Proposition 4.14**.**

Let $A$ be a full rank $d\times n$ integer matrix and let $C_{1},\dots,C_{n-d}$ be a collection of circuits of the matroid $\mathcal{M}(A)$, each containing the same basis $S$. Using Definition 4.4 to obtain the ideal $J=\langle\overline{\Phi(C_{1})},\dots,\overline{\Phi(C_{n-d})}\rangle\subseteq R=\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x_{1},\dots,x_{n}]$ and denoting the variables present in a circuit as $x(C_{1})\subseteq\{x_{1},\dots,x_{n}\}$, then the intersection ideal $J_{C_{i}}=J\cap\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x(C_{i})]$ is principal (generated by one element) with generator $\overline{\Phi(C_{i})}$.

Proof.

By construction both, the numerator and the denominator of $\Phi(C_{i})$ contain only variables in $C_{i}$ and, hence, after clearing the denominator the resulting polynomial $\overline{\Phi(C_{i})}$ also only contains variables in $C_{i}$. Since $C_{i}$ is a circuit, the ideal $J_{C_{i}}$ has codimension 1 in $R\cap\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x(C_{i})]$ and, hence, by [52, I.,§7, Proposition 4] it is principal. Further, $J\cap\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x(C_{i})]=\langle\overline{\Phi(C_{i})}\rangle$ and, therefore, $J_{C_{i}}=\langle\overline{\Phi(C_{i})}\rangle$. ∎

Proposition 4.14 shows that the positive part of the steady state variety can easily be projected onto coordinate subspaces by dropping circuits. Hence, the PSS matroid allows for some freedom to “pick and mix” variables according to measurements. The picking and mixing corresponds to the geometric operation of projection of the variety $X=V(\langle\overline{\Phi(C_{1})},\dots,\overline{\Phi(C_{n-d})}\rangle)\subseteq\mathbb{R}(\kappa_{1},\dots,\kappa_{m})^{n}$. Next, suppose there exists a measurement $\xi=\{\xi_{i_{1}},\dots,\xi_{i_{\ell}+d}\}$ containing measurements for the species in a basis $S$ of the PSS and the species in the circuits $C_{1},\dots,C_{\ell}$, all containing $S$. Combining the idea of projection and measurement (evaluation) leads to the definition of a parameter variety.

Definition 4.15**.**

Keeping the same notation as above and, by choosing an appropriate set of generators, i.e. “clearing the denominators”, let $J=\langle\overline{\Phi(C_{1})},\dots,\overline{\Phi(C_{\ell})}\rangle\subseteq\mathbb{R}[\kappa_{1},\dots,\kappa_{m},x_{1},\dots,x_{n}]$. Hence, $V(J)\subset\mathbb{R}^{m}\times\mathbb{R}^{n}$. The parameter variety, $X_{m}$, is obtained from $V(J)$ by the evaluation

[TABLE]

The parameter variety is obtained by the selective projection and evaluation of the binomials obtained from Definition 4.5 and is a variety in the space of parameters only. Every parameter vector compatible with the measurement $\xi$ is on the parameter variety. Every set of measurements gives rise to a parameter variety and, in order to be compatible with a model, their intersection needs to be non empty, in fact, their intersection needs to be non empty in the positive orthant. In order to uniquely identify a model based on a given measurement the positive orthant of $X_{m}$ needs to consist of a single point only. Various algebraic techniques can be applied to show when this is the case, e.g. **[45, 53]**.

The parameter varieties of families of toric networks can be related by applications of projections (selective application of Definition 4.5) and the partial evaluation map $\pi$. Suppose $J_{N}$ and $J_{\text{int}}$ are the ideals of the $N^{\text{th}}$ member of a family and the intermediate model between the $N^{\text{th}}$ and $(N+1)^{\text{th}}$ member, respectively. Let both ideals (or a projection of them) contain the same circuits $C_{1},\dots,C_{\ell}$. Then, by Proposition 3.15, $J_{N}=\pi(J_{\text{int}})$.

5 Inheritance of Multistationarity for Toric Families

In this section we investigate the inheritance of multistationarity among members of families of toric chemical reaction networks. Our main result is Theorem 5.2, where we apply results of **[24, 53]** to show that if we can find a multistationary member of a family (satisfying certain conditions), then every larger member of the same family is multistationary for some parameter values. We begin by introducing some notation.

Definition 5.1**.**

Let $I_{B}\subset\mathbb{R}[x_{1},\dots,x_{n}]$ be a prime binomial ideal defining a complete intersection of codimension $n-d$ and let the matrix $B$ be an exponent matrix of a minimal generating set of $I_{B}$. Define

[TABLE]

where $Z$ is an integer matrix. Further, let

[TABLE]

where $\left(B^{T}\right)^{\lambda}=(b_{1}\lambda_{1},\dots,b_{n}\lambda_{n})$ for the columns $b_{i}$ of $B^{T}$. We call $J^{\lambda}$ regular if $\det(J^{\lambda})\neq 0$ for some values of $(\lambda_{1},\dots,\lambda_{n})$.

We now state the main theorem of this section.

Theorem 5.2**.**

Fix a family of toric chemical reaction networks. Suppose that the family obeys the following conditions.

1. (C1)

The family has toric intermediates. 2. (C2)

The number of conservation relations stays constant. 3. (C3)

The matrix $J^{\lambda}$ for the $N^{\text{th}}$ network is regular.

Then, if the $N^{\text{th}}$ member of the family is capable of multistationarity, every member of the family for which $M\geq N$ is also capable of multistationarity.

Before proving the theorem above we prove the following lemma.

Lemma 5.3**.**

Fix a family of reaction networks, where each member of the family and each intermediate network have the same number of conservation relations. Let $Z_{N}$ be a conservation relation matrix of the $N^{\text{th}}$ network. Then, $Z_{N+1}$ is obtained by adding columns to $Z_{N}$, i.e. $Z_{N}$ is a submatrix of $Z_{N+1}$.

Proof.

First, build the stoichiometric matrix of the intermediate network by adding $n^{\prime}$ species and $m^{\prime}$ reactions to the $N^{\text{th}}$ network. Hence, the stoichiometric matrix of the intermediate network, $\Gamma^{\text{int}}_{N}$, has the form

[TABLE]

where $\Gamma_{N}=(r_{1},\dots,r_{m})$ is the $n\times m$ stoichiometric matrix of the $N^{\text{th}}$ network and, $\tilde{\Gamma}_{N}=(r_{m+1},\dots,r_{m+m^{\prime}})$ is the $(n+n^{\prime})\times m^{\prime}$ matrix of new reactions. Hence, any basis of the left kernel of $\Gamma^{\text{int}}_{N}$ can be written as $\tilde{z}_{i}=(z_{i}\>,\>z^{\prime}_{i})$, where $z_{i}$ is an element of a basis of the left kernel of $\Gamma_{N}$ or zero and $z^{\prime}_{i}$ is a vector to be determined by $\tilde{\Gamma}_{N}$. The assumption of equidimensionality of the left kernels of $\Gamma_{N}$ and $\Gamma_{N}^{\text{int}}$ can be satisfied only if $\text{rank}(\tilde{\Gamma}_{N})=n^{\prime}$ and, in the same fashion as Corollary 3.18, it follows, that $Z_{N}$ is a submatrix of $Z_{N}^{\text{int}}$. Now, delete $m^{\prime\prime}$ of the first $m$ columns of $\Gamma_{N}^{\text{int}}$ to form $\Gamma_{N+1}$. By the equidimensionality assumption it follows, that the rows of $Z_{N}^{\text{int}}$ form a basis of the left kernel of $\Gamma_{N+1}$. Hence, $Z_{N}^{\text{int}}=Z_{N+1}$ and $Z_{N}$ is a submatrix of $Z_{N+1}$. ∎

From Lemma 5.3 and Theorem 4.8 it becomes apparent why (C2) and toric steady states are required. Namely, toric steady states are needed to construct a matrix $B$ and, hence, $J_{N}^{\lambda}$ and (C2) is required to ensure that the conservation relation matrices are submatrices. The proof of Theorem 5.2 is now stated.

Proof.

Fix a vector $x^{*}\in\mathbb{R}_{>0}^{n}$ as in Definition 4.4 and choose a basis $S$ and $n-d$ circuits $\{C_{1},\dots,C_{n-d}\}$ of the PSS matroid associated to the $N^{\text{th}}$ member of a family such that $\bigcap_{i}C_{i}=S$. Consider the polynomial system

[TABLE]

Further, let $B_{N}$ denote the exponent matrix of the binomials $\overline{\Phi(C_{1})},\dots,\overline{\Phi(C_{n-d})}$. Hence, following the construction of Definition 5.1 we obtain the square matrix

[TABLE]

By construction (see Theorem 4.8), $V(\overline{\Phi(C_{1})}(x),\cdots,\overline{\Phi(C_{n-d})}(x))\cap\mathbb{R}_{>0}^{n}\neq\emptyset$. By [24, Theorem 2.7] the system (15) is multistationary if and only if either $\det(J^{\lambda}_{N})=0$ or $\det(J^{\lambda}_{N})\neq 0$ and the polynomial $\det(J^{\lambda}_{N})$ in $\lambda_{1},\dots,\lambda_{n}$ has a positive and a negative term. Suppose the $N^{\text{th}}$ network is multistationary and, by condition (C3), $\det(J^{\lambda}_{N})\neq 0$. Next, build the $(N+1)^{\text{th}}$ network by adding $n^{\prime}$ new species and consider its matrix $J_{N+1}^{\lambda}$. To obtain this matrix consider the new circuit binomials $\overline{\Phi(C_{n-d+1})}(x)=\cdots=\overline{\Phi(C_{n+n^{\prime}-d})}(x)=0$, where each circuit contains the basis $S$ and one new variable $x_{n+i}$ with its corresponding positive exponent $b_{n+i}$, where $i=1,\dots,n^{\prime}$. The exponent vectors will therefore have exactly one non-zero element (namely the positive integer $b_{n+i}$) in the last $n^{\prime}$ rows. Therefore, the matrix $J_{N+1}^{\lambda}$ has the block form

[TABLE]

where $A|_{[y_{1}\dots y_{m}\times a_{1}\dots a_{n}]}$ denotes the restriction of a matrix $A$ to the rows $y_{1}\dots y_{m}$ and columns $a_{1}\dots a_{n}$. The matrix ${\rm diag}(b_{n+1}\lambda_{n+1},\dots,b_{n+n^{\prime}}\lambda_{n+n^{\prime}})$ is a diagonal matrix with diagonal elements $b_{n+1}\lambda_{n+1},\dots,b_{n+n^{\prime}}\lambda_{n+n^{\prime}}$. Hence, the expression $T=\left(\prod_{i=n+1}^{n+n^{\prime}}b_{i}\lambda_{i}\right)\det(J^{\lambda}_{N})\neq 0$ must appear in $\det(J^{\lambda}_{N+1})$ and, in particular, no term in $T$ can be cancelled by any other term appearing in $\det(J^{\lambda}_{N+1})$. Hence, if $\det(J^{\lambda}_{N})$ has coefficients of opposite sign so does $\det(J^{\lambda}_{N+1})$ and, therefore the network is multistationary. The proof is completed by induction. ∎

Remark 5.4**.**

Note that, since we are allowed to choose different bases of the PSS matroid, and hence, different binomial systems it may be possible to satisfy condition (C3) of Theorem 5.2 for one particular choice of basis but not for a different choice of basis.

We illustrate our results using the two and three site distributive phosphorylation networks.

Example 5.5**.**

The PSS matroid of the one-site and two-site distributive phosphorylation networks are represented by the $A$-matrices

[TABLE]

respectively. Hence, choosing a basis of $a_{1}=(1,1,0)^{T}$, $a_{2}=(0,1,0)^{T}$ and $a_{3}=(0,0,1)^{T}$ we find a parameterization for the one-site network as

[TABLE]

The two-site model has three additional equations, namely

[TABLE]

Hence, we get the $B$-matrices

[TABLE]

The conservation relations can be chosen as

[TABLE]

This leads to

[TABLE]

Note that this determinant is square-free and homogeneous as investigated in [23] and that it has only coefficients equal to $+1$. Indeed, it is a well known fact the the one-site distributive network is monostationary [54]. The determinant of $J^{\lambda}_{2}$ is equal to

[TABLE]

and, therefore, contains a term $T=\lambda_{7}\lambda_{8}\lambda_{9}\text{det}(J^{\lambda}_{1})$. The determinant of the two-site network has coefficients of opposite signs and, hence, the network is multistationary and so are all larger networks in the family. In particular, it is known that the $N$-site distributive network has a maximum of $2N-1$ positive steady states [54].

6 Conclusion

In this paper we studied families of chemical reaction networks with toric steady states, which we called toric families. First, we investigated the dimensions and parameterizations of toric steady state varieties and connected them to network properties whenever possible. In particular, if there is a finite number of positive steady states for all choices of $c_{1},\dots,c_{d}$, then the number of conservation relations determines the dimension of the toric steady state variety and, with certain restrictions, the monomial parameterization of a chemical species $X_{i}$ is preserved throughout the family.

We next studied the PSS matroid defined by the parameterization of the positive steady states. In particular, we proved that the algebraic matroid defined by the binomial steady state ideal is equivalent to the PSS matroid. We showed how the PSS matroid can be decorated with binomials and how they can be used for model selection, experimental design or even parameter identification.

*The final section investigated the multistationarity structure of toric families. The main result of the section showed that, under some mild restrictions, if a member of a family is capable of multistationarity then all larger members are too. This result was proved using the binomials constructed from the PSS matroids. We illustrated our results on the multisite distributive phosphorylation network. *

Further research could include applying the results of this paper to other meaningful biochemical families such as different models for immune system reactions, e.g. **[1, 55]**. Another direction could be to study the parameter varieties defined in this paper in the context of previous identifiability research and aim to include noisy data.

Acknowledgements

The authors would like to thank the American Institute of Mathematics (AIM) for initiating this collaboration, the Erwin Schödinger Institute, Vienna (ESI) and the Institute for Computational and Experimental Research in Mathematics (ICERM) where part of the research was carried out. MFA would like to thank the EPSRC for supporting this research through grant EP/G03706X/1 and Merton College Oxford. MH was partially supported by the Independent Research Fund of Denmark. The authors would like to than Elisenda Feliu and Heather Harrington for helpful comments on the manuscript.

7 Appendix A: Additional Proofs

7.1 Proof of Lemma 3.13

Proof.

Note that $\overline{V(L)-V(x_{1}\cdots x_{n+n^{\prime}})}=V(L)$. Hence,

[TABLE]

which proves the first statement. For the second statement note that, by construction of the evaluation map, $X_{N+1}\cap V(L)=X_{N}\subset\mathbb{C}^{m}_{\kappa}\times\mathbb{C}^{n}_{x}\times\mathbb{C}^{n^{\prime}}$. Hence,

[TABLE]

∎

7.2 Proof of Proposition 3.15

Proof.

By assumption $I_{N+1}:(x_{1}\cdots x_{n+n^{\prime}})^{\infty}$ is a binomial ideal, hence,

[TABLE]

where $\Xi_{j}^{\pm}$ are polynomial functions of the rate constants. For each term $\beta_{j}=\Xi_{j}^{+}x^{b_{j}^{+}}-\Xi_{j}^{-}x^{b_{j}^{-}}$ it holds that $\pi(\beta_{j})=\pi(\Xi_{j}^{+})x^{b_{j}^{+}}-\pi(\Xi_{j}^{-})x^{b_{j}^{-}}$, hence either:

• •

$\pi(\beta_{j})=0$, this happens if $\pi(\Xi_{j}^{+})=\pi(\Xi_{j}^{-})=0$;

• •

$\pi(\beta_{j})$ is a monomial, this happens if one of $\pi(\Xi_{j}^{+})$ or $\pi(\Xi_{j}^{-})$ evaluates to zero and the other does not;

• •

$\pi(\beta_{j})=\beta^{\prime}_{j}$ and still binomial, this happens when both $\pi(\Xi^{+}_{j})\neq 0$ and $\pi(\Xi^{-}_{j})\neq 0$.

If $\pi(\beta_{j})$ is a monomial for any $j$ then the ideal $\pi(I_{N+1}):(x_{1}\cdots x_{n+n^{\prime}})^{\infty}=(1)$, otherwise, it must be a binomial ideal generated by the $\beta_{j}$ such that $\pi(\beta_{j})=\beta^{\prime}_{j}$. By Lemma 3.13 this ideal is the binomial ideal $I_{N}:(x_{1}\cdots x_{n})^{\infty}$.

Since for a fixed choice of rate constants $I_{N}:(x_{1}\cdots x_{n})^{\infty}$ is a binomial ideal in $\mathbb{R}[x_{1},\dots,x_{n}]$ and it is assumed that $V(I_{N})$ has a nonempty intersection with the positive orthant $\mathbb{R}_{>0}^{n}$ then the last statement follows by [46, Proposition 5.2]. This Proposition states that if $I$ is a binomial ideal with variety $V(I)$ at most one of its irreducible components intersects $\mathbb{R}^{n}_{>0}$. ∎

7.3 Proof of Proposition 3.17

Proof.

First statement (i) is proved. By assumption $I_{N+1}:(x_{1}\cdots x_{n+n^{\prime}})^{\infty}$ is a binomial ideal, hence,

[TABLE]

By Proposition 3.15 it follows that if $I_{N}:(x_{1}\cdots x_{n})^{\infty}$ is not empty, then it is also a binomial ideal and its set of generators must appear in $\pi(\Xi_{1}^{+}x^{b_{1}^{+}}-\Xi_{1}^{-}x^{b_{1}^{-}},\dots,\Xi_{\nu}^{+}x^{b_{\nu}^{+}}-\Xi_{\nu}^{-}x^{b_{\nu}^{-}})$. Hence, its matrix of exponents can be obtained by choosing a subset (of size $\mu\leq\nu$) of the exponent vector pairs $({b_{1}^{+}},{b_{1}^{-}}),\dots,({b_{\nu}^{+}},{b_{\nu}^{-}})$. By definition of the matrix $B$, it follows that $B_{N}$ is a submatrix of $B_{N+1}$, proving (i).

Note that no variables in $x_{n+1},\dots,x_{n+n^{\prime}}$ appear in $I_{N}$ and, hence, the subset of exponent vectors $b_{i}=b_{i}^{+}-b_{i}^{-}$ defining $I_{N}:(x_{1}\cdots x_{n})^{\infty}$ have the general form $b_{i}=(b_{i_{1}},\dots,b_{i_{n}},0,\dots,0)^{T}$. Further, the entries of the remaining $\nu-\mu$ exponent vectors have at least one non-zero entry in the last $n^{\prime}$ rows and these vectors can be collected into a matrix $\tilde{B}$. Hence, $B_{N+1}$ may be written in block form as

[TABLE]

It is clear that the first $n$ entries of any basis vector of the left kernel of $B_{N+1}$ must form a vector in the left kernel of $B_{N}$. Hence, the left kernel of $B_{N+1}$ has a basis of the form $\tilde{a}_{i}=(a_{i},\bar{a}_{i})$, where $a_{i}$ is an element of a given basis of the left kernel of $B_{N}$ or the zero vector and $\bar{a}_{i}\in\mathbb{Z}^{n^{\prime}}$ is a vector of variables to be determined. From the rank condition it is clear that the matrix $\tilde{B}$ has a rank of at most $n^{\prime}$ and, for any given $i$, a system of linear equations is obtained, $\tilde{a}_{i}\tilde{b}=0$, for all columns $\tilde{b}$ of $\tilde{B}$. Due to the rank condition, this system has always at least one (rational) solution and, hence, the collection of all solutions form a basis of the left kernel of $B_{N+1}$ to give the matrix $A_{N+1}$ (after clearing the denominators by multiplications of rows with scalars). Therefore, $A_{N}$ can be chosen as a submatrix of $A_{N+1}$, proving (ii). ∎

7.4 Proof of Theorem 3.19

Proof.

Part (i) follows immediately from Theorem 3.17. To prove (ii), suppose $I_{N}$ and $I_{N+1}$ are the (binomial) steady state ideals of two members of a family and $I^{\prime}_{N+1}$ is the steady state ideal of their intermediate network. It follows from Theorem 3.17 that the $A$-matrices of $I_{N}$ and $I_{N+1}$ are submatrices of the $A$-matrix of the intermediate. Each column of the $A$-matrix corresponds to a variable $x_{i}$. Since both $I^{\prime}_{N+1}$ and $I_{N+1}$ are ideals with species $\{x_{1},\dots,x_{n+n^{\prime}}\}$ then their associated $A$-matrices have the same number of columns. Further, by assumption, the number of conservation relations stays constant. Therefore, by Proposition 2.14 the $A$ matrices associated to $I^{\prime}_{N+1}$ and $I_{N+1}$ are the same, since the ideals have the same dimension, meaning the $A$-matrices must have the same number of rows. Hence, the $A$-matrices associated to $I_{N}$ and $I_{N+1}$ are submatrices; i.e. $A_{N}$ is a submatrix of $A_{N+1}$. ∎

7.5 Proof of Theorem 4.3

Proof.

The variables of the ground set of $\mathcal{M}(A)$ have algebraic dependencies as defined in Proposition 4.1. Hence, their algebraic dependencies can be found by solving the implicitization problem $I=J\cap\mathbb{R}(\kappa_{1},\dots,\kappa_{m})[x_{1},\dots,x_{n}]$ where $J=\langle x_{1}-x^{*}_{1}t^{a_{1}},\dots,x_{n}-x^{*}_{n}t^{a_{n}}\rangle$. However, the same ideal defines the implicit equations of a toric variety defined by $\psi_{A}:t\to(x^{*}_{1}t^{a_{1}},\dots,x^{*}_{n}t^{a_{n}})$. Hence, the algebraic relations between the monomials $\phi_{A}(x)$ are identical to the algebraic relations defined by the binomial ideal $I_{b}$. This implies that $\mathcal{M}(A)=\mathcal{M}(I_{b})$. ∎