# Adding species to chemical reaction networks: preserving rank preserves
nondegenerate behaviours

Murad Banaji  Balázs Boros  Josef Hofbauer  Department of Design Engineering
and Mathematics, Middlesex University London  Department of Mathematics,
University of Vienna

###### Abstract

We show that adding new chemical species into the reactions of a chemical
reaction network (CRN) in such a way that the rank of the network remains
unchanged preserves its capacity for multiple nondegenerate equilibria and/or
periodic orbits. One consequence is that any bounded nondegenerate behaviours
which can occur in a CRN can occur in a CRN with bounded stoichiometric
classes. The main result adds to a family of theorems which tell us which
enlargements of a CRN preserve its capacity for nontrivial dynamical
behaviours. It generalises some earlier claims, and complements similar claims
involving the addition of reactions into CRNs. The result gives us information
on how ignoring some chemical species, as is common in biochemical modelling,
might affect the allowed dynamics in differential equation models of CRNs. We
demonstrate the scope and limitations of the main theorem via several
examples. These illustrate how we can use the main theorem to predict
multistationarity and oscillation in CRNs enlarged with additional species;
but also how the enlargements can introduce new behaviours such as additional
periodic orbits and new bifurcations.

###### keywords:

Oscillation; multistationarity; chemical reaction networks; bifurcations MSC.
92E20; 37C25; 34D10

##  1 Introduction and outline of the main result

An important theme in the mathematical study of chemical reaction networks
(CRNs) relates to how network structure influences network dynamics. The
results in this direction sometimes allow us to infer detailed information on
dynamical behaviours of reaction networks using only graph theory and linear
algebra. We may, for example, be able to conclude from basic computations that
a CRN has very simple behaviour, such as the convergence of all initial
conditions to an equilibrium, regardless of parameter values. In the opposite
direction, we may be able to state, without numerical simulation, that a CRN
admits some interesting behaviour such as stable oscillation, and even know a
priori how to choose parameter values to obtain this behaviour.

Amongst more complicated behaviours which occur in ordinary differential
equation models of CRNs, the most well-studied are multistationarity and
oscillation. The history of study of these behaviours in the context of
biological modelling is reviewed in [40]. Crucially, multistationarity and
oscillation are not just of abstract interest, but may be of functional
importance in biological switching and signalling processes [24, 28, 31, 33,
29, 12, 30]. For this reason, results which tell us which structures in a CRN
guarantee such behaviours are of considerable interest.

A family of theorems termed “inheritance results” tell us when a CRN is
guaranteed to exhibit some dynamical behaviour simply because of the presence
of a certain subnetwork. Examples of such results can be found in [23, 18, 8,
3, 4]. Such conditions which guarantee nontrivial behaviours based on CRN
structure are dual to claims which rule out nontrivial behaviours in a CRN.
There is a large classical and modern literature on conditions which preclude
certain dynamical behaviours in CRNs. Examples include [22, 16, 2, 13, 37, 6,
17, 1]. Each new result in either direction helps to narrow the gap between
conditions guaranteeing, and conditions ruling out, nontrivial behaviours in
CRNs.

The main result of this paper is a new inheritance result, which is both
natural and relatively straightforward to prove. The theorem is most simply
stated in terms of the rank of a CRN. Each reaction of a CRN on $n$ chemical
species defines a reaction vector, a real (often integer) $n$-vector whose
$k$th entry tells us the net production of species $k$ in the reaction. The
span of these vectors is the stoichiometric subspace of the CRN, whose
dimension is defined to be the rank of the CRN. In any model of the CRN,
numbers or concentrations of the chemical species are confined to affine
subspaces parallel to the stoichiometric subspace. The nonnegative portions of
these affine subspaces are termed the stoichiometric classes of the system.

The rank of a CRN figures in various aspects of the theory. For example, it
plays a crucial role in the original results of deficiency theory [16], and
indeed in the definition of the deficiency of a reaction network. More
recently, various results have been proved for CRNs of sufficiently low rank,
regardless of how many species or reactions are involved. Characterisations of
rank-1 CRNs admitting multistationarity are given in [27]. In rank-2 CRNs, the
Poincaré-Bendixson theorem [41, Chapter 9] can be used to rule out chaos or
guarantee oscillation; and the same can sometimes be extended to rank-3 CRNs,
a fact exploited in the analysis in [39] and [11]. In [32], the famous “global
attractor conjecture” is proved for rank-3 CRNs.

Stated informally, we have the following complementary inheritance results
involving the rank of a CRN. Both are simple applications of regular
perturbation theory:

  1. 1.

Adding new reactions into a CRN without changing its rank preserves its
capacity for nontrivial behaviours including nondegenerate multistationarity
and oscillation. This was proved in previous work (Theorem 1 in [8] and
Theorem 1 in [3]).

  2. 2.

Adding new species into a CRN without changing its rank preserves its capacity
for nontrivial behaviours including nondegenerate multistationarity and
oscillation. This is the content of Theorem 1 here.

These claims demonstrate that in CRN theory we often have complementary
results where we can interchange “species” and “reactions”. Note, however,
that adding species into a CRN while preserving its rank can result in some
fairly fundamental changes to the CRN. For example, stoichiometric classes of
the enlarged CRN may be bounded even if those of the original were unbounded.
This is in contrast to adding linearly dependent reactions, a process which
leaves stoichiometric classes unchanged.

Taken together, the above claims imply that building a CRN without altering
its rank preserves its capacity for nondegenerate dynamical behaviours. This
claim is stated more formally as Corollary 1 later. Note that the enlarged
CRNs may, of course, admit more complicated and interesting behaviours than
the original, as we shall see by example.

After some preliminary definitions, we will present the statement and proof of
the main theorem and several remarks on its implications and generalisations.
This is followed by several examples which demonstrate both the main result
and its limitations.

##  2 Preliminaries

We present some key notions briefly. A much more expansive treatment of the
main background can be found in previous work [7, 8, 3]. We consider a CRN to
be an ordered set of chemical species and an ordered set of reactions, where
the orderings are arbitrary but fixed. Each reaction is an ordered pair of
complexes, namely formal linear combinations of chemical species. The
coefficient of each species in a complex is taken to be nonnegative, and often
to be an integer, although the latter is not required here. The pair of
complexes which define a reaction are termed the “reactant complex” and
“product complex” of the reaction.

Positive sets in Euclidean space. The positive orthant in $\mathbb{R}^{n}$ is
denoted by $\mathbb{R}^{n}_{+}$ and defined as
$\\{x\in\mathbb{R}^{n}\,:\,x_{i}>0\,\,\mbox{for}\,\,i=1,\ldots,n\\}$. A set in
$\mathbb{R}^{n}$ is termed positive if it lies in $\mathbb{R}^{n}_{+}$. The
closure of $\mathbb{R}^{n}_{+}$ is denoted by $\mathbb{R}^{n}_{\geq 0}$ and
referred to as the nonnegative orthant in $\mathbb{R}^{n}$, namely
$\mathbb{R}^{n}_{\geq 0}=\\{x\in\mathbb{R}^{n}\,:\,x_{i}\geq
0\,\,\mbox{for}\,\,i=1,\ldots,n\\}$.

Stoichiometric matrix and stoichiometric classes. Each reaction is associated
with a reaction vector whose $k$th entry is the net production of the $k$th
species in the reaction: this is just the stoichiometric coefficient of the
$k$th species in the product complex minus its stoichiometric coefficient in
the reactant complex. The stoichiometric matrix of a CRN is the matrix whose
$j$th column is the $j$th reaction vector of the CRN. Given a CRN
$\mathcal{R}$ on $n$ species with stoichiometric matrix $\Gamma$, the span of
the columns of $\Gamma$ (a linear subspace of $\mathbb{R}^{n}$) is denoted by
$\mathrm{im}\,\Gamma$, and is termed the stoichiometric subspace of
$\mathcal{R}$. The rank of $\Gamma$ is termed the rank of $\mathcal{R}$. The
nonnegative parts of cosets of $\mathrm{im}\,\Gamma$, namely sets of the form
$(x+\mathrm{im}\,\Gamma)\cap\mathbb{R}^{n}_{\geq 0}$
($x\in\mathbb{R}^{n}_{\geq 0}$), are the stoichiometric classes of
$\mathcal{R}$. The positive parts of stoichiometric classes, namely sets of
the form $(x+\mathrm{im}\,\Gamma)\cap\mathbb{R}^{n}_{+}$
($x\in\mathbb{R}^{n}_{+}$), are the positive stoichiometric classes of
$\mathcal{R}$.

Ordinary differential equations (ODEs) and rate functions. Any system of ODEs
describing the evolution of a CRN with stoichiometric matrix $\Gamma$ takes
the form $\dot{x}=\Gamma v(x)$. The function $v$ is termed the rate function
of the CRN and its $j$th component tells us how the rate of the $j$th reaction
depends on the concentrations of the chemical species. In this paper, $v$ is
assumed, as a minimum, to be defined and continously differentiable on
$\mathbb{R}^{n}_{+}$.

Nondegenerate and linearly stable limit sets. Consider a CRN of rank $r$ with
stoichiometric matrix $\Gamma$. Let $\mathcal{S}$ be some coset of
$\mathrm{im}\,\Gamma$ containing a positive equilibrium (resp., periodic
orbit) $\mathcal{O}$. Since the positive part of $\mathcal{S}$ is locally
invariant, $\mathcal{O}$ has exactly $r$ eigenvalues (resp., Floquet
multipliers) relative to $\mathcal{S}$. If none of these are equal to zero
(resp., exactly one of these is equal to one), then we say that $\mathcal{O}$
is nondegenerate. (Note that this terminology differs from that in [3], where
an invariant set was referred to as “nondegenerate” only if it was hyperbolic
relative to its stoichiometric class.) In an abuse of terminology, we refer to
$\mathcal{O}$ as hyperbolic if it is hyperbolic relative to its stoichiometric
class, and linearly stable if it is linearly stable relative to its
stoichiometric class.

Continuation of nondegenerate limit sets. Suppose $U$ is some open region in
$\mathbb{R}^{n}$, $a>0$, and $f:U\times(-a,a)\to\mathbb{R}^{n}$ is $C^{1}$. If
$\dot{x}=f(x,0)$ has a nondegenerate equilibrium (resp., periodic orbit) on
$U$, then the same is true for $\dot{x}=f(x,\varepsilon)$ for all
$\varepsilon$ sufficiently small. The same conclusion holds if we replace
“nondegenerate” by “hyperbolic” or “linearly stable”. The results follow from
the implicit function theorem and the fact that eigenvalues of a matrix depend
continuously on its entries. The claim for equilibria is an immediate
consequence of the implicit function theorem (see, for example, Corollary 6.9
in [8]), while the details for periodic orbits are laid out in Section IV in
[20]. An immediate consequence is that nondegenerate equilibria and periodic
orbits of $\dot{x}=\Gamma v(x)$ on some positive stoichiometric class survive
sufficiently small $C^{1}$ perturbations to $v(x)$ (here $v$ is assumed to be
at least $C^{1}$). In fact, if we restrict attention to hyperbolic equilibria
and periodic orbits, then the results for equilibria and periodic orbits are
special cases of more general results on the persistence of normally
hyperbolic invariant manifolds in [19] and [21]. Thus, the main result here
generalises naturally to the case of invariant manifolds which are normally
hyperbolic relative to their stoichiometric class.

Entrywise products and generalised monomials. The notation $a\circ b$ denotes
the entrywise product of matrices or vectors $a$ and $b$ assumed to have the
same dimensions. Given $x=(x_{1},\ldots,x_{n})^{\mathrm{t}}$ and
$a=(a_{1},\ldots,a_{n})$, $x^{a}$ is an abbreviation for the generalised
monomial $x_{1}^{a_{1}}x_{2}^{a_{2}}\cdots x_{n}^{a_{n}}$. Let
$A_{1},\ldots,A_{m}$ be the rows of an $m\times n$ matrix $A$. Each
$x^{A_{i}}$ is then a generalised monomial, and $x^{A}$ denotes the vector of
these monomials, namely,
$(x^{A_{1}},x^{A_{2}},\ldots,x^{A_{m}})^{\mathrm{t}}$.

Kinetics and admitted behaviours. When we restrict the rate function of a CRN
$\mathcal{R}$ to some class of functions $\mathcal{K}$, the pair
$(\mathcal{R},\mathcal{K})$ is referred to as a CRN with kinetics. We can
think of $(\mathcal{R},\mathcal{K})$ as a set of allowed ODE models of the
CRN. We say that $(\mathcal{R},\mathcal{K})$ admits some particular dynamical
behaviour if this behaviour occurs in some allowed model, i.e., for some
choice of rate function from $\mathcal{K}$ and on some stoichiometric class.
Otherwise the CRN with kinetics forbids this behaviour. Different classes of
kinetics for CRNs are discussed in detail in [3].

Power-law kinetics. Let $\mathsf{X}_{1},\ldots,\mathsf{X}_{n}$ denote the
chemical species of a CRN and $x_{1},\ldots,x_{n}$ denote the concentrations
of these species. If the $i$th reaction has power-law kinetics, this means
that the $i$th rate function takes the form $v_{i}(x)=\kappa_{i}x^{a}$, where
$\kappa_{i}$ is a positive constant (termed the rate constant of the $i$th
reaction), and $a$ is a real (row) vector, termed the vector of exponents for
the reaction. If all reactions of a CRN have power-law kinetics, we can stack
these row vectors into a matrix $A$, termed the matrix of exponents of the
reaction network, whose $ij$th entry tells us the exponent of species $j$ in
the rate function for reaction $i$. In this case, the rate function can be
written briefly $\kappa\circ x^{A}$. If $A$ is fixed in advance, we say that
the CRN has fixed power-law kinetics.

Mass action kinetics. Mass action kinetics is a special case of fixed power-
law kinetics where $a_{ij}$, the $ij$th entry in the matrix of exponents, is
precisely the stoichiometric coefficient of species $\mathsf{X}_{j}$ in the
reactant complex of reaction $i$.

Enlarging CRNs and inheritance. We are often interested in claiming that
whenever a CRN with kinetics admits some behaviour, then so does an enlarged
CRN with kinetics in some related class. Inheritance results relate dynamical
behaviours admitted in the enlarged CRNs to those admitted in the original
CRNs.

Adding linearly dependent species. In this paper we are interested in an
enlargement where $\mathcal{R}^{\prime}$ is obtained from $\mathcal{R}$ by
adding a linearly dependent species into the reactions of $\mathcal{R}$. In
such a modification, the stoichiometric matrix of the CRN is unchanged except
for the addition of a new row, and this new row is a linear combination of the
existing rows of the stoichiometric matrix. We refer to the process of
enlarging a CRN $\mathcal{R}$ by adding in some new, linearly dependent,
species as lifting. The terminology is motivated by the fact that the addition
of a new species increases the dimension of the state space by $1$, with the
original, lower dimensional, state space naturally embedded in the new state
space.

Derived power-law kinetics. Suppose that $(\mathcal{R},\mathcal{K})$ and
$(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ are two CRNs, with
$\mathcal{R}^{\prime}$ obtained by adding new species and/or reactions to
$\mathcal{R}$, and $\mathcal{K}$, $\mathcal{K}^{\prime}$ being fixed power-law
kinetics with matrices of exponents $A$ and $A^{\prime}$ respectively. We then
say that $\mathcal{K}^{\prime}$ is derived from $\mathcal{K}$, if the
submatrix of $A^{\prime}$ corresponding to the original species and reactions
of $\mathcal{R}$ is precisely $A$.

Permanence. Consider a system of ODEs on some subset of $\mathbb{R}^{n}$ and
suppose that $\mathcal{X}\subseteq\mathbb{R}^{n}_{\geq 0}$ is forward
invariant for the system. The system is permanent on $\mathcal{X}$ if there
exists a forward invariant, compact, positive set
$\mathcal{Z}\subseteq\mathcal{X}$ such that the forward trajectory of every
positive initial condition in $\mathcal{X}$ eventually enters $\mathcal{Z}$.
In the context of CRNs we may think of $\mathcal{X}$ as a stoichiometric
class: we may be interested in permanence on some or all stoichiometric
classes.

##  3 The main result

Given a parameterised family of compact sets $\mathcal{X}_{\varepsilon}$ in
Euclidean space, with $\varepsilon\in(0,a)$ for some $a>0$,
“$\mathcal{X}_{\varepsilon}$ is close to $\mathcal{X}$” will mean that given
any $\delta>0$, there exists $\varepsilon_{1}\in(0,a]$ such that for all
$\varepsilon\in(0,\varepsilon_{1})$ the Hausdorff distance between
$\mathcal{X}_{\varepsilon}$ and $\mathcal{X}$ is less than $\delta$.

######  Theorem 1.

Let $(\mathcal{R},\mathcal{K})$ be a CRN with fixed power-law kinetics. Let
$\mathcal{R}^{\prime}$ be derived by adding to $\mathcal{R}$ a new linearly
dependent species, and let $\mathcal{K}^{\prime}$ be any fixed power-law
kinetics for $\mathcal{R}^{\prime}$ derived from $\mathcal{K}$. Suppose that
for some choice of kinetics from $\mathcal{K}$, $\mathcal{R}$ has, on some
stoichiometric class, at least $0\leq r_{1}<\infty$ positive, nondegenerate
(resp., hyperbolic, resp., linearly stable) equilibria and at least $0\leq
r_{2}<\infty$ positive, nondegenerate (resp., hyperbolic, resp., linearly
stable) periodic orbits. Then, for some choice of kinetics from
$\mathcal{K}^{\prime}$, $\mathcal{R}^{\prime}$ has, on some stoichiometric
class, at least $r_{1}$ positive, nondegenerate (resp., hyperbolic, resp.,
linearly stable) equilibria and at least $r_{2}$ positive, nondegenerate
(resp., hyperbolic, resp., linearly stable) periodic orbits.

###### Proof.

Let $\mathcal{R}$ include $n$ chemical species and $m$ irreversible reactions
with stoichiometric matrix $\Gamma\in\mathbb{R}^{n\times m}$. By hypothesis,
there exists $v\in\mathcal{K}$, such that the associated ODE system

$\dot{x}=\Gamma v(x)\,$ (1)

has $r_{1}$ nondegenerate equilibria and $r_{2}$ nondegenerate periodic orbits
on some positive stoichiometric class $\mathcal{S}_{+}$. The rate function $v$
is analytic on $\mathbb{R}^{n}_{+}$, and hence certainly $C^{1}$ on
$\mathbb{R}^{n}_{+}$ which is, in fact, the only assumption we will need about
$v$. Let $\mathcal{S}$ be the affine hull of $\mathcal{S}_{+}$, i.e., the
coset of $\mathrm{im}\,\Gamma$ containing $\mathcal{S}_{+}$.

Let $\mathcal{O}$ refer to a nondegenerate positive equilibrium (resp.,
periodic orbit) on $\mathcal{S}$. Choose $\mathcal{Z}\subseteq\mathcal{S}$ to
be compact and positive, with $\mathcal{O}\subseteq\mathcal{Z}^{o}$, the
relative interior of $\mathcal{Z}$ in $\mathcal{S}$. We may assume that
similar sets are constructed around each of the $r_{1}+r_{2}$ nondegenerate
equilibria and periodic orbits of (1), and that these sets are pairwise
disjoint.

The hypothesis that the added species is linearly dependent implies that there
exists $c\in\mathbb{R}^{n}$ such that the stoichiometric matrix of
$\mathcal{R}^{\prime}$ takes the form

$\Gamma^{\prime}:=\left(\begin{array}[]{c}\Gamma\\\
c^{t}\Gamma\end{array}\right)\,.$

Define $\varepsilon$ to be a positive parameter to be controlled, and let
$\alpha_{j}\in\mathbb{R}\,\,(j=1,\ldots,m)$ be any real numbers. Denoting the
concentration of the new species by $y$, set the rate of the $j$th reaction in
$\mathcal{R}^{\prime}$ to be

$v^{\prime}_{j}(x,y,\varepsilon)=\varepsilon^{\alpha_{j}}\,v_{j}(x)\,y^{\alpha_{j}}\,.$

This choice of reaction rate corresponds to giving the new species exponent
$\alpha_{j}$ in reaction $j$, and multiplying the original rate constant of
the $j$th reaction by $\varepsilon^{\alpha_{j}}$. In brief notation,
$\mathcal{R}^{\prime}$ now has rate function
$v^{\prime}(x,y,\varepsilon):=\varepsilon^{\alpha}\circ v(x)\circ y^{\alpha}$
where $\alpha:=(\alpha_{1},\ldots,\alpha_{m})^{t}$. The evolution of
$\mathcal{R}^{\prime}$ is governed by

$\left(\begin{array}[]{c}\dot{x}\\\
\dot{y}\end{array}\right)=\Gamma^{\prime}\,v^{\prime}(x,y,\varepsilon)\,.$ (2)

Note that $\mathcal{R}^{\prime}$ has a new conservation law of the form
$-c^{t}x+y=\mathrm{constant}$. For any fixed $\varepsilon>0$, we will focus
our attention on the invariant set
$\mathcal{H}_{\varepsilon}\subseteq\mathbb{R}^{n+1}$ defined by setting this
constant to be equal to $\frac{1}{\varepsilon}$, namely,

$\mathcal{H}_{\varepsilon}=\\{(x,y)\in\mathbb{R}^{n}\times\mathbb{R}\,\colon\,y=\frac{1}{\varepsilon}+c^{t}x\\}\,.$

The map

$h_{\varepsilon}\colon x\mapsto\left(\begin{array}[]{c}x\\\
\frac{1}{\varepsilon}+c^{t}x\end{array}\right)$

is an affine bijection between $\mathbb{R}^{n}$ and
$\mathcal{H}_{\varepsilon}$. If $c=0$, set $\varepsilon_{1}=1$, and otherwise
set

$\varepsilon_{1}=\frac{1}{\sup_{x\in\mathcal{Z}}\,|c^{t}x|}\,.$

Then, for $\varepsilon\in(0,\varepsilon_{1})$, $h_{\varepsilon}(\mathcal{Z})$
is a compact, positive subset of $\mathcal{H}_{\varepsilon}$. The map
$h_{\varepsilon}$ defines local coordinates on $\mathcal{H}_{\varepsilon}$
which evolve according to

$\dot{x}=\Gamma\,v^{\prime}(x,\frac{1}{\varepsilon}+c^{t}x,\varepsilon)=\Gamma\,(\varepsilon^{\alpha}\circ
v(x)\circ(\frac{1}{\varepsilon}+c^{t}x)^{\alpha})=\Gamma\,(v(x)\circ(1+\varepsilon\,c^{t}x)^{\alpha})\,.$
(3)

Here $x$ refers to the local coordinate on $\mathcal{H}_{\varepsilon}$, rather
than the original coordinate on $\mathbb{R}^{n}$: this should cause no
confusion as we are identifying $\mathcal{H}_{\varepsilon}$ with
$\mathbb{R}^{n}$ via $h_{\varepsilon}$. Note that the right hand side of (3)
is well-defined (and $C^{1}$) provided $x$ is positive, and
$\varepsilon<\frac{1}{|c^{t}x|}$, which certainly holds on
$\mathcal{Z}^{o}\times(-\varepsilon_{1},\varepsilon_{1})$.

We wish to restrict our attention to $\mathcal{Z}^{o}$. We can, if desired,
pass to local coordinates on $\mathcal{S}$ in a standard way (see the proofs
of several results in [3]), but here this is unnecessary: we simply bear in
mind that we are considering the restriction of (3) to $\mathcal{Z}^{o}$, with
$\varepsilon\in(-\varepsilon_{1},\varepsilon_{1})$.

Since the vector field in (3) is a $C^{1}$ perturbation of that of (1), by
regular perturbation theory, there exists
$\varepsilon_{2}\in(0,\varepsilon_{1}]$ such that for each
$\varepsilon\in(0,\varepsilon_{2})$, (3) has an equilibrium (resp., periodic
orbit) $\mathcal{O}_{\varepsilon}$ in $\mathcal{Z}^{o}$, which is
nondegenerate, and close to $\mathcal{O}$. (The details in the harder case
where $\mathcal{O}$ is a periodic orbit are in Section IV in [20], for
example.) If $\mathcal{O}$ is hyperbolic relative to $\mathcal{S}$, then we
can choose $\varepsilon_{2}$ to ensure that the linear stability type of
$\mathcal{O}_{\varepsilon}$ relative to $\mathcal{S}$ is the same as that of
$\mathcal{O}$. More precisely, (i) if $\mathcal{O}$ is an equilibrium with,
relative to $\mathcal{S}$, $k_{1}$ eigenvalues with positive real part,
$k_{2}$ eigenvalues with negative real part, and no eigenvalues on the
imaginary axis, then the same holds for $\mathcal{O}_{\varepsilon}$; (ii) if
$\mathcal{O}$ is a periodic orbit with, relative to $\mathcal{S}$, $k_{1}$
Floquet multipliers inside the unit circle, $k_{2}$ outside the unit circle,
and precisely one multiplier on the unit circle, then the same holds for
$\mathcal{O}_{\varepsilon}$. As a special case, if $\mathcal{O}$ was linearly
stable relative to $\mathcal{S}$, then the same holds for
$\mathcal{O}_{\varepsilon}$.

Since $\mathcal{O}_{\varepsilon}$ lies in $\mathcal{Z}^{o}$, it is positive.
Let $\mathcal{S}^{\prime}_{\varepsilon}=h_{\varepsilon}(\mathcal{S})$.
Explicitly,

$\mathcal{S}^{\prime}_{\varepsilon}:=\left(\begin{array}[]{c}x_{0}\\\
\frac{1}{\varepsilon}+c^{t}x_{0}\end{array}\right)+\mathrm{im}\,\Gamma^{\prime}\,,$

where $x_{0}$ is any element of $\mathcal{S}$. Note that
$\mathcal{S}^{\prime}_{\varepsilon}$ has the same dimension as $\mathcal{S}$.
Clearly,
$\mathcal{O}^{\prime}_{\varepsilon}:=h_{\varepsilon}(\mathcal{O}_{\varepsilon})$
is an equilibrium (resp., periodic orbit) of (2), and we have ensured (via the
choice of $\varepsilon_{1}$) that $h_{\varepsilon}(\mathcal{Z}^{o})$, and
hence $\mathcal{O}^{\prime}_{\varepsilon}\subseteq
h_{\varepsilon}(\mathcal{Z}^{o})$, are positive. As $h_{\varepsilon}$ is an
affine bijection between $\mathcal{S}$ and
$\mathcal{S}^{\prime}_{\varepsilon}$, the choice of $\varepsilon_{2}$ ensures
that the linear stability type of $\mathcal{O}^{\prime}_{\varepsilon}$
relative to $\mathcal{S}^{\prime}_{\varepsilon}$ is the same as that of
$\mathcal{O}_{\varepsilon}$ relative to $\mathcal{S}$.

We can repeat the same argument in a neighbourhood of each of a finite number
of nondegenerate equilibria or periodic orbits of (1) on $\mathcal{S}_{+}$. By
choosing $\varepsilon_{2}^{*}$ to be the minimum of the values of
$\varepsilon_{2}$ associated with each limit set, we can ensure that provided
$\varepsilon\in(0,\varepsilon_{2}^{*})$, $\mathcal{R}^{\prime}$ has at least
$r_{1}$ positive, nondegenerate equilibria and at least $r_{2}$ positive,
nondegenerate periodic orbits on $\mathcal{S}^{\prime}_{\varepsilon}$.
Moreover, whenever one of the original limit sets was hyperbolic relative to
$\mathcal{S}$, we can ensure that the lifted limit set is of the same linear
stability type relative to $\mathcal{S}^{\prime}_{\varepsilon}$. This
completes the proof. ∎

Several remarks are in order.

######  Remark 1 (Mass action kinetics).

The result clearly holds if we insist that both $\mathcal{R}$ and
$\mathcal{R}^{\prime}$ have mass action kinetics which is simply a special
case of fixed power-law kinetics. In this case, in the proof of Theorem 1,
$\alpha_{j}$ is the stoichiometric coefficient of the new species in the
reactant complex of the $j$th reaction of $\mathcal{R}^{\prime}$.

######  Remark 2 (The proof is constructive).

The proof of Theorem 1, as with other inheritance results based on
perturbation theory, is constructive. It tells us how to set rate constants
and how to choose a stoichiometric class in order to find the desired
behaviour in the enlarged CRN $\mathcal{R}^{\prime}$.

######  Remark 3 (The projected dynamics are close to the original).

Consider some “lifted” bounded orbit of $\mathcal{R}^{\prime}$ such as
$\mathcal{O}^{\prime}_{\varepsilon}$ in the proof of Theorem 1. Its projection
$\mathcal{O}_{\varepsilon}$ onto $x$ coordinates, can be made as close as we
desire to the original orbit of $\mathcal{R}$ (namely, $\mathcal{O}$) by
choosing $\varepsilon$ to be small. But this comes at the cost of large values
of the new species concentration $y$ on the lifted orbit, and small rate
constants. The next remark indicates the limitations of the lifting process.

######  Remark 4 (We cannot always control the lifted dynamics over an entire
stoichiometric class).

The proof of Theorem 1 tells us the following: given any positive
stoichiometric class, say $\mathcal{S}_{+}$, of the original CRN
$\mathcal{R}$, fixing the perturbation parameter $\varepsilon$ at any positive
value selects a positive stoichiometric class, say
$(\mathcal{S}^{\prime}_{\varepsilon})_{+}$, of the lifted CRN
$\mathcal{R}^{\prime}$. Assume that rate constants are fixed and let $V$ and
$V^{\prime}_{\varepsilon}$ refer to the original and lifted vector fields on
$\mathcal{S}_{+}$ and $(\mathcal{S}^{\prime}_{\varepsilon})_{+}$ respectively.
Choosing $\varepsilon$ to be small ensures that the projection of
$V^{\prime}_{\varepsilon}$ onto $\mathcal{S}_{+}$ is close to the original
vector field $V$ on $\mathcal{S}_{+}$ on that portion of
$\mathcal{S}^{\prime}_{\varepsilon}$ where the concentration of the added
species (denoted by $y$ in the proof) is large. But, regardless of how small
we choose $\varepsilon$ to be, if there are regions of
$(\mathcal{S}^{\prime}_{\varepsilon})_{+}$ where $y$ is small, then in these
regions $V^{\prime}_{\varepsilon}$ need not be close to $V$. The consequences
are illustrated in the example of the Brusselator in Section 4.4, where the
lifting process leads to a loss of permanence on every stoichiometric class.
Note, however, that if $\mathcal{S}_{+}$, the original positive stoichiometric
class of $\mathcal{R}$, is itself bounded, then we can control the lifted
vector field over the entirety of $(\mathcal{S}^{\prime}_{\varepsilon})_{+}$.

######  Remark 5 (Normally hyperbolic invariant manifolds persist).

Although Theorem 1 is phrased in terms of equilibria and periodic orbits, the
result admits generalisation. Indeed, with the assumptions of the theorem, if
$\mathcal{O}$ is any positive, compact, invariant manifold admitted by
$\mathcal{R}$ and normally hyperbolic relative to its stoichiometric class
$\mathcal{S}$, then it survives $C^{1}$ perturbations [19, 21], and hence is
admitted by $\mathcal{R}^{\prime}$. If, for example, $\mathcal{R}$ admits a
$k$-dimensional torus on some positive stoichiometric class, and the torus is
normally hyperbolic relative to this class, then the same holds for
$\mathcal{R}^{\prime}$.

######  Remark 6 (Bifurcations persist).

Suppose $\mathcal{R}$ admits, on some positive stoichiometric class, a
nondegenerate local bifurcation of an equilibrium or periodic orbit, unfolded
nondegenerately by the rate constants. Then, for sufficiently small, fixed,
$\varepsilon>0$, $\mathcal{R}^{\prime}$ admits the same nondegenerate local
bifurcation on some positive stoichiometric class as we vary the same
combination of rate constants (note that rate constants of $\mathcal{R}$ and
$\mathcal{R}^{\prime}$ are in natural one-to-one correspondence). Essentially,
the nondegeneracy and transversality conditions associated with the
bifurcation, allow us to continue the bifurcation as we vary $\varepsilon$.
Moreover, these conditions continue to hold for sufficiently small
$\varepsilon$. For a concrete example demonstrating the pesistence of a
bifurcation, see Section 4.2.

######  Remark 7 (A generalisation of previous claims).

For CRNs with power-law kinetics, Theorem 1 generalises the claims in Theorems
3 in [8] and [3], which treat the very special case where the added species
figures only trivially in reactions, i.e., adds only a row of zeros to the
stoichiometric matrix of the network. Note that, in that case, the new
stoichiometric classes were bounded if and only if the original stoichiometric
classes were bounded.

######  Remark 8 (Generalisations to other classes of kinetics).

Phrasing the result in terms of power-law kinetics simplifies the proof, but
is not key to it. The broad template of the proof can be applied to CRNs with
other classes of kinetics.

######  Remark 9 (CRNs with bounded stoichiometric classes do not have greatly
restricted dynamics).

One immediate consequence of Theorem 1 is that insisting a CRN has bounded
stoichiometric classes, does not greatly restrict its behaviour. If a given
CRN with unbounded stoichiometric classes admits some finite set of bounded
nondegenerate limit sets on one of its stoichiometric classes, then we can
always construct, by adding in a dependent species, a CRN with bounded
stoichiometric classes which admits the same bounded nondegenerate limit sets
on one of its stoichiometric classes. We see several instances of this in the
examples presented in Section 4.

The next corollary tells us that it may be helpful to examine full-rank
subnetworks of a CRN: finding nontrivial behaviours in these subnetworks is
sufficient to ensure that they occur in the original CRN. We recall the
definition of an induced subnetwork of a CRN from [3]: this is a CRN obtained
by removing some reactions from a CRN, and/or some species from all the
reactions in which they figure. In terms of the Petri-net graph of the CRN
this corresponds to removing some vertices from the graph along with all their
incident arcs.

######  Corollary 1.

Let $\mathcal{R}$ be a CRN of rank $r$, and let $\mathcal{R}_{0}$ be any
rank-$r$ induced subnetwork of $\mathcal{R}$. If, for some fixed power-law
kinetics, $\mathcal{R}_{0}$ admits $k_{1}$ positive nondegenerate (resp.,
hyperbolic, resp., linearly stable) equilibria and $k_{2}$ positive
nondegenerate (resp., hyperbolic, resp., linearly stable) periodic orbits,
then the same holds for $\mathcal{R}$ with any derived power-law kinetics.

###### Proof.

Clearly $\mathcal{R}$ can be built from $\mathcal{R}_{0}$ by adding linearly
dependent reactions and linearly dependent species to $\mathcal{R}_{0}$. The
result is thus an immediate consequence of Theorem 1 above which deals with
the case of adding linearly dependent species, and Theorems 1 in [8] and [3]
which deal with the case of adding linearly dependent reactions. (Although
Theorems 1 in [8] and [3] are stated in more restricted terms, the
generalisations required are immediate.) ∎

##  4 Examples

First, we introduce some terminology relevant to the examples below.

Homogeneous CRNs. The _molecularity_ of a complex in a CRN is the sum of its
stoichiometric coefficients. We call a CRN _homogeneous_ if, for every
reaction, the molecularities of the reactant complex and the product complex
are equal. Clearly, this condition is equivalent to
$(1,1,\ldots,1)^{\mathrm{t}}$ being an element of the kernel of
$\Gamma^{\mathrm{t}}$. In particular, the stoichiometric classes of a
homogeneous CRN are bounded: for the ODE associated with a homogeneous CRN,
regardless of the precise nature of the kinetics, the quantity
$x_{1}+x_{2}+\cdots+x_{n}$ is conserved. A partial converse is also true: it
is easily shown that if a CRN is endowed with any fixed power-law kinetics,
and $x_{1}+x_{2}+\cdots+x_{n}$ is constant along trajectories for some open
set of rate constants, then the CRN is homogeneous.

Homogenisation of CRNs. Starting with an arbitrary CRN $\mathcal{R}$, one can
make it homogeneous by adding a new species with appropriate stoichiometric
coefficients to the reactant or product complex of each reaction [14, Exercise
4 on page 29]. This operation preserves the rank of $\mathcal{R}$, and the
homogenisation can be carried out in multiple ways. By Theorem 1, if
$\mathcal{R}$ has mass action kinetics, and nondegenerate multistationarity
(resp., oscillation) occurs in $\mathcal{R}$, then it also occurs in the
homogenised CRN.

In all of the examples below, we start with a network that is not homogeneous,
and then homogenise it. The examples in Sections 4.1 and 4.2 illustrate our
main result, while those in Sections 4.3 and 4.4 demonstrate its limitations.
The rank of each of these networks is one or two. For the homogenisation of
some rank-three mass-action systems, consult [11]. For applications of Theorem
1 where the enlarged network is not homogeneous, but nevertheless has bounded
stoichiometric classes, see [10].

Mathematica code for the analysis of all the examples is available on GitHub
[9].

###  4.1 Schlögl model: a single-species CRN with multiple equilibria

We use the reversible version of the Schlögl model [34, 35] to demonstrate the
use of Theorem 1 for guaranteeing the existence of multiple nondegenerate
equilibria in an enlarged reaction network. It also provides some insight into
the proof of Theorem 1.

Consider the following single-species mass action system and its associated
differential equation:

$\mathsf{0}$$\mathsf{X}$$2\mathsf{X}$$3\mathsf{X}$$6$$11$$6$$1$$\dot{x}=-x^{3}+6x^{2}-11x+6.$

It has $3$ positive nondegenerate equilibria: at $x=1$, $x=2$, and $x=3$. The
first and third are linearly stable, while the one at $x=2$ is linearly
unstable.

We now homogenise the network. Note that the simplest choice of homogenisation
gives the enlarged network
$\mathsf{Y}\rightleftharpoons\mathsf{X},\,\,2\mathsf{X}+\mathsf{Y}\rightleftharpoons
3\mathsf{X}$. However, we choose the slightly more complicated homogenised
network
$2\mathsf{Y}\rightleftharpoons\mathsf{X}+\mathsf{Y},\,\,2\mathsf{X}+\mathsf{Y}\rightleftharpoons
3\mathsf{X}$ in order to demonstrate two points: that there are many ways to
homogenise a CRN; and that the modified rate constants given in the proof of
Theorem 1 can have nonlinear dependence on a perturbation parameter
$\varepsilon$ as seen below.

Theorem 1 now tells us that this network must admit three nondegenerate
equilibria on some stoichiometric class, two linearly stable and one linearly
unstable. To see why, we follow the proof of Theorem 1 and obtain the
following mass action system and its associated differential equation,
dependent on a new parameter $\varepsilon$:

$2\mathsf{Y}$$\mathsf{X}+\mathsf{Y}$$2\mathsf{X}+\mathsf{Y}$$3\mathsf{X}$$6\varepsilon^{2}$$11\varepsilon$$6\varepsilon$$1$$\begin{aligned}
\dot{x}&=-x^{3}+6\varepsilon x^{2}y-11\varepsilon xy+6\varepsilon^{2}y^{2},\\\
\dot{y}&=+x^{3}-6\varepsilon x^{2}y+11\varepsilon
xy-6\varepsilon^{2}y^{2}.\end{aligned}$

In the homogenised system, the stoichiometric subspace remains
$1$-dimensional, and the quantity $x+y$ is conserved. We now restrict
attention to the stoichiometric class defined by $x+y=\frac{1}{\varepsilon}$
and replace $y$ by $\frac{1}{\varepsilon}-x$. The dynamics of $x$ for
$0<x<\frac{1}{\varepsilon}$ is then given by

$\displaystyle\dot{x}=-x^{3}+6x^{2}(1-\varepsilon x)-11x(1-\varepsilon
x)+6(1-\varepsilon x)^{2}\,.$

On any compact subinterval of $(0,\infty)$, the vector field
$-x^{3}+6x^{2}(1-\varepsilon x)-11x(1-\varepsilon x)+6(1-\varepsilon x)^{2}$
converges uniformly to $-x^{3}+6x^{2}-11x+6$ as $\varepsilon\to 0$. It is not
hard to see that for all sufficiently small $\varepsilon$ the lifted system
with the scaled rate constants has $3$ positive equilibria in the
stoichiometric class $x+y=\frac{1}{\varepsilon}$, two of which are stable,
while one is unstable. Note, however, that for any fixed rate constants, each
stoichiometric class defined by $x+y=c$ with $c>0$ being large enough, has a
unique positive equilibrium. This illustrates that in order to obtain the
desired behaviour we must simultaneously choose the rate constants and a
stoichiometric class of the lifted system.

###  4.2 Lotka-Volterra-Autocatalator

Consider the Lotka-Volterra-Autocatalator [15, (8)] with one reversible
reaction [38, Section 3] (henceforth abbreviated to LVA), and its associated
mass action differential equation:

$2\mathsf{X}$$3\mathsf{X}$$\mathsf{X}+\mathsf{Y}$$2\mathsf{Y}$$\mathsf{Y}$$\mathsf{0}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\kappa_{4}$$\begin{aligned}
\dot{x}&=\kappa_{1}x^{2}-\kappa_{2}x^{3}-\kappa_{3}xy,\\\
\dot{y}&=\kappa_{3}xy-\kappa_{4}y.\end{aligned}$

The system has exactly one positive equilibrium if and only if
$\kappa_{1}\kappa_{3}>\kappa_{2}\kappa_{4}$ (otherwise there is no positive
equilibrium). We can confirm that at
$\kappa_{1}\kappa_{3}=2\kappa_{2}\kappa_{4}$ the positive equilibrium
undergoes a supercritical Andronov-Hopf bifurcation, and thus, for
$\kappa_{1}\kappa_{3}-2\kappa_{2}\kappa_{4}$ positive and sufficiently small,
the system has a stable periodic orbit.

Let us now homogenise the LVA. The resulting network and its associated mass
action differential equation take the form

$2\mathsf{X}+\mathsf{Z}$$3\mathsf{X}$$\mathsf{X}+\mathsf{Y}$$2\mathsf{Y}$$\mathsf{Y}$$\mathsf{Z}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\kappa_{4}$$\begin{aligned}
\dot{x}&=\kappa_{1}x^{2}z-\kappa_{2}x^{3}-\kappa_{3}xy,\\\
\dot{y}&=\kappa_{3}xy-\kappa_{4}y,\\\
\dot{z}&=\kappa_{2}x^{3}-\kappa_{1}x^{2}z+\kappa_{4}y.\end{aligned}$

Theorem 1 tells us that the homogenised LVA must admit a linearly stable
periodic orbit on some stoichiometric class. Let us use these systems to
demonstrate the arguments in the proof of Theorem 1 and some subsequent
remarks. Consider again the ODE systems associated with the LVA,

$\displaystyle\dot{x}$
$\displaystyle=\kappa_{1}x^{2}-\kappa_{2}x^{3}-\kappa_{3}xy,$ (4)
$\displaystyle\dot{y}$ $\displaystyle=\kappa_{3}xy-\kappa_{4}y,$

and with the homogenised LVA, where we have now added a prime ′ to the rate
constants to distinguish them from those of the original LVA,

$\displaystyle\dot{x}$
$\displaystyle=\kappa^{\prime}_{1}x^{2}z-\kappa^{\prime}_{2}x^{3}-\kappa^{\prime}_{3}xy,$
(5) $\displaystyle\dot{y}$
$\displaystyle=\kappa^{\prime}_{3}xy-\kappa^{\prime}_{4}y,$
$\displaystyle\dot{z}$
$\displaystyle=\kappa^{\prime}_{2}x^{3}-\kappa^{\prime}_{1}x^{2}z+\kappa^{\prime}_{4}y.$

If we set
$(\kappa^{\prime}_{1},\kappa^{\prime}_{2},\kappa^{\prime}_{3},\kappa^{\prime}_{4})=(\varepsilon\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})$,
and restrict attention to the stoichiometric class on which
$x+y+z=\frac{1}{\varepsilon}$, then, using this equation to eliminate $z$, we
find that $x$ and $y$ in the homogenised LVA evolve according to

$\displaystyle\dot{x}$
$\displaystyle=\kappa_{1}x^{2}-\kappa_{2}x^{3}-\kappa_{3}xy-\varepsilon\kappa_{1}(x^{3}+x^{2}y),$
(6) $\displaystyle\dot{y}$ $\displaystyle=\kappa_{3}xy-\kappa_{4}y\,.$

We see immediately that the vector field in (6) converges on compact sets to
that in (4) as $\varepsilon\to 0$.

In particular, suppose we fix some values of
$\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4}$ such that
$\kappa_{1}\kappa_{3}-2\kappa_{2}\kappa_{4}$ is positive and sufficiently
small to ensure that (4) has a linearly stable, positive periodic orbit
$\mathcal{O}$. Then, for sufficiently small $\varepsilon>0$, (6) has a
linearly stable, positive, periodic orbit $\mathcal{O}_{\varepsilon}$.
Moreover, provided $\varepsilon$ is sufficiently small,
$z=\frac{1}{\varepsilon}-x-y$ remains positive as $x$ and $y$ vary along this
periodic orbit, and so the lifted system (5) has a positive periodic orbit
$\mathcal{O}^{\prime}_{\varepsilon}$ on the stoichiometric class defined by
$x+y+z=\frac{1}{\varepsilon}$. Clearly, since $\mathcal{O}_{\varepsilon}$ is
linearly stable for (6), $\mathcal{O}^{\prime}_{\varepsilon}$ is linearly
stable relative to this class.

In fact, we can go further and argue that the bifurcation which gives rise to
stable periodic orbits in (4) itself survives lifting (see Remark 6). Suppose
we fix some path in parameter space, $P\colon\mathbb{R}\to\mathbb{R}^{4}_{+}$,
$s\mapsto(\kappa_{1}(s),\kappa_{2}(s),\kappa_{3}(s),\kappa_{4}(s))$ such that
$P(0)$ lies on the surface in $\mathbb{R}^{4}_{+}$ defined by
$\kappa_{1}\kappa_{3}-2\kappa_{2}\kappa_{4}=0$, and $P$ is transverse to this
surface at $0$. This ensures that (4) undergoes a nondegenerate supercritical
Andronov-Hopf bifurcation at $s=0$. Then we can be sure that for any
sufficiently small $\varepsilon>0$, a nondegenerate supercritical Andronov-
Hopf bifurcation occurs in (5) on the stoichiometric class defined by
$x+y+z=\frac{1}{\varepsilon}$ as we traverse the corresponding path
$P^{\prime}$ in
$(\kappa^{\prime}_{1},\kappa^{\prime}_{2},\kappa^{\prime}_{3},\kappa^{\prime}_{4})$
space obtained by noting that
$(\kappa^{\prime}_{1},\kappa^{\prime}_{2},\kappa^{\prime}_{3},\kappa^{\prime}_{4})=(\varepsilon\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})$.
This bifurcation occurs for some value of $s$ close to $0$.

###  4.3 Lotka reactions

We now provide an example where the assumptions of Theorem 1 are violated, and
periodic solutions are not preserved by lifting.

Consider the Lotka reactions and the associated mass action differential
equation:

$\mathsf{X}$$2\mathsf{X}$$\mathsf{X}+\mathsf{Y}$$2\mathsf{Y}$$\mathsf{Y}$$\mathsf{0}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\begin{aligned}
\dot{x}&=\kappa_{1}x-\kappa_{2}xy,\\\
\dot{y}&=\kappa_{2}xy-\kappa_{3}y.\end{aligned}$

In this case, the unique positive equilibrium
$\left(\frac{\kappa_{3}}{\kappa_{2}},\frac{\kappa_{1}}{\kappa_{2}}\right)$ is
surrounded by a continuum of periodic orbits; these are level sets of the
nonlinear first integral $x^{-\kappa_{3}}y^{-\kappa_{1}}e^{\kappa_{2}(x+y)}$.
Since these periodic orbits are degenerate, Theorem 1 does not apply.

Indeed, adding a new species, $\mathsf{Z}$, to some of the reactions while
preserving the rank of the network can lead to the destruction of all of the
periodic orbits. Consider the following network and its associated mass action
differential equation:

$\mathsf{X}+\mathsf{Z}$$2\mathsf{X}$$\mathsf{X}+\mathsf{Y}$$2\mathsf{Y}$$\mathsf{Y}$$\mathsf{Z}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\begin{aligned}
\dot{x}&=\kappa_{1}xz-\kappa_{2}xy,\\\ \dot{y}&=\kappa_{2}xy-\kappa_{3}y,\\\
\dot{z}&=-\kappa_{1}xz+\kappa_{3}y.\end{aligned}$

The set of positive equilibria is
$\left\\{\left(\frac{\kappa_{3}}{\kappa_{2}},\frac{\kappa_{1}}{\kappa_{2}}t,t\right)\colon
t>0\right\\}$. Thus, the stoichiometric classes $x+y+z=c$ with
$c\leq\frac{\kappa_{3}}{\kappa_{2}}$ have no positive equilibria, while those
with $c>\frac{\kappa_{3}}{\kappa_{2}}$ have a unique positive equilibrium.
Note that the divergence of the vector field after division by $xyz$ equals
$-\frac{\kappa_{3}}{xz^{2}}$. Since this quantity is negative on
$\mathbb{R}^{3}_{+}$, there is no periodic orbit that lies entirely in the
positive orthant [36, Satz 1] (see also Remark (v) following the proof of
Theorem 2.3 in [26]). In fact, every positive equilibrium is globally
asymptotically stable within its positive stoichiometric class: this follows
via the Poincaré-Bendixson Theorem, since it can be shown that on any
stoichiometric class which includes a positive equilibrium no positive initial
condition has omega limit set intersecting the boundary of the nonnegative
orthant.

###  4.4 Brusselator

Our final example demonstrates that while the lifted CRN must admit the
nondegenerate behaviours of the original CRN, it may also allow other
behaviours not seen in the original CRN, such as multiple periodic orbits and
homoclinic orbits. Looked at from another angle, omitting a single species
from a CRN, even without changing its rank, can result in the loss of many
different nontrivial behaviours.

Consider the Brusselator and its associated mass action differential equation:

$\mathsf{0}$$\mathsf{X}$$\mathsf{Y}$$2\mathsf{X}+\mathsf{Y}$$3\mathsf{X}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\kappa_{4}$$\begin{aligned}
\dot{x}&=\kappa_{1}-\kappa_{2}x-\kappa_{3}x+\kappa_{4}x^{2}y,\\\
\dot{y}&=\kappa_{3}x-\kappa_{4}x^{2}y.\end{aligned}$

At any fixed values of the rate constants the system has a unique positive
equilibrium,
$(x^{*},y^{*})=\left(\frac{\kappa_{1}}{\kappa_{2}},\frac{\kappa_{2}\kappa_{3}}{\kappa_{1}\kappa_{4}}\right)$.
It can be shown that the system is permanent. Moreover, it is known that
$(x^{*},y^{*})$ is globally asymptotically stable for
$\kappa_{3}\leq\kappa_{2}+\frac{\kappa_{1}^{2}\kappa_{4}}{\kappa_{2}^{2}}$,
while it is repelling for
$\kappa_{3}>\kappa_{2}+\frac{\kappa_{1}^{2}\kappa_{4}}{\kappa_{2}^{2}}$ and is
surrounded by a periodic orbit that is born via a supercritical Andronov-Hopf
bifurcation. Moreover, this periodic orbit is unique and attracts every
positive initial condition, except $(x^{*},y^{*})$ (see [42, Example 5 on page
135]).

Let us now homogenise the Brusselator. The resulting CRN and its mass action
differential equation take the form

$\mathsf{Z}$$\mathsf{X}$$\mathsf{Y}$$2\mathsf{X}+\mathsf{Y}$$3\mathsf{X}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\kappa_{4}$$\begin{aligned}
\dot{x}&=\kappa_{1}z-\kappa_{2}x-\kappa_{3}x+\kappa_{4}x^{2}y,\\\
\dot{y}&=\kappa_{3}x-\kappa_{4}x^{2}y,\\\
\dot{z}&=-\kappa_{1}z+\kappa_{2}x.\end{aligned}$

By Theorem 1, there exist rate constants and a stoichiometric class where the
new system has a stable periodic orbit. In fact, given any rate constants
satisfying
$\kappa_{3}>\kappa_{2}+\frac{\kappa_{1}^{2}\kappa_{4}}{\kappa_{2}^{2}}$ and
such that the Brusselator has a linearly stable periodic orbit, the proof of
Theorem 1 gives us a procedure for finding a linearly stable periodic orbit in
the lifted system. Moreover, by Remark 6, a supercritical Andronov-Hopf
bifurcation must occur on some stoichiometric class as rate constants are
varied in the lifted system.

Note, however, that the global behaviour of the homogenised Brusselator is
quite different from that of the Brusselator. The first interesting difference
is that whereas the original system was permanent, the lifted system is not
permanent on any stoichiometric class intersecting the positive orthant; we
show in A that for all rate constants and all positive values of the parameter
$c$, the stoichiometric class satisfying $x+y+z=c$ includes a boundary
equilibrium $(0,c,0)$ which is asymptotically stable relative to this class.
Thus we cannot control the vector field at all points on any of these
stoichiometric classes, even though these classes are compact (see Remark 4
following the proof of Theorem 1).

Aside from the loss of permanence, the homogeneous Brusselator admits several
nontrivial behaviours forbidden in the original. In A we show that for the
homogenised network both supercritical and subcritical Andronov-Hopf
bifurcations can occur on some stoichiometric classes as we vary rate
constants. Furthermore, various interesting co-dimension two bifurcations take
place: a generic Bautin bifurcation [25, Section 8.3] and a generic Bogdanov-
Takens bifurcation [25, Section 8.4] can both occur. This implies that the
lifted system admits, for various choices of the rate constants, the following
behaviours not seen in the original Brusselator:

  * 1.

a fold bifurcation of equilibria;

  * 2.

an unstable positive equilibrium surrounded by a stable and an unstable
periodic orbit;

  * 3.

a fold bifurcation of periodic orbits;

  * 4.

a stable positive equilibrium surrounded by an unstable periodic orbit or a
homoclinic orbit;

  * 5.

a homoclinic bifurcation.

It is not unexpected that the lifted system admits richer dynamical behaviour;
but it is worth noting that if we restrict attention to any stoichiometric
class and consider the evolution of concentrations of $\mathsf{X}$ and
$\mathsf{Y}$, then the (2-dimensional) lifted vector field differs from the
original only by the addition of linear terms (see equation (7) in A). The
addition of linear terms to a 2D differential equation can thus quite
dramatically increase the complexity of its behaviour.

##  5 Discussion and conclusions

Inheritance results tell us how we might enlarge a CRN while preserving its
capacity for various dynamical behaviours. The main result in this paper adds
a simple but important inheritance result relevant to the study of both
multistationarity and oscillation in CRNs. Some inheritance results, including
the main result of this paper, are gathered in [5]. Taken together, these
results provide a powerful tool for predicting nontrivial behaviours in a CRN
based on its subnetworks.

It is useful to think of inheritance results in terms of partial orders. For
example, suppose a CRN $\mathcal{R}$ with mass action kinetics admits linearly
stable oscillation. The same then holds for all CRNs greater than
$\mathcal{R}$ in the partial order defined by available inheritance results.
We may then look for CRNs admitting linearly stable oscillation which are
minimal with respect to this partial order in order to gain insight into the
capacity for stable oscillation in larger and more complex CRNs [3]. Such a
program provides a rigorous basis for claims about “motifs”, namely small
subnetworks which are at the root of certain behaviours in biological systems.
Systematically identifying minimal CRNs with prescribed behaviours, followed
by the development of algorithms to test for their presence in larger CRNs, is
a natural avenue for future work.

While inheritance results are often phrased in terms of enlarging CRNs, they
can also be seen in terms of modelling choices. The main result here is
relevant to choices which might affect conservation laws. In physically
realistic systems of chemical reactions occurring in a closed environment we
expect numbers of atoms of each element to be conserved: this is the so-called
law of atomic balance [14]. However, it is common in modelling CRNs to omit
some species from reactions, particularly when they are considered to be
present in abundance, or their concentration is subject to external control.
As an example, when reactions involving ATP and ADP occur in biochemical
models, inorganic phosphate and water are often omitted from the equations.
Such omissions tend to destroy physical conservation laws.

It is natural to worry that omitting species from reactions might introduce
fundamentally new, and unrealistic, behaviours into the system. Theorems 4 in
[8] and [3] provide some reassurance that this will not occur if we omit
species whose concentration is sufficiently strongly controlled by external
processes. In a similar way, Theorem 1 here provides some reassurance in the
case where the omitted species are linearly dependent on the others. In this
case, the omissions cannot introduce the capacity for behaviours such as
stable oscillation or multistability. In fact, any behaviour of the simplified
CRN occurring on a compact set, and which is robust in the sense that it
survives $C^{1}$ perturbation, can also be obtained on some stoichiometric
class of the larger CRN for appropriate choices of parameters. We can, of
course, still lose interesting behaviours by omitting dependent species, as
illustrated by the example of the Brusselator in Section 4.4.

The examples demonstrated both the power and limitations of the result.
Crucially, we only expect nondegenerate behaviours to survive the addition of
new species as illustrated by the example in Section 4.3. And the result is
not global: introducing new linearly dependent species may introduce new
global behaviours regardless of how small we make the perturbation parameter
$\varepsilon$. In particular, we saw, in the example of the Brusselator
(Section 4.4), how the lifted system could allow positive trajectories to
converge to the boundary of state space, even though the original system had a
positive, globally attracting, compact set for all choices of parameters. The
details are laid out in A.

##  Appendix A Dynamics of the homogenised Brusselator

Consider the homogenised Brusselator and its associated mass action
differential equation:

$\mathsf{Z}$$\mathsf{X}$$\mathsf{Y}$$2\mathsf{X}+\mathsf{Y}$$3\mathsf{X}$$\kappa_{1}$$\kappa_{2}$$\kappa_{3}$$\kappa_{4}$$\begin{aligned}
\dot{x}&=\kappa_{1}z-\kappa_{2}x-\kappa_{3}x+\kappa_{4}x^{2}y,\\\
\dot{y}&=\kappa_{3}x-\kappa_{4}x^{2}y,\\\
\dot{z}&=-\kappa_{1}z+\kappa_{2}x.\end{aligned}$

The stoichiometric classes which intersect the positive orthant are
$\mathcal{P}_{c}=\\{(x,y,z)\in\mathbb{R}^{3}_{\geq 0}\colon x+y+z=c\\}$ for
$c>0$. Note that for each positive $c$, the corner $(0,c,0)$ of the triangle
$\mathcal{P}_{c}$ is an equilibrium of this system (and this is the only
boundary equilibrium in $\mathcal{P}_{c}$).

###  A.1 The homogenised Brusselator is not permanent

Let us now restrict the dynamics to $\mathcal{P}_{c}$. After elimination of
$z$ by the conservation law $x+y+z=c$, we obtain the following ODE system on
$\\{(x,y)\in\mathbb{R}^{2}_{\geq 0}\,:\,x+y\leq c\\}$:

$\displaystyle\begin{split}\dot{x}&=\kappa_{1}(c-x-y)-\kappa_{2}x-\kappa_{3}x+\kappa_{4}x^{2}y,\\\
\dot{y}&=\kappa_{3}x-\kappa_{4}x^{2}y.\end{split}$ (7)

From here on, we will focus on this system, which can be seen as describing,
in local coordinates, the dynamics of the original ODEs restricted to a
particular stoichiometric class, parameterised by $c$. Note that parameters of
the system are now the four original rate constants
$\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4}$, along with $c$. For the
purposes of bifurcation analysis we assume, however, that $c$ is fixed.

Observe that the Jacobian matrix of (7) at the corner equilibrium $(0,c)$,
denoted by $J_{0}$, equals

$\displaystyle\begin{pmatrix}-\kappa_{1}-\kappa_{2}-\kappa_{3}&-\kappa_{1}\\\
\kappa_{3}&0\end{pmatrix}.$

Since $\det J_{0}>0$, $\operatorname{tr}J_{0}<0$, and

$\displaystyle(\operatorname{tr}J_{0})^{2}-4\det
J_{0}=(\kappa_{1}+\kappa_{2}+\kappa_{3})^{2}-4\kappa_{1}\kappa_{3}=(\kappa_{1}+\kappa_{2}-\kappa_{3})^{2}+4\kappa_{2}\kappa_{3}>0,$

both eigenvalues are real and negative. Therefore, the boundary equilibrium
$(0,c)$ is asymptotically stable. In particular, for all
$(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4},c)\in\mathbb{R}^{5}_{+}$, (7) is
_not_ permanent, and since $c>0$ was arbitrary, the same clearly holds for the
homogenised Brusselator on each positive stoichiometric class. By contrast,
the Brusselator was permanent for all positive values of the rate constants.

###  A.2 Local bifurcation analysis

We proceed as follows. We first parameterise the set of equilibria which
simplifies many calculations. The parameter used, denoted by $t$, can replace
the parameter $c$ in many calculations: along the branch of equilibria we
consider, $t$ and $c$ are in one-to-one correpondence. We write down necessary
conditions for fold and Andronov-Hopf bifurcations to occur, and confirm that
both supercritical and subcritical Andronov-Hopf bifurcations can occur. We
also confirm that two codimension-$2$ bifurcations can occur on stoichiometric
classes as we vary the rate constants: a generalised Andronov-Hopf
bifurcation, also known as a Bautin bifurcation; and a Bogdanov-Takens
bifurcation. We can check that apart from along exceptional sets,
nondegeneracy and transversality conditions hold for all the bifurcations. See
[9] for full details of all the calculations.

We first note that the set of positive equilibria of the homogenised
Brusselator is

$\displaystyle
E_{+}=\left\\{\left(t,\frac{\kappa_{3}}{\kappa_{4}}\frac{1}{t},\frac{\kappa_{2}}{\kappa_{1}}t\right)\colon
t>0\right\\}.$

The parameters $c$ and $t$ are related via
$t\frac{\kappa_{1}+\kappa_{2}}{\kappa_{1}}+\frac{\kappa_{3}}{\kappa_{4}t}=c$.
Defining

$c^{*}:=2\sqrt{\frac{(\kappa_{1}+\kappa_{2})\kappa_{3}}{\kappa_{1}\kappa_{4}}},$

we find that $|E_{+}\cap\mathcal{P}_{c}|=0$ if $c<c^{*}$, while
$|E_{+}\cap\mathcal{P}_{c}|=1$ if $c=c^{*}$, and
$|E_{+}\cap\mathcal{P}_{c}|=2$ if $c>c^{*}$. The equilbrium set intersects the
stoichiometic class with parameter $c^{*}$ when $t=t^{*}$ where

$t^{*}:=\sqrt{\frac{\kappa_{1}\kappa_{3}}{(\kappa_{1}+\kappa_{2})\kappa_{4}}}.$

The Jacobian matrix of (7) along $E_{+}$, denoted by $J(t)$, equals

$\displaystyle\begin{pmatrix}\kappa_{3}-\kappa_{1}-\kappa_{2}&\kappa_{4}t^{2}-\kappa_{1}\\\
-\kappa_{3}&-\kappa_{4}t^{2}\end{pmatrix}.$

Its determinant and trace equal

$\det
J(t)=(\kappa_{1}+\kappa_{2})\kappa_{4}t^{2}-\kappa_{1}\kappa_{3}\quad\mbox{and}\quad\operatorname{tr}J(t)=-\kappa_{4}t^{2}+(\kappa_{3}-\kappa_{1}-\kappa_{2}).$

In the positive stoichiometric classes $\mathcal{P}_{c}$ with $2$ positive
equilibria, the equilibrium corresponding to the value of $t<t^{*}$, which is
closer to the corner equilibrium $(0,c)$, is a saddle for all values of the
rate constants. For the purposes of local bifurcation analysis, we thus focus
attention on equilibria satisfying $t\geq t^{*}$. As values of $t$ and $c$ are
in one-to-one correspondence for $t\geq t^{*}$ and $c\geq c^{*}$, we can thus
pass back and forth between parameters $c$ and $t$, and do so, sometimes
without comment, in calculations to follow.

Fold bifurcations of equilibria. These potentially occur when
$\mathrm{det}\,J(t)=0$, namely on the set $\mathrm{T}$ in parameter space
where $t=t^{*}$:

$\mathrm{T}:=\\{(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4},t)\in\mathbb{R}^{5}_{+}\colon
t^{2}=\frac{\kappa_{1}\kappa_{3}}{(\kappa_{1}+\kappa_{2})\kappa_{4}}\\}\,.$

Equivalently, fold bifurcations potentially occur when $c=c^{*}$, namely
along:

$\widetilde{\mathrm{T}}:=\\{(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4},c)\in\mathbb{R}^{5}_{+}\colon
c^{2}=4\frac{(\kappa_{1}+\kappa_{2})\kappa_{3}}{\kappa_{1}\kappa_{4}}\\}\,.$

For any fixed value of $c$, $\widetilde{\mathrm{T}}$ tells us combinations of
rate constants for which we expect a fold bifurcation to occur on the
stoichiometric class parameterised by $c$. It is straightforward to confirm
that, provided $\kappa_{3}\neq(\kappa_{1}+\kappa_{2})^{2}/\kappa_{2}$, the
fold bifurcation is nondegenerate, and is unfolded nondegenerately by the rate
constants $\kappa_{i}$ (any of $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$ or
$\kappa_{4}$ serves to unfold the bifurcation nondegenerately for all rate
constants). The degenerate case where
$\kappa_{3}=(\kappa_{1}+\kappa_{2})^{2}/\kappa_{2}$ will be of importance
later.

Andronov-Hopf bifurcations. These potentially occur when $\mathrm{tr}\,J(t)=0$
and $\mathrm{det}\,J(t)>0$, namely, when

$t^{2}=\frac{\kappa_{3}-\kappa_{1}-\kappa_{2}}{\kappa_{4}},\quad
t^{2}>\frac{\kappa_{1}\kappa_{3}}{(\kappa_{1}+\kappa_{2})\kappa_{4}}\,.$

Combining these conditions, Andronov-Hopf bifurcations potentially occur along

$\mathrm{H}:=\\{(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4},t)\in\mathbb{R}^{5}_{+}\colon
t^{2}=\frac{\kappa_{3}-\kappa_{1}-\kappa_{2}}{\kappa_{4}},\,\quad\kappa_{2}\kappa_{3}>(\kappa_{1}+\kappa_{2})^{2}\\}\,.$

(We may, if desired, write the bifurcation set in terms of the parameters
$\kappa_{i}$ and $c$ instead of $\kappa_{i}$ and $t$.) To check nondegeneracy
of the Andronov-Hopf bifurcations, we compute the first focal value, $L_{1}$,
along $\mathrm{H}$. Defining $a=\frac{\kappa_{2}}{\kappa_{1}}$ and
$b=\frac{\kappa_{3}}{\kappa_{1}}$, we obtain
$L_{1}=\frac{1}{t^{2}}\frac{P(a,b)}{Q(a,b)}$, where

$\displaystyle P(a,b)$
$\displaystyle=b^{3}[2-a]-b^{2}[5+3a-a^{2}]+b[5+12a+8a^{2}+a^{3}]-[4+13a+15a^{2}+7a^{3}+a^{4}],$
$\displaystyle Q(a,b)$ $\displaystyle=(ab-(1+a)^{2})^{\frac{3}{2}}(b-a-2).$

Note that the necessary condition for Andronov-Hopf bifurcation
$\kappa_{2}\kappa_{3}>(\kappa_{1}+\kappa_{2})^{2}$ is equivalent to
$ab>(1+a)^{2}$. Since $Q(a,b)$ is positive whenever this condition is
satisfied, it suffices to investigate the sign of $P(a,b)$. As $a$ and $b$
vary along $\mathrm{H}$, $P(a,b)$, and hence $L_{1}$, can be positive,
negative or zero along the bifurcation set, corresponding to subcritical,
supercritical, and degenerate Andronov-Hopf bifurcations respectively, as
shown in Figure 1.

Figure 1: Analysis of the Andronov-Hopf bifurcation for the mass-action system
(7) on $\mathrm{H}$. The figure shows the sign of the first focal value
$L_{1}$ as a function of $a=\frac{\kappa_{2}}{\kappa_{1}}$ and
$b=\frac{\kappa_{3}}{\kappa_{1}}$. The black curve is where $ab=(1+a)^{2}$,
namely the boundary of $\mathrm{H}$ where the determinant of the Jacobian
matrix vanishes and we expect Bogdanov-Takens bifurcations to occur; notice
that it touches only the $L_{1}>0$ region, and we can indeed confirm that all
Bogdanov-Takens bifurcations are subcritical. Along the red curve, where
$L_{1}=0$ we can confirm that $L_{2}>0$. This is where we expect Bautin
bifurcations to occur and, in fact, we find these bifurcations to be
nondegenerate away from an exceptional set.

We can also confirm that wherever $L_{1}\neq 0$, the parameters unfold the
bifurcation nondegenerately.

Bautin bifurcations. Generically, points along $\mathrm{H}$ where $L_{1}$
vanishes correspond to a generalised Andronov-Hopf bifurcation, also known as
a Bautin bifurcation [25]. We define the parameter set corresponding to
potential Bautin bifurcations as

$\mathrm{GH}:=\\{(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4},t)\subseteq\mathrm{H}\colon
L_{1}=0\\}\,.$

Indeed, the nondegeneracy conditions for the Bautin bifurcation can be checked
and found to hold on $\mathrm{GH}$ apart from along an exceptional set:

  1. 1.

Along $\mathrm{GH}$, we can compute the second focal value, $L_{2}$, and find
it is always positive. The calculations are lengthy and are omitted, but can
be found in the supporting documentation [9].

  2. 2.

The parameters $(\kappa_{3},\kappa_{4})$ unfold the bifurcation
nondegenerately. Choosing these as our bifurcation parameters, then for fixed
choices of the remaining parameters, we can confirm that the regularity
condition (B.2) in [25, Theorem 8.2] is almost always satisfied on
$\mathrm{GH}$, apart from along an exceptional set where
$\kappa_{2}/\kappa_{1}=\overline{a}=\sqrt{\frac{73+22\sqrt{11}}{20}}-1\approx
1.70$. On Figure 1, a vertical line at $a=\overline{a}$ would touch the
$L_{1}=0$ curve. The calculations are again omitted here but available in the
supporting documentation [9].

Thus, away from the exceptional set, the Bautin bifurcation is nondegenerate.
For parameter values near to $\mathrm{GH}$, on certain stoichiometric classes
an unstable positive equilibrium is surrounded by a pair of periodic orbits,
one stable and one unstable. As we vary rate constants, a fold bifurcation of
periodic orbits must occur on some stoichiometric class.

Bogdanov-Takens (B-T) bifurcations. These potentially occur where $\mathrm{T}$
meets the closure of $\mathrm{H}$, namely along the set where
$\mathrm{det}\,J(t)=0$ and $\mathrm{tr}\,J(t)=0$, defined by

$\mathrm{BT}:=\\{(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4},t)\in\mathbb{R}^{5}_{+}\colon\kappa_{2}\kappa_{3}=(\kappa_{1}+\kappa_{2})^{2},\quad
t^{2}=\frac{\kappa_{1}\kappa_{3}}{(\kappa_{1}+\kappa_{2})\kappa_{4}}\\}\,.$

We can confirm that the B-T bifurcation is always nondegerate, subcritical,
and unfolded nondegenerately by parameters $(\kappa_{3},\kappa_{4})$. This
involves checking conditions (BT.0), (BT.1), (BT.2), and (BT.3) in [25,
Theorem 8.4]. The computations can be carried out explicitly, but are lengthy,
and so are omitted here. They can be found in the supporting documentation
[9]. Thus for parameter values close to $\mathrm{BT}$, on some stoichiometric
classes, the system can have a stable equilibrium surrounded by an unstable
periodic orbit. Furthermore, homoclinic bifurcations can occur.

In Figure 2, we fix $\kappa_{1}=2$, $\kappa_{2}=4$, and $c=6$, and depict the
curves corresponding to $\mathrm{T}$ and $\mathrm{H}$, and the points
corresponding to $\mathrm{BT}$ and $\mathrm{GH}$, in the
$(\kappa_{3},\kappa_{4})$ plane. This completes our analysis of the local
bifurcations of equilibria in the homogeneous Brusselator.

Figure 2: Bifurcation diagram in the $(\kappa_{3},\kappa_{4})$-plane of the
mass-action system (7) with $\kappa_{1}$, $\kappa_{2}$ and $c$ fixed at
$\kappa_{1}=2$, $\kappa_{2}=4$, $c=6$. On the solid ochre curve (T) a fold
bifurcation of equilibria occurs. Along the dashed blue curve (H) there is an
Andronov-Hopf bifurcation. A Bogdanov-Takens bifurcation occurs at (BT), where
(H) touches (T). On the curve (H), the bifurcation is subcritical close to the
(BT) point, but becomes supercritical beyond the point (GH), where the first
focal value $L_{1}$ vanishes. At (GH), there is a Bautin bifurcation.

## References

## References

  * [1] D. Angeli, M. Banaji, and C. Pantea.  Combinatorial approaches to Hopf bifurcations in systems of interacting elements.  Commun. Math. Sci., 12:1101–1133, 2014.
  * [2] D. Angeli, P. De Leenheer, and E. D. Sontag.  Graph-theoretic characterizations of monotonicity of chemical reaction networks in reaction coordinates.  J. Math. Biol., 61(4):581–616, 2010.
  * [3] M. Banaji.  Inheritance of oscillation in chemical reaction networks.  Appl. Math. Comput., 325:191–209, 2018.
  * [4] M. Banaji.  Building oscillatory chemical reaction networks by adding reversible reactions.  SIAM J. Appl. Math., 80(4):1751–1777, 2020.
  * [5] M. Banaji.  Splitting reactions preserves nondegenerate behaviours in chemical reaction networks, 2022.  https://arxiv.org/abs/2201.13105.pdf.
  * [6] M. Banaji and G. Craciun.  Graph-theoretic approaches to injectivity and multiple equilibria in systems of interacting elements.  Commun. Math. Sci., 7(4):867–900, 2009.
  * [7] M. Banaji and C. Pantea.  Some results on injectivity and multistationarity in chemical reaction networks.  SIAM J. Appl. Dyn. Syst., 15(2):807–869, 2016.
  * [8] M. Banaji and C. Pantea.  The inheritance of nondegenerate multistationarity in chemical reaction networks.  SIAM J. Appl. Math., 78(2):1105–1130, 2018.
  * [9] B. Boros.  Reaction Networks GitHub repository, 2022.  https://github.com/balazsboros/reaction_networks.
  * [10] B. Boros and J. Hofbauer.  Limit cycles in mass-conserving deficiency-one mass-action systems, 2022\.  https://arxiv.org/abs/2202.10406.pdf.
  * [11] B. Boros and J. Hofbauer.  Some minimal bimolecular mass-action systems with limit cycles, 2022.  https://arxiv.org/abs/2202.11034.pdf.
  * [12] R. Cheong and A. Levchenko.  Oscillatory signaling processes: the how, the why and the where.  Curr. Opin. Genet. Dev., 20(6):665–669, 2010.
  * [13] G. Craciun and M. Feinberg.  Multiple equilibria in complex chemical reaction networks: I. The injectivity property.  SIAM J. Appl. Math., 65(5):1526–1546, 2005.
  * [14] P. Érdi and J. Tóth.  Mathematical models of chemical reactions. Theory and applications of deterministic and stochastic models.  Princeton University Press, Princeton, 1989.
  * [15] H. Farkas and Z. Noszticzius.  Generalized Lotka-Volterra schemes and the construction of two-dimensional explodator cores and their Liapunov functions via ‘critical’ Hopf bifurcations.  J. Chem. Soc., Faraday Trans. 2, 81(10):1487–1505, 1985.
  * [16] M. Feinberg.  Chemical reaction network structure and the stability of complex isothermal reactors - I. The deficiency zero and deficiency one theorems.  Chem. Eng. Sci., 42(10):2229–2268, 1987.
  * [17] E. Feliu and C. Wiuf.  Preclusion of switch behavior in networks with mass-action kinetics.  Appl. Math. Comput., 219:1449–1467, 2012.
  * [18] E. Feliu and C. Wiuf.  Simplifying biochemical models with intermediate species.  J. Roy. Soc. Interface, 10:20130484, 2013.
  * [19] N. Fenichel.  Persistence and smoothness of invariant manifolds for flows.  Indiana Univ. Math. J., 21(3):193–226, 1971.
  * [20] N. Fenichel.  Geometric singular perturbation theory for ordinary differential equations.  J. Differ. Equ., 31(1):53–98, 1979.
  * [21] M. W. Hirsch, C. C. Pugh, and M. Shub.  Invariant Manifolds, volume 583 of Lecture Notes in Mathematics.  Springer-Verlag, 1977.
  * [22] F. Horn and R. Jackson.  General mass action kinetics.  Arch. Ration. Mech. Anal., 47(2):81–116, 1972.
  * [23] B. Joshi and A. Shiu.  Atoms of multistationarity in chemical reaction networks.  J. Math. Chem., 51(1):153–178, 2013.
  * [24] B. N. Kholodenko.  Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades.  Eur. J. Biochem., 267(6):1583–1588, 2000.
  * [25] Y. A. Kuznetsov.  Elements of Applied Bifurcation Theory, volume 112 of Applied Mathematical Sciences.  Springer-Verlag, New York, third edition, 2004.
  * [26] M. Y. Li.  Dulac criteria for autonomous systems having an invariant affine manifold.  J. Math. Anal. Appl., 199(2):374–390, 1996.
  * [27] K. Lin, X. Tang, and Z. Zhang.  Multistationarity of reaction networks with one-dimensional stoichiometric subspaces, 2021.  https://arxiv.org/abs/2108.09695.
  * [28] N. I. Markevich, J. B. Hoek, and B. N. Kholodenko.  Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades.  J. Cell Biol., 164(3):353–359, 2004.
  * [29] B. Novák and J. J. Tyson.  Design principles of biochemical oscillators.  Nature Rev. Mol. Cell Biol., 9:981–991, 2008.
  * [30] N. Obatake, A. Shiu, X. Tang, and A. Torres.  Oscillations and bistability in a model of ERK regulation.  J. Math. Biol., 79:1515–1549, 2019.
  * [31] F. Ortega, J. L. Garcés, F. Mas, B. N. Kholodenko, and M. Cascante.  Bistability from double phosphorylation in signal transduction.  FEBS Journal, 273(17):3915–3926, 2006.
  * [32] C. Pantea.  On the persistence and global stability of mass-action systems.  SIAM J. Math. Anal., 44(3):1636–1673, 2012.
  * [33] L. Qiao, R. B. Nachbar, I. G. Kevrekidis, and S. Y. Shvartsman.  Bistability and oscillations in the Huang-Ferrell model of MAPK signaling.  PLoS Comput. Biol., pages 1819–1826, 2007.
  * [34] F. Schlögl.  On thermodynamics near a steady state.  Z. Physik, 248(5):446–458, 1971.
  * [35] F. Schlögl.  Chemical reaction models for non-equilibrium phase transitions.  Z. Physik, 253(2):147–161, 1972.
  * [36] K. R. Schneider.  Über die periodischen Lösungen einer Klasse nichtlinearer autonomer Differentialgleichungssysteme dritter Ordnung.  Z. Angew. Math. Mech., 49(7):441–443, 1969.
  * [37] G. Shinar and M. Feinberg.  Concordant chemical reaction networks.  Math. Biosci., 240(2):92–113, 2012.
  * [38] P. L. Simon.  The reversible LVA model.  J. Math. Chem., 9(4):307–322, 1992.
  * [39] H. Smith.  Global dynamics of the smallest chemical reaction system with Hopf bifurcation.  J. Math. Chem., 50(4):989–995, 2012.
  * [40] J. J. Tyson, R. Albert, A. Goldbeter, P. Ruoff, and J. Sible.  Biological switches and clocks.  J. R. Soc. Interface, 5:S1–S8, 2008.
  * [41] S. Wiggins.  Introduction to Applied Nonlinear Dynamics and Chaos.  Springer, 2003.
  * [42] Y.-Q. Ye et al.  Theory of Limit Cycles.  Translations of Mathematical Monographs, Vol. 66, American Mathematical Society, Providence, R. I., 1986.