Inheritance of oscillation in chemical reaction networks

TL;DR

This paper investigates how oscillatory behavior in chemical reaction networks (CRNs) can be inherited when networks are enlarged, providing conditions under which oscillations are preserved, especially in fully open CRNs with mass action kinetics.

Contribution

It establishes four key methods for enlarging CRNs while maintaining oscillation capacity and extends results to various kinetics, highlighting the importance of small oscillatory motifs.

Findings

Enlargement methods preserve oscillation in CRNs.

Nondegenerate and stable oscillations are inherited in fully open CRNs.

Small oscillatory motifs are significant sources of oscillation.

Abstract

Some results are presented on how oscillation is inherited by chemical reaction networks (CRNs) when they are built in natural ways from smaller oscillatory networks. The main results describe four important ways in which a CRN can be enlarged while preserving its capacity for oscillation. The results are for general CRNs, not necessarily fully open, but lead to an important corollary for fully open networks: if a fully open CRN R with mass action kinetics admits a nondegenerate (resp., linearly stable) periodic orbit, then so do all such CRNs which include R as an induced subnetwork. This claim holds for other classes of kinetics, but fails, in general, for CRNs which are not fully open. Where analogous results for multistationarity can be proved using the implicit function theorem alone, the results here call on regular and singular perturbation theory. Equipped with these results and with the help of some analysis and numerical simulation, lower bounds are put on the proportion of small fully open CRNs capable of stable oscillation under various assumptions on the kinetics. This exploration suggests that small oscillatory motifs are an important source of oscillation in CRNs.

Tables

Table 1. Table 2: Simulations were carried out on the 444 ( 3 , 3 ) 3 3 (3,3) CRNs known to admit SPPOs with mass action kinetics to determine the effect of number of simulations on the likelihood of observing SPPOs in numerical experiment. For example, even when 100,000 simulations were carried out on each CRN, only 298 (67%) of the CRNs were identified as oscillatory.

parameter sets	CRNs found to admit SPPOs
100	17
1,000	63
10,000	174
100,000	298

Equations

$$\dot{x}=F(x)$$

$$\dot{z}=DF(\theta (t))z.$$

$$Z(t)=A(t)e^{tB}$$

$$\dot{x}=F(x,ϵ).$$

$$\dot{x}=\Gamma v(x).$$

$$\dot{y}=Qv(x_0+{\Gamma }_{0}y).$$

$$\dot{z}=QDv(\theta (t)){\Gamma }_{0}z$$

$$\dot{y}=R^{-1}Qv(x_0^{\prime }+{\Gamma }_0^{\prime }y)$$

$$\dot{z}=R^{-1}[QDv(\theta (t)){\Gamma }_0]Rz.$$

$$X+Y\to 2Y,Y+Z\to X\to W+Z,W\to X.$$

$$X+Y\to 2Y,Y+Z\to X,$$

$$\dot{x}=\Gamma v(x).$$

$$({\mathcal{R}}_1,{\mathcal{K}}_1){\le }_S({\mathcal{R}}^{\prime },{\mathcal{K}}^{\prime }){\le }_R({\mathcal{R}}_2,{\mathcal{K}}_2).$$

$$\dot{x}=\Gamma v(x)+ϵ\alpha f(x)={\Gamma }_0(Qv(x)+ϵcf(x)).$$

$$\dot{z}=Qv(x_0+{\Gamma }_{0}z)+ϵcf(x_0+{\Gamma }_{0}z).$$

$$\dot{x}=\Gamma v(x)+ϵI_n(x_0-x).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{cc}\hat{z}\\ \hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{z}}}\end{array}\right)=\left(\begin{array}[]{c}Q\\ 0\end{array}\right)v(x_{0}+\Gamma_{0}^{\prime}z)-\epsilon\left(\begin{array}[]{cc}\hat{z}\\ \hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{z}}}\end{array}\right)\,.$$

$$\frac{\mathrm{d}\hat{z}}{\mathrm{d}t}=Qv(x_0+{\Gamma }_0\hat{z})-ϵ\hat{z}.$$

$$\frac{\mathrm{d}\hat{\zeta }}{\mathrm{d}t}=[QDv(x_0+{\Gamma }_0{\psi }_{ϵ}(t)){\Gamma }_0-ϵI_r]\hat{\zeta }.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{cc}\hat{\zeta}\\ \hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\zeta}}}\end{array}\right)=\left(\begin{array}[]{cc}QDv(x_{0}+\Gamma_{0}\psi_{\epsilon}(t))\Gamma_{0}-\epsilon I_{r}&QDv(x_{0}+\Gamma_{0}\psi_{\epsilon}(t))\Gamma_{1}\\ 0&-\epsilon I_{n-r}\end{array}\right)\left(\begin{array}[]{cc}\hat{\zeta}\\ \hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\zeta}}}\end{array}\right)\,.$$

$$\frac{\mathrm{d}\hat{\zeta }}{\mathrm{d}t}=[QDv(x_0+{\Gamma }_0{\psi }_{ϵ}(t)){\Gamma }_0-ϵI_r]\hat{\zeta }+e^{-ϵt}QDv(x_0+{\Gamma }_0{\psi }_{ϵ}(t)){\Gamma }_1\hat{\hat{\zeta }}(0).$$

$$Z_{\epsilon}(t)=\left(\begin{array}[]{cc}\hat{Z}_{\epsilon}(t)&A(t)\\ 0&e^{-\epsilon t}I_{n-r}\end{array}\right),\quad\mbox{and hence,}\quad Z_{\epsilon}(T_{\epsilon})=\left(\begin{array}[]{cc}\hat{Z}_{\epsilon}(T_{\epsilon})&A(T_{\epsilon})\\ 0&e^{-\epsilon T_{\epsilon}}I_{n-r}\end{array}\right).$$

$$\left(\begin{array}[]{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}[]{c}\Gamma\\ 0\end{array}\right)w(x,y).$$

$$\dot{z}=Qv(x_0+{\Gamma }_{0}z).$$

$$\begin{array}[]{rcl}\dot{x}&=&\Gamma w(x,y)\\ \epsilon\dot{y}&=&\epsilon sw(x,y)+(1-f(y)).\end{array}$$

$$\begin{array}[]{rcl}\dot{x}&=&\epsilon\Gamma w(x,y)\\ \dot{y}&=&\epsilon sw(x,y)+(1-f(y)).\end{array}$$

$$\begin{array}[]{rcl}\dot{z}&=&Qw(x_{0}+\Gamma_{0}z,y)\\ \epsilon\dot{y}&=&\epsilon sw(x_{0}+\Gamma_{0}z,y)+(1-f(y)).\end{array}$$

$$\begin{array}[]{rcl}\dot{z}&=&\epsilon Qw(x_{0}+\Gamma_{0}z,y)\\ \dot{y}&=&\epsilon sw(x_{0}+\Gamma_{0}z,y)+(1-f(y)).\end{array}$$

$$\dot{z}=Qv(x_0+{\Gamma }_{0}z)=:H(z),y=1,$$

Full text

1 Introduction and context of the paper

At the heart of many biological systems are chemical reaction networks (CRNs), and the question of when these admit oscillation is of both theoretical and practical interest. Oscillation is known to occur – and play a key role – in a great variety of biological contexts. Examples include the natural rhythms of body clocks and ovulation, biochemical oscillations in cellular signalling, cyclic behaviour of various diseases, and periodic fluctuations in Lotka-Volterra-type models of interacting populations. Several chapters of [1] and [2] detail mathematical models of oscillation in biological settings. Some general biological principles underlying biological oscillation are discussed in [3]. Once a network admitting oscillation is identified, we might naturally wonder whether this network occurs as a “motif” in other larger networks and, if so, whether the larger networks must themselves admit oscillation. The desire to phrase this question precisely and provide some simple and partial answers motivates this work.

Several papers have treated analogous questions about the inheritance of multistationarity in CRNs [4, 5, 6, 7]. In a recent contribution it was shown that a great deal can be done in this direction using the implicit function theorem [8]. An (incomplete) list of network modifications proven to preserve the property of admitting nondegenerate multistationarity were listed; these collectively define a partial order $\preceq$ on the set of all CRNs such that if a CRN $\mathcal{R}$ admits nondegenerate multistationarity, then so do all CRNs $\succeq\mathcal{R}$ in this partial order. Although it is likely that most, if not all, of the results in [8] can be restated with “nondegenerate oscillation” replacing “nondegenerate multistationarity”, only part of this task is undertaken here: we prove four results about general CRNs, Theorems 1 to 4, which are analogues of related results about multistationarity in [8], also numbered Theorems 1 to 4. An example of what these tell us is the following corollary about fully open CRNs:

Proposition 1.1**.**

If a fully open CRN $\mathcal{R}$ with mass action kinetics admits nondegenerate (resp., stable) oscillation, then so does any fully open CRN with mass action kinetics which includes $\mathcal{R}$ as an induced subnetwork.

The definitions required to make this result precise will follow. Proposition 1.1 is the specialisation for mass action kinetics of a result with more general kinetic assumptions, Proposition 4.8, (see Remark 4.9) which is a natural starting point for some computational exploration on small fully open CRNs admitting oscillation. It is worth noting at the outset that Proposition 1.1 fails if the CRNs are not assumed to be fully open. An example is provided in the concluding section (Example 6.1).

Much of the mathematical literature on oscillation in CRNs has focussed on conditions which forbid oscillation, or forbid stable oscillation of the kind which might be observed in numerical simulations, or forbid bifurcations leading to oscillation. For CRNs with mass action kinetics, there are the original results of deficiency theory [9, 10, 11, 12]; for CRNs with more general kinetics there are results based on the theory of monotone dynamical systems ([13, 14, 15, 16] for example), and algebraic approaches ([17] for example). Various papers which do not directly treat CRNs also have natural applications to forbidding oscillation or stable oscillation in CRNs, including the work of Angeli, Hirsch and Sontag on “coherent” systems [18], and of Li and Muldowney on generalised Bendixson’s criteria [19, 20, 21]. On the other hand oscillation has been shown to occur in numerical studies of various CRNs of interest (for example, [22, 23, 24, 25]). Aside from numerical work, there exists an important strand of theory drawing on approaches in convex and toric geometry which provides sufficient conditions for Hopf bifurcations in CRNs with mass action and generalised mass action kinetics [26, 27, 28, 29]. These approaches lead to algorithms for the determination of parameter regions where Hopf bifurcation occurs. Other papers treating the question of sufficient conditions for oscillation in chemical reaction networks include [30] and [31].

The work here is aimed at closing the gap between theory which forbids oscillation and examples of oscillatory networks or particular sufficient conditions for oscillation. It is likely that many examples of CRNs admitting oscillation in fact oscillate because they inherit this property from a smaller CRN which admits oscillation, and the goal is then to identify an appropriate notion of inheritance, and minimal oscillatory CRNs in some sense. The importance of inheritance approaches is increasingly recognised. In [32], Conradi and Shiu pose a question closely related to the main question in this paper, namely whether Hopf bifurcation is preserved when CRNs are modified in natural ways. The problem of identifying a “minimal” oscillatory subnetwork was tackled for the biologically important MAPK cascade in [33].

Computational work on fully open CRNs towards the end of the paper confirms the practical usefulness of inheritance approaches. As oscillation may occur in very small regions of parameter space, it may be hard to find by brute force in numerical simulations, even where it is straightforward to predict its occurrence by inheritance results. Finding a single small oscillatory CRN on the other hand immediately gives us knowledge of a large number of CRNs which inherit this oscillation. Ultimately, the hope is that examining CRNs which can neither be proven to forbid oscillation nor be shown to oscillate (using numerics, known sufficient conditions for oscillation, or inheritance results such as here) may lead to new theorems about necessary conditions for oscillation.

1.1 Notational preliminaries

Notation 1.2** (Nonnegative and positive vectors).**

A real vector $x=(x_{1},\ldots,x_{n})^{\mathrm{t}}$ is nonnegative (resp., positive) if $x_{i}\geq 0$ (resp., $x_{i}>0$) for each $i$, and we refer to the nonnegative (resp., positive) orthant in $\mathbb{R}^{n}$ as $\mathbb{R}^{n}_{\geq 0}$ (resp., $\mathbb{R}^{n}_{\gg 0}$). Subsets of $\mathbb{R}^{n}_{\gg 0}$ are referred to as positive.

Notation 1.3** (Vector of ones).**

$\mathbf{1}$* denotes a vector of ones whose length is inferred from the context.*

Notation 1.4** (Identity matrix).**

$I_{n}$* is the $n\times n$ identity matrix.*

Notation 1.5** (Set theoretic inverse).**

Given sets $X,Y$ and a function $f\colon X\to Y$, not necessarily invertible, $f^{-1}$ will generally refer to the set theoretic inverse, namely, given $Y_{0}\subseteq Y$, $f^{-1}(Y_{0})=\{x\in X\colon f(x)\in Y_{0}\}$.

Notation 1.6** (Monomials, vector of monomials).**

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 $\prod_{i}x_{i}^{a_{i}}$. If $A$ is an $m\times n$ matrix with rows $A_{1},\ldots,A_{m}$, then $x^{A}$ means the vector of (generalised) monomials $(x^{A_{1}},x^{A_{2}},\ldots,x^{A_{m}})^{\mathrm{t}}$.

Notation 1.7** (Entrywise product).**

Given two matrices $A$ and $B$ with the same dimensions, $A\circ B$ will refer to the entrywise (or Hadamard) product of $A$ and $B$, namely $(A\circ B)_{ij}=A_{ij}B_{ij}$.

2 Periodic orbits

We remind the reader of some standard results from Floquet theory (Chapters 3 and 4 of [34] for example) as needed here. Let $X\subseteq\mathbb{R}^{r}$ be open, $F\colon X\to\mathbb{R}^{r}$ be $C^{1}$, and consider the ODE

[TABLE]

on $X$. Assume that (2.1) has a nontrivial periodic solution $\theta\colon\mathbb{R}\to X$ with smallest positive period $T$, and with corresponding periodic orbit $\mathcal{O}:=\mathrm{im}\,\theta$. The variational equation about $\theta$ is

[TABLE]

$DF(\theta(t))$ is an $r\times r$ $T$-periodic matrix and Floquet theory tells us that any fundamental matrix solution $Z(t)$ of (2.2) can be written in the form

[TABLE]

where $A$ is a nonsingular $T$-periodic matrix, and $B$ is a constant matrix. The eigenvalues of $e^{TB}$ are termed the characteristic multipliers (or Floquet multipliers) of $\mathcal{O}$. If $Z(0)=I$, then $A(T)=A(0)=I$, in which case the characteristic multipliers are the eigenvalues of $Z(T)$. $\mathcal{O}$ is termed hyperbolic (resp., linearly stable) if $r-1$ of its characteristic multipliers are disjoint from (resp., inside) the unit circle in $\mathbb{C}$. Hyperbolicity (resp., linear stability) of a periodic orbit is precisely hyperbolicity (resp., linear stability) of the associated fixed point of any Poincaré map constructed on a section transverse to the periodic orbit: see Chapter 10 onwards of [35], for example. Hyperbolic periodic orbits survive under sufficiently small perturbations of vector fields in a sense made precise in Lemma 2.1 below. Linear stability of a periodic orbit implies asymptotic orbital stability, namely that forward trajectories of all sufficiently nearby initial conditions converge to the periodic orbit (Theorem 4.2 in [34]).

The following is a well-known result of regular perturbation theory. $d_{\mathrm{H}}(\cdot,\cdot)$ denotes the Hausdorff distance between nonempty compact subsets of Euclidean space.

Lemma 2.1**.**

Let $X\subseteq\mathbb{R}^{r}$ be open, $\epsilon^{\prime}>0$ and $F\colon X\times(-\epsilon^{\prime},\epsilon^{\prime})\to\mathbb{R}^{r}$ be $C^{1}$. Consider the $\epsilon$-dependent family of ODEs on $X$

[TABLE]

Suppose that (2.3) has a nontrivial hyperbolic (resp., linearly stable) $T$-periodic orbit $\mathcal{O}\subseteq X$. Then there exists $\epsilon_{0}>0$ s.t. for $\epsilon\in(-\epsilon_{0},\epsilon_{0})$ (2.3) has a hyperbolic (resp., linearly stable) periodic orbit $\mathcal{O}_{\epsilon}$ satisfying $\lim_{\epsilon\to 0}d_{\mathrm{H}}(\mathcal{O}_{\epsilon},\mathcal{O})=0$ and with period $T_{\epsilon}$ satisfying $\lim_{\epsilon\to 0}T_{\epsilon}=T$.

Proof 1*.*

These claims are proved, for example, by constructing a family of Poincaré maps $\Pi_{\epsilon}$ for (2.3) and applying the implicit function theorem at the fixed point of $\Pi_{0}$ corresponding to $\mathcal{O}$ as described in Section IV of [36]. $\square$

We next consider some specialisations of Floquet theory to systems with linear first integrals relevant to the study of CRNs. Let $x\in\mathbb{R}^{n}_{\gg 0}$, $v\colon\mathbb{R}^{n}_{\gg 0}\to\mathbb{R}^{m}$ be $C^{1}$, $\Gamma$ be an $n\times m$ real matrix of rank $r$, and consider the ODE

[TABLE]

Assume that (2.4) has a nontrivial positive periodic orbit $\mathcal{O}$ (see Notation 1.2), namely there exists some periodic solution $\theta\colon\mathbb{R}\to\mathbb{R}^{n}_{\gg 0}$ of (2.4) with smallest period $T>0$ and with $\mathcal{O}:=\mathrm{im}\,\theta$. Clearly, $S_{\mathcal{O}}:=(\mathcal{O}+\mathrm{im}\,\Gamma)\cap\mathbb{R}^{n}_{\gg 0}$ is locally invariant under (2.4). If $r\neq n$, then $\mathcal{O}$ cannot be hyperbolic or linearly stable in the senses defined above. However, our interest is in whether it is hyperbolic (resp., linearly stable) relative to $S_{\mathcal{O}}$. Associated with $S_{\mathcal{O}}$ are $r$ characteristic multipliers and we would like to know whether $r-1$ of these are disjoint from (resp., inside) the unit circle. The single remaining multiplier associated with $S_{\mathcal{O}}$, corresponding to travel along the periodic orbit, is $1$, while the additional $n-r$ multipliers associated with directions transverse to $S_{\mathcal{O}}$ are also easily shown all to be $1$.

An explicit calculation of the multipliers of $\mathcal{O}$ relative to $S_{\mathcal{O}}$ is needed in certain proofs to follow. This proceeds as follows. Choose $x_{0}\in S_{\mathcal{O}}$ and choose $\Gamma_{0}$ to be any matrix whose columns form a basis for $\mathrm{im}\,\Gamma$. Define $Q$ by $\Gamma=\Gamma_{0}Q$, and define the bijection $h\colon\mathbb{R}^{r}\to x_{0}+\mathrm{im}\,\Gamma$ by $h(y)=x_{0}+\Gamma_{0}y$. Note that $W:=h^{-1}(S_{\mathcal{O}})$ is an open subset of $\mathbb{R}^{r}$, and $\left.h\right|_{W}$ is an affine bijection between $W$ and $S_{\mathcal{O}}$. Setting $x=h(y)$ we get, for the evolution of $y$:

[TABLE]

(2.5) has a $T$-periodic solution $\psi\colon\mathbb{R}\to W$ defined by $\psi(t)=h^{-1}(\theta(t))$. Let $\mathcal{O}^{\prime}:=\mathrm{im}\,\psi=h^{-1}(\mathcal{O})$ be the corresponding periodic orbit. The multipliers of $\mathcal{O}^{\prime}$ are precisely the multipliers of $\mathcal{O}$ relative to $S_{\mathcal{O}}$. By definition $\mathcal{O}^{\prime}$ is hyperbolic (resp., linearly stable) if it has $r-1$ characteristic multipliers disjoint from (resp., inside) the unit circle in $\mathbb{C}$. This motivates the following definitions:

Definition 2.2** (NPPO, SPPO).**

Let $\mathcal{O}$ be a positive periodic orbit of (2.4). With $h$ defined as above, $\mathcal{O}$ is a nondegenerate positive periodic orbit (NPPO) of (2.4) if $\mathcal{O}^{\prime}:=h^{-1}(\mathcal{O})$ is a hyperbolic periodic orbit of (2.5). $\mathcal{O}$ is a linearly stable positive periodic orbit (SPPO) of (2.4) if $\mathcal{O}^{\prime}:=h^{-1}(\mathcal{O})$ is a linearly stable periodic orbit of (2.5). An SPPO is clearly also an NPPO.

The use of the transformation $h$ to define new coordinates on $S_{\mathcal{O}}$ is illustrated schematically in Figure 1.

Remark 2.3*.*

Note that the overloading of the term “linearly stable” in Definition 2.2 is an abuse of terminology which should cause no confusion: if $\mathrm{im}\,\Gamma=\mathbb{R}^{n}$, then linearly stable has its usual meaning; if $\mathrm{im}\,\Gamma\neq\mathbb{R}^{n}$, then no periodic orbit of (2.4) can truly be linearly stable, and linear stability is taken to mean linear stability relative $\mathrm{im}\,\Gamma$.

We can easily verify that Definition 2.2 makes sense: different choices of $x_{0}$ or $\Gamma_{0}$ lead to the same characteristic multipliers. To see this, recall that according to Floquet theory the variational equation of (2.5) about $\psi(t)=h^{-1}(\theta(t))$, namely,

[TABLE]

has a fundamental matrix solution $Z(t)$ which can be written $Z(t)=A(t)e^{tB}$ with $A$ a nonsingular $T$-periodic matrix and $B$ a constant matrix. The characteristic multipliers associated with $\psi$ are the eigenvalues of $e^{TB}$. Now suppose we make some different choices $x_{0}^{\prime}\in S_{\mathcal{O}}$ and $\Gamma_{0}^{\prime}$ and let $h^{\prime}\colon\mathbb{R}^{r}\to S_{\mathcal{O}}$ be defined by $h^{\prime}(y):=x_{0}^{\prime}+\Gamma_{0}^{\prime}y$. As the columns of $\Gamma_{0}^{\prime}$ are a basis for $\mathrm{im}\,\Gamma$, $\Gamma_{0}^{\prime}=\Gamma_{0}R$ where $R$ is a nonsingular $r\times r$ matrix. Thus $\Gamma=\Gamma_{0}Q=\Gamma_{0}^{\prime}R^{-1}Q$. With $x=h^{\prime}(y)$, we get the evolution on $W^{\prime}:=h^{\prime-1}(S_{\mathcal{O}})$

[TABLE]

with $T$-periodic solution $\psi^{\prime}\colon\mathbb{R}\to W^{\prime}$ defined by $\psi^{\prime}(t)=h^{\prime-1}(\theta(t))$. The variational equation of (2.7) about $\psi^{\prime}(t)$ is

[TABLE]

Then $Z^{\prime}(t):=R^{-1}Z(t)=R^{-1}A(t)e^{tB}$ is a fundamental matrix solution of (2.8) with $R^{-1}A(t)$ clearly a $T$-periodic matrix. Thus the characteristic multipliers associated with the solution $\psi^{\prime}(t)$ of (2.7) are again the eigenvalues of $e^{TB}$, i.e., those associated with the solution $\psi(t)$ of (2.5).

3 Background on CRNs

As the framework and terminology closely follow that of [37], the reader is referred to this paper for some of the detail. The goal is to remain precise while minimising the extensive preamble on basic notation, terminology and definitions which accompanies many papers on CRNs. We consider a CRN involving $n$ chemical species $X_{1},\ldots,X_{n}$.

Definition 3.1** (Complexes, the zero complex, stoichiometry).**

A complex is a formal linear combination of species. If $a=(a_{1},\ldots,a_{n})^{\mathrm{t}}$ is a nonnegative integer vector, then $a\cdot X:=a_{1}X_{1}+a_{2}X_{2}+\cdots+a_{n}X_{n}$ is a complex. $a_{i}$ is the stoichiometry of $X_{i}$ in the complex $a\cdot X$. The zero complex $0X_{1}+\cdots+0X_{n}$ is denoted [math].

An irreversible reaction is an ordered pair of complexes, termed the source complex (or left hand side) and the target complex (or right hand side). We always assume that the source and target complexes are distinct. A reversible reaction may be considered either as two irreversible reactions or, equivalently, as an unordered pair of (distinct) complexes. A CRN is a set of species and a set of reactions. We adopt the common convention that the reactions of a CRN are distinct. However, for technical reasons, we do not forbid a priori the possibility that some chemical species occurs in a CRN but participates in none of its reactions.

Definition 3.2** (Flow reaction, fully open CRN, fully open extension of a CRN).**

For the purposes of this paper, reactions of the form $0\rightarrow A$ or $A\rightarrow 0$ are referred to as flow reactions, while all others are non-flow reactions (even where these clearly violate any conservation laws: for example $2A\rightarrow 0$ or $0\rightarrow A+B$ are referred to as non-flow reactions). A CRN involving species $X_{1},\ldots,X_{n}$ is fully open if it includes all flow reactions $0\rightleftharpoons X_{i}$ ($i=1,\ldots,n$). Note that if, for example, a CRN includes all reactions of the form $0\rightleftharpoons 2X_{i}$, but not all the reactions $0\rightleftharpoons X_{i}$, we do not refer to it as fully open. The fully open extension of a CRN $\mathcal{R}$ is the smallest fully open CRN containing all the reactions of $\mathcal{R}$, namely the CRN created by adjoining to $\mathcal{R}$ any flow reactions which are absent from $\mathcal{R}$.

3.1 Combinatorial representations of CRNs

CRNs are combinatorial objects which give rise to dynamical systems in different ways depending on various modelling choices. The most common combinatorial representation of a CRN is via its complex graph [38], a digraph whose vertices are complexes and whose arcs correspond to (irreversible) reactions. For example, the reaction $X_{1}+2X_{2}\rightarrow X_{3}$ is an ordered pair of complexes naturally represented as an arc from source complex $X_{1}+2X_{2}$ to target complex $X_{3}$. The set of species and the complex graph together make up a formal description of the CRN.

An alternative representation, particularly useful when discussing isomorphism of CRNs, is a Petri net (PN) graph [39], an edge-weighted bipartite digraph, defined in the form used here in [8]. The PN graph of a CRN $\mathcal{R}$, denoted $PN(\mathcal{R})$, has two vertex sets $V_{S}$ (species vertices) and $V_{R}$ (reaction vertices) identified with the species and the reactions of $\mathcal{R}$. Given $X_{i}\in V_{S}$ and $R_{j}\in V_{R}$, there exists an arc $X_{i}R_{j}$ (resp., $R_{j}X_{i}$) with weight $w$ if and only if the species corresponding to $X_{i}$ occurs with stoichiometry $w>0$ in the source complex (resp., target complex) of the reaction corresponding to $R_{j}$. Arc weights of $1$ are omitted from drawings for neatness. An unlabelled PN graph is referred to as a motif.

CRNs $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are isomorphic if $PN(\mathcal{R}_{1})$ and $PN(\mathcal{R}_{2})$ are isomorphic in a natural sense, namely there exists a relabelling of the vertices of $PN(\mathcal{R}_{1})$ which preserves the bipartition and gives $PN(\mathcal{R}_{2})$. Given CRNs $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$, we say that $\mathcal{R}_{1}$ is an induced subnetwork of $\mathcal{R}_{2}$, and write $\mathcal{R}_{1}\leq\mathcal{R}_{2}$, if $PN(\mathcal{R}_{1})$ is a vertex-induced subgraph of $PN(\mathcal{R}_{2})$. Clearly, the induced subnetwork relationship induces a partial order on the set of CRNs as discussed in [40, 8]. Note that if $\mathcal{R}_{1}\leq\mathcal{R}_{2}$, the occurrence of a reaction $R$ in both $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ does not mean that $R$ is, physically speaking, the same reaction with the same source and target complexes in $\mathcal{R}_{2}$ as in $\mathcal{R}_{1}$: identifying reactions with (labelled) vertices in a PN graph means that they maintain their identity as graph theoretic modifications are carried out equivalent to inserting or deleting species. If $\mathcal{R}_{1}\leq\mathcal{R}_{2}$, and both have the same set of species, we say that $\mathcal{R}_{1}$ is a reaction-induced subnetwork of $\mathcal{R}_{2}$, and write $\mathcal{R}_{1}\leq_{R}\mathcal{R}_{2}$. If $\mathcal{R}_{1}\leq\mathcal{R}_{2}$, and both have the same set of reactions, we say that $\mathcal{R}_{1}$ is a species-induced subnetwork of $\mathcal{R}_{2}$, and write $\mathcal{R}_{1}\leq_{S}\mathcal{R}_{2}$. Some of the definitions are illustrated in the following example.

Example 3.3**.**

Consider the following CRN:

[TABLE]

$\mathcal{R}$* involves $4$ species $\{W,X,Y,Z\}$, $6$ complexes $\{W,X,2Y,X+Y,Y+Z,W+Z\}$ and four (irreversible) reactions. The complex graph of $\mathcal{R}$ is shown below to the left and the PN graph in the centre. Removing the highlighted vertices and their incident arcs leads to the induced subnetwork*

[TABLE]

represented in unlabelled form with species vertices as open circles and reaction vertices as filled circles to the right. Note that two reactions and a species were removed from $\mathcal{R}$ to obtain $\mathcal{R}_{1}$, and so the subnetwork is neither species-induced nor reaction-induced.

\,X+Y\,$$\,2Y\,$$\,Y+Z\,$$\,W\,$$\,X\,$$\,W+Z\,

Z$$\scriptstyle{4}$$\scriptstyle{3}$$\scriptstyle{1}$$\scriptstyle{2}$$W$$X$$Y$$\scriptstyle{2}

\scriptstyle{2}$${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\scriptstyle{2}}

To preview the nature of results to follow, the motif on the right leads to stable periodic behaviour in fully open CRNs with mass action kinetics, and so the fully open extension of $\mathcal{R}$ with mass action kinetics admits an SPPO as a consequence of the presence of this motif.

3.2 ODE models of CRNs: basic definitions

We take the concentrations of chemical species to be nonnegative real numbers. Consider a CRN $\mathcal{R}$ involving $n$ chemical species $X_{1},\ldots,X_{n}$ with corresponding concentration vector $x=(x_{1},\ldots,x_{n})^{\mathrm{t}}$, and $m$ irreversible reactions between the species. Orderings on the species and reactions are arbitrary but assumed fixed. Define nonnegative $n\times m$ matrices $\Gamma_{l}$ and $\Gamma_{r}$ as follows: $(\Gamma_{l})_{ij}$ (resp., $(\Gamma_{r})_{ij}$) is the stoichiometry of species $X_{i}$ on the left (resp., right) of reaction $j$. The stoichiometric matrix of $\mathcal{R}$ is $\Gamma=\Gamma_{r}-\Gamma_{l}$. The $j$th column of $\Gamma$ is termed the reaction vector for the $j$th reaction.

If the reactions of $\mathcal{R}$ proceed with rates $v_{1}(x),v_{2}(x),\ldots,v_{m}(x)$, we define the rate function of $\mathcal{R}$ to be $v(x)=(v_{1}(x),v_{2}(x),\ldots,v_{m}(x))^{\mathrm{t}}$. The evolution of the species concentrations is then governed by the ODE:

[TABLE]

If $v$ is defined and $C^{1}$ on $\mathbb{R}^{n}_{\gg 0}$ then (3.8) defines a local flow on $\mathbb{R}^{n}_{\gg 0}$, while if $v$ is defined and $C^{1}$ on $\mathbb{R}^{n}_{\geq 0}$ (namely, on an open subset of $\mathbb{R}^{n}$ containing $\mathbb{R}^{n}_{\geq 0}$) then, under physically reasonable assumptions on $v$ which ensure $\mathbb{R}^{n}_{\geq 0}$ is forward invariant, (3.8) defines a local semiflow on $\mathbb{R}^{n}_{\geq 0}$. See the introductory chapter of [41] for definitions of local flows (there termed “local dynamical systems”) and local semiflows (there termed “local semi-dynamical systems”).

$\mathrm{im}\,\Gamma$ is referred to as the stoichiometric subspace of the CRN. The nonempty intersection of a coset of $\mathrm{im}\,\Gamma$ with $\mathbb{R}^{n}_{\geq 0}$ (resp., $\mathbb{R}^{n}_{\gg 0}$) is a stoichiometry class (resp., positive stoichiometry class) of the CRN. If $\mathbb{R}^{n}_{\gg 0}$ (resp., $\mathbb{R}^{n}_{\geq 0}$) is forward invariant under the evolution defined by (3.8), then positive stoichiometry classes (resp., stoichiometry classes) are invariant under (3.8).

3.3 Kinetics

In order to state the results to follow with maximum applicability, we need some discussion of the rate functions of CRNs, namely the allowed functions $v$ in (3.8). The reader familiar with and primarily interested in mass action kinetics can skip directly to Proposition 4.1 below.

Given a CRN $\mathcal{R}$ with evolution governed by (3.8) we may assume that $v(x)$ belongs to some set of functions $\mathcal{K}$ with domain $\mathbb{R}^{n}_{\gg 0}$ and codomain $\mathbb{R}^{m}$. We refer to $\mathcal{K}$ as the kinetics of $\mathcal{R}$ and to the pair $(\mathcal{R},\mathcal{K})$ as a “CRN with kinetics”. $\mathcal{K}$ may be finitely parameterised or a larger class of functions. Given a CRN with kinetics $(\mathcal{R},\mathcal{K})$, and a given reaction $R$ in $\mathcal{R}$, the set of reaction rates for $R$ allowed by $\mathcal{K}$ is denoted $\mathcal{K}^{(R)}$.

When discussing kinetics it is assumed that a CRN consists of irreversible reactions: the allowed rates of a reversible reaction are derived by considering it as a pair of irreversible reactions. In each case below we assume that the CRN $\mathcal{R}$ involves $n$ species and $m$ (irreversible) reactions, and the $n\times m$ matrices $\Gamma_{l}$, $\Gamma_{r}$ and $\Gamma$ are defined as above. The following is a very large class of kinetics.

Definition 3.4** (Positive general kinetics).**

A rate function $v$ for $\mathcal{R}$ belongs to the class of positive general kinetics if and only if $v(x)$ is defined, positive-valued, and $C^{1}$ on $\mathbb{R}^{n}_{\gg 0}$ and satisfies for each $x\in\mathbb{R}^{n}_{\gg 0}$: (i) $\frac{\partial v_{j}}{\partial x_{i}}>0$ if species $X_{i}$ occurs on the left of reaction $j$, (ii) $\frac{\partial v_{j}}{\partial x_{i}}=0$ if species $X_{i}$ does not occur on the left of reaction $j$. Conditions (i) and (ii) can together be rephrased as “the matrix $Dv(x)$ of partial derivatives of $v$ has the same sign pattern as $\Gamma_{l}$”. (Note that in [37] condition (ii) was not spelled out explicitly, although it is implicit throughout.)

Definition 3.5** (General kinetics).**

A rate function $v$ for $\mathcal{R}$ belongs to the class of general kinetics if and only if $v(x)$ is defined and $C^{1}$ on $\mathbb{R}^{n}_{\geq 0}$, satisfies all the restrictions of positive general kinetics on $\mathbb{R}^{n}_{\gg 0}$, and $v_{j}(x)=0$ if and only if $x_{i}=0$ for some species $X_{i}$ occurring on the left of reaction $j$. $\mathbb{R}^{n}_{\geq 0}$ can easily be shown to be positively invariant for (3.8) under the assumption of general kinetics.

Definition 3.6** (Power-law kinetics, physical power-law kinetics, mass action kinetics).**

A rate function $v$ for $\mathcal{R}$ belongs to the class of power-law kinetics if there exist $K\in\mathbb{R}^{m}_{\gg 0}$ and $M\in\mathbb{R}^{m\times n}$ such that $v(x)=K\circ x^{M}$ (recall Notation 1.6 and 1.7 above). $K$ is the vector of rate constants and $M$ is the matrix of exponents. $v$ belongs to the class of physical power-law kinetics if, additionally, $M$ has the same sign pattern as $\Gamma_{l}^{\mathrm{t}}$, and of mass action kinetics if $M=\Gamma_{l}^{\mathrm{t}}$.

Remark 3.7* (Fixed power-law kinetics).*

Stating only that $\mathcal{R}$ has power-law kinetics, or physical power-law kinetics, implies that the entries of both $M$ and of $K$ are parameters which may vary. Stating that $\mathcal{R}$ has fixed power-law kinetics means that $M$ is fixed, while only the entries of $K$ are parameters which may vary.

Remark 3.8* (Relationships between kinetic classes).*

It is easily seen that physical power-law kinetics is a subclass of positive general kinetics, and that mass action kinetics is a particular case of fixed, physical power-law kinetics, and also of general kinetics. Further inclusions amongst classes of kinetics are detailed in [37].

When we refer to $(\mathcal{R},\mathcal{K})$ as a “CRN with mass action kinetics”, or more briefly a “mass action CRN”, this means that the set of allowed rate functions $\mathcal{K}$ is precisely that given by the assumption of mass action kinetics. A similar comment applies to other classes of kinetics.

Definition 3.9** (Derived power-law kinetics).**

Let $(\mathcal{R}_{1},\mathcal{K}_{1})$ and $(\mathcal{R}_{2},\mathcal{K}_{2})$ be CRNs with fixed power-law kinetics and corresponding matrices of exponents $M_{1}$ and $M_{2}$. Let $\mathcal{R}_{2}$ have $n$ species and $m$ reactions and let $\mathcal{R}_{1}$ be an induced subnetwork of $\mathcal{R}_{2}$ with species indexed from $\alpha\subseteq\{1,\ldots,n\}$ and reactions indexed from $\beta\subseteq\{1,\ldots,m\}$. Then $\mathcal{K}_{2}$ is derived from $\mathcal{K}_{1}$ if $M_{1}=M_{2}(\beta|\alpha)$, where $M_{2}(\beta|\alpha)$ is the submatrix of $M_{2}$ with rows from $\beta$ and columns from $\alpha$.

Definition 3.10** (Scaling invariant kinetics).**

Let $(\mathcal{R},\mathcal{K})$ be a CRN with kinetics. Then $\mathcal{K}$ is scaling invariant if, for each reaction $R$ of $\mathcal{R}$ and each $\epsilon>0$, $F\in\mathcal{K}^{(R)}$ implies that $\epsilon F\in\mathcal{K}^{(R)}$.

Remark 3.11* (Scaling invariant kinetics).*

A CRN with any reasonable kinetics, including positive general kinetics, power-law kinetics, physical power-law kinetics, or any fixed power-law kinetics (including mass action) has kinetics which is scaling invariant: if $v_{j}(x)$ is an allowed reaction rate from one of these classes for reaction $j$, then so is $\epsilon v_{j}(x)$ for each $\epsilon>0$.

3.4 Extending the kinetics of an induced subnetwork

Consider CRNs $\mathcal{R}_{1}\leq\mathcal{R}_{2}$ with $\mathcal{R}_{1}$ given kinetics $\mathcal{K}_{1}$. Are there natural ways of “extending” $\mathcal{K}_{1}$ to a kinetics $\mathcal{K}_{2}$ for $\mathcal{R}_{2}$? For example, it is reasonable and often mathematically convenient to assume that:

• •

Where a reaction of $\mathcal{R}_{1}$ occurs with the same source and target complexes in $\mathcal{R}_{2}$, the rates for this reaction allowed by $\mathcal{K}_{2}$ should include those allowed by $\mathcal{K}_{1}$.

• •

Where a reaction of $\mathcal{R}_{1}$ occurs in $\mathcal{R}_{2}$ with some new species involved, fixing the concentrations of the new species at some positive values should give back (at least) all the rate functions allowed by $\mathcal{K}_{1}$.

These notions are formalised in the following two definitions.

Definition 3.12** (Reaction-extensions).**

Consider CRNs with kinetics $(\mathcal{R}_{1},\mathcal{K}_{1})$ and $(\mathcal{R}_{2},\mathcal{K}_{2})$. Then $(\mathcal{R}_{2},\mathcal{K}_{2})$ is a reaction-extension of $(\mathcal{R}_{1},\mathcal{K}_{1})$, written $(\mathcal{R}_{1},\mathcal{K}_{1})\leq_{R}(\mathcal{R}_{2},\mathcal{K}_{2})$, if $\mathcal{R}_{1}$ is a reaction-induced subnetwork of $\mathcal{R}_{2}$ and $\mathcal{K}_{1}^{(R)}\subseteq\mathcal{K}_{2}^{(R)}$ for each reaction $R$ occurring in both $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$. In other words, reactions which $\mathcal{R}_{2}$ inherits from $\mathcal{R}_{1}$ are allowed (at least) all the rate functions allowed by $\mathcal{K}_{1}$.

Definition 3.13** (Species-extensions).**

Consider CRNs with kinetics $(\mathcal{R}_{1},\mathcal{K}_{1})$ and $(\mathcal{R}_{2},\mathcal{K}_{2})$ and suppose that $\mathcal{R}_{1}$ is a species-induced subnetwork of $\mathcal{R}_{2}$. Let $\mathcal{R}_{2}$ have $n_{2}$ species $X_{1},\ldots,X_{n_{2}}$ and assume, without loss of generality, that the species of $\mathcal{R}_{1}$ are $X_{1},\ldots,X_{n_{1}}$ where $n_{1}\leq n_{2}$. Let $\hat{x}=(x_{1},\ldots,x_{n_{1}})$ and $\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{x}}}=(x_{n_{1}+1},\ldots,x_{n_{2}})$. Then $(\mathcal{R}_{2},\mathcal{K}_{2})$ is a species-extension of $(\mathcal{R}_{1},\mathcal{K}_{1})$, written $(\mathcal{R}_{1},\mathcal{K}_{1})\leq_{S}(\mathcal{R}_{2},\mathcal{K}_{2})$, if for each $v(\hat{x})\in\mathcal{K}_{1}$ there exists $w(\hat{x},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{x}}})\in\mathcal{K}_{2}$ such that $w(\hat{x},\mathbf{1})=v(\hat{x})$, and such that if $v$ is $C^{k}$ on $\mathbb{R}^{n_{1}}_{\gg 0}$ (resp., $\mathbb{R}^{n_{1}}_{\geq 0}$), then $w$ is $C^{k}$ on $\mathbb{R}^{n_{2}}_{\gg 0}$ (resp., $\mathbb{R}^{n_{2}}_{\geq 0}$). (In the case $n_{1}=n_{2}$ and $\mathcal{K}_{1}\subseteq\mathcal{K}_{2}$ we take $(\mathcal{R}_{1},\mathcal{K}_{1})\leq_{S}(\mathcal{R}_{2},\mathcal{K}_{2})$ to be trivially true.)

Lemma 3.14** (Species-extensions).**

Let $\mathcal{R}_{1}\leq_{S}\mathcal{R}_{2}$. Then the CRNs with kinetics $(\mathcal{R}_{1},\mathcal{K}_{1})$ and $(\mathcal{R}_{2},\mathcal{K}_{2})$ satisfy $(\mathcal{R}_{1},\mathcal{K}_{1})\leq_{S}(\mathcal{R}_{2},\mathcal{K}_{2})$ if any of the following hold:

$\mathcal{K}_{1}$* and $\mathcal{K}_{2}$ are both given by positive general kinetics.* 2. 2.

$\mathcal{K}_{1}$* and $\mathcal{K}_{2}$ are both given by power-law kinetics.* 3. 3.

$\mathcal{K}_{1}$* and $\mathcal{K}_{2}$ are both given by physical power-law kinetics.* 4. 4.

$\mathcal{K}_{1}$* and $\mathcal{K}_{2}$ are both given by mass action kinetics.* 5. 5.

$\mathcal{K}_{1}$* and $\mathcal{K}_{2}$ are both given by fixed power-law kinetics with $\mathcal{K}_{2}$ derived from $\mathcal{K}_{1}$ (see Definition 3.9).*

Proof 2*.*

Using the notation in Definition 3.13, for each rate function $v\in\mathcal{K}_{1}$, set $w(\hat{x},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{x}}})=v(\hat{x})\circ\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{x}}}^{M}$ where, in cases (1) to (4), $M$ consists of the final $n_{2}-n_{1}$ columns of $\Gamma_{l}^{\mathrm{t}}$, while in case (5) $M$ consists of the final $n_{2}-n_{1}$ columns of $M_{2}$, the matrix of exponents of $\mathcal{R}_{2}$. In each case it is easily seen that $w\in\mathcal{K}_{2}$, that $w(\hat{x},\mathbf{1})=v(\hat{x})$, and that if $v$ is $C^{k}$ on $\mathbb{R}^{n_{1}}_{\gg 0}$ (resp., $\mathbb{R}^{n_{1}}_{\geq 0}$ in case 4), then $w$ is $C^{k}$ on $\mathbb{R}^{n_{2}}_{\gg 0}$ (resp., $\mathbb{R}^{n_{2}}_{\geq 0}$ in case 4).

Definition 3.15** (Species-reaction-extensions).**

Let $\mathcal{R}_{1}\leq\mathcal{R}_{2}$ and consider CRNs with kinetics $(\mathcal{R}_{1},\mathcal{K}_{1})$ and $(\mathcal{R}_{2},\mathcal{K}_{2})$. Observe that there is a uniquely defined CRN $\mathcal{R}^{\prime}$ satisfying $\mathcal{R}_{1}\leq_{S}\mathcal{R}^{\prime}\leq_{R}\mathcal{R}_{2}$. ($\mathcal{R}^{\prime}$ is obtained by inserting all missing species into the reactions of $\mathcal{R}_{1}$, but without adding any new reactions.) Then $(\mathcal{R}_{2},\mathcal{K}_{2})$ is a species-reaction-extension of $(\mathcal{R}_{1},\mathcal{K}_{1})$ if there exists $\mathcal{K}^{\prime}$ such that

[TABLE]

Intuitively, we first add in missing species and extend the kinetics of any modified reactions consistent with the species-extension condition, and then add in any remaining missing reactions.

4 Results on the inheritance of NPPOs and SPPOs

A CRN with kinetics $(\mathcal{R},\mathcal{K})$ admits an NPPO (resp., SPPO) if there exists some rate function $v\in\mathcal{K}$ s.t. the associated ODE system (3.8) has an NPPO (resp., SPPO). A broad question is when, given CRNs with kinetics $(\mathcal{R}_{1},\mathcal{K}_{1})$ and $(\mathcal{R}_{2},\mathcal{K}_{2})$ related in some natural way, knowledge that one admits an NPPO (resp., SPPO) allows us to predict the same for the other. Four “inheritance” theorems in this direction will be proved below under varying kinetic assumptions. For the reader primarily interested in mass action kinetics, these can be summarised in a single corollary:

Proposition 4.1**.**

Let $\mathcal{R}$ and $\mathcal{R}^{\prime}$ be CRNs, and suppose that $\mathcal{R}$ admits an NPPO (resp., SPPO) with mass action kinetics. Suppose that we create $\mathcal{R}^{\prime}$ from $\mathcal{R}$ by

adding to $\mathcal{R}$ a new reaction with reaction vector in the span of reaction vectors of $\mathcal{R}$; or 2. 2.

taking the fully open extension of $\mathcal{R}$; or 3. 3.

adding into some reactions of $\mathcal{R}$ a new species $Y$ which occurs with the same stoichiometry on both sides of each reaction in which it participates; or 4. 4.

adding into reactions of $\mathcal{R}$ a new species $Y$ in any way, while also adding the new reaction $0\rightleftharpoons Y$.

Then, with mass action kinetics, $\mathcal{R}^{\prime}$ admits an NPPO (resp., SPPO).

Proof 3*.*

Claims 1 to 4 are immediate corollaries of Theorems 1 to 4 below and the surrounding remarks. In order to apply the results we need only note that mass action kinetics is polynomial and hence certainly $C^{2}$, is scaling invariant, and that the assumptions imply that $\mathcal{R}^{\prime}$ with mass action kinetics is a species-reaction extension of $\mathcal{R}$ with mass action kinetics. $\square$

Theorems 1 and 2 require only basic regular perturbation theory to prove: in Theorem 1 the application is almost trivial while in Theorem 2 it takes a little more work to set up the problem. Theorem 3 requires essentially no machinery to prove: the proof is almost immediate from the definitions. Theorem 4 requires some results from singular perturbation theory. Theorems 1 and 4 together imply an important corollary about fully open networks spelled out as Proposition 4.8.

In each of the following theorems, $\mathcal{R}$ is a CRN with $m$ reactions involving $n$ species $X_{1},\ldots,X_{n}$ with concentrations $x_{1},\ldots,x_{n}$. $\Gamma$, the stoichiometric matrix of $\mathcal{R}$, has rank $r$, $\Gamma_{0}$ is a matrix whose columns are a basis for $S:=\mathrm{im}\,\Gamma$, and $Q$ is defined by $\Gamma=\Gamma_{0}Q$. Given a periodic orbit $\mathcal{O}$, $S_{\mathcal{O}}:=(\mathcal{O}+S)\cap\mathbb{R}^{n}_{\gg 0}$ is the positive stoichiometry class of $\mathcal{O}$, and $x_{0}$ is some point on $S_{\mathcal{O}}$ (recall Figure 1).

Theorem 1** (Adding a dependent reaction).**

Let $(\mathcal{R},\mathcal{K})$ be a CRN with $C^{1}$ kinetics admitting an NPPO (resp., SPPO). Let $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ be a reaction-extension of $(\mathcal{R},\mathcal{K})$ created by adding to $\mathcal{R}$ a new irreversible reaction with $C^{1}$, scaling invariant, kinetics, and with reaction vector in the span of reaction vectors of $\mathcal{R}$. Then $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ admits an NPPO (resp., SPPO).

Proof 4*.*

Fix the rate function $v\in\mathcal{K}$ such that $\mathcal{R}$ has an NPPO (resp., SPPO) $\mathcal{O}$. Let the new reaction of $\mathcal{R}^{\prime}$ be $a\cdot X\rightarrow a^{\prime}\cdot X$. Define $\alpha=a^{\prime}-a$ and define $c$ by $\alpha=\Gamma_{0}c$. Consistent with the kinetic assumptions, set the rate of the new reaction to be $\epsilon f(x)$ where $f\colon\mathbb{R}^{n}_{\gg 0}\to\mathbb{R}$ is $C^{1}$ and $\epsilon$ is a parameter to be controlled (for example, with mass action kinetics the rate would be $\epsilon x^{a}$). The evolution of $\mathcal{R}^{\prime}$ is then governed by:

[TABLE]

Define $h\colon\mathbb{R}^{r}\to x_{0}+S$ by $h(z):=x_{0}+\Gamma_{0}z$. $h$ is an affine bijection between $h^{-1}(S_{\mathcal{O}})$ and $S_{\mathcal{O}}$ and defines local coordinates on $S_{\mathcal{O}}$ via $x=h(z)$. $z$ evolves according to

[TABLE]

By definition, (4.10) has the hyperbolic (resp., linearly stable) periodic orbit $\mathcal{O}^{\prime}:=h^{-1}(\mathcal{O})$. By Lemma 2.1 there exists $\epsilon_{0}>0$ s.t. for $\epsilon\in(-\epsilon_{0},\epsilon_{0})$ (4.10) has a hyperbolic (resp., linearly stable) periodic orbit $\mathcal{O}^{\prime}_{\epsilon}$. Thus, for $\epsilon\in(-\epsilon_{0},\epsilon_{0})$, (4.9) has the NPPO (resp., SPPO) $\mathcal{O}_{\epsilon}:=h(\mathcal{O}^{\prime}_{\epsilon})$. $\square$

Remark 4.2* (Adding the reverse of a reaction).*

Clearly, by Theorem 1, given a CRN $\mathcal{R}$ with kinetics from any $C^{1}$ class admitting an NPPO (resp., SPPO), adding the reverse of any existing reaction to $\mathcal{R}$ with $C^{1}$, scaling invariant, kinetics preserves this property. Thus if a CRN with, say, mass action kinetics admits an NPPO (resp., SPPO), then so does the corresponding reversible CRN with mass action kinetics.

Remark 4.3* (Preservation of bifurcations when dependent reactions are added).*

In [32] Conradi and Shiu posed the question of whether Hopf bifurcations in CRNs are preserved when some irreversible reactions are made reversible. In fact, any generic bifurcation ([35] or [42] for example) survives the addition of dependent reactions with sufficiently smooth, scaling-invariant, kinetics. Although Theorem 1 is not about bifurcation per se, the key idea in its proof is the construction of local coordinates on a stoichiometry class $S$ so that the vector field of $\mathcal{R}^{\prime}$ in these local coordinates is a perturbation of the original vector field of $\mathcal{R}$. Suppose that some $C^{r}$ $k$-parameter family of vector fields $\mathcal{F}_{\lambda}$ on $S$ associated with $\mathcal{R}$ admits a nondegenerate codimension-$k$ bifurcation at $(x_{0},\lambda_{0})$. Then, as we see from (4.10), addition of a new dependent reaction with $C^{r}$, scaling-invariant, kinetics gives rise, for each fixed $\epsilon$, to a new $C^{r}$, $k$-parameter, family $\mathcal{F}^{\epsilon}_{\lambda}$ of vector fields for $\mathcal{R}^{\prime}$, $C^{r}$ close to $\mathcal{F}_{\lambda}$; for $r$ sufficiently large, and $\epsilon$ sufficiently small, the family $\mathcal{F}^{\epsilon}_{\lambda}$ will admit the same nondegenerate bifurcation. Analogous remarks apply to the other network modifications detailed in the theorems to follow. As a practical note, confirming that a given CRN does indeed admit a generic Hopf bifurcation at some parameter values is not always entirely straightfoward, as it may involve approximation of a parameter-dependent center manifold in order to confirm the nondegeneracy conditions.

Theorem 2** (Adding inflows and outflows of all species).**

Let $(\mathcal{R},\mathcal{K})$ be a CRN with $C^{1}$ kinetics admitting an NPPO (resp., SPPO). Suppose that $\mathcal{R}$ includes no flow reactions (i.e., no reactions of the form $0\rightarrow X_{i}$ or $X_{i}\rightarrow 0$). Let $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ be a reaction-extension of $(\mathcal{R},\mathcal{K})$ created by adding to $\mathcal{R}$ all the reactions $0\rightleftharpoons X_{i}$ ($i=1,\ldots,n$) with kinetics from a class including mass action kinetics. Then $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ admits an NPPO (resp., SPPO).

Proof 5*.*

Fix the rate function $v\in\mathcal{K}$ such that $\mathcal{R}$ has an NPPO (resp., SPPO) $\mathcal{O}$. Treat the $i$th inflow-outflow reaction as a single reversible reaction with mass action kinetics and forward and backwards rate constants $\epsilon(x_{0})_{i}$ and $\epsilon$ respectively. The evolution of $\mathcal{R}^{\prime}$ is then governed by:

[TABLE]

Let $\Gamma_{0}^{\prime}=[\Gamma_{0}|\Gamma_{1}]$ where $\Gamma_{1}$ is any $n\times(n-r)$ matrix chosen so that $\Gamma_{0}^{\prime}$ has rank $n$. Observe that \Gamma=\Gamma_{0}^{\prime}\left(\begin{array}[]{c}Q\\ 0\end{array}\right). Define new coordinates $z=(\hat{z},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{z}}})\in\mathbb{R}^{r}\times\mathbb{R}^{n-r}$ by $x=h(z):=x_{0}+\Gamma_{0}^{\prime}z$. $h$ is an affine bijection between $W:=h^{-1}(\mathbb{R}^{n}_{\gg 0})\subseteq\mathbb{R}^{r}\times\mathbb{R}^{n-r}$ and $\mathbb{R}^{n}_{\gg 0}$, and $z$ evolves according to

[TABLE]

Define $W_{1}:=W\cap(\mathbb{R}^{r}\times\{0\})$, so that $h(W_{1})=S_{\mathcal{O}}$. Define $\mathcal{O}^{\prime}:=h^{-1}(\mathcal{O})\subseteq W_{1}$ and define $\overline{\mathcal{O}}\subseteq\mathbb{R}^{r}$ by $\overline{\mathcal{O}}\times\{0\}=\mathcal{O}^{\prime}$. $W_{1}$ is locally invariant for (4.11), and restricting (4.11) to $W_{1}$ gives the differential equation:

[TABLE]

By definition, $\mathcal{O}$ is an NPPO (resp., SPPO) of $\mathcal{R}$ if and only if $\overline{\mathcal{O}}$ is a hyperbolic (resp., linearly stable) periodic orbit of (4.12). In this case, by Lemma 2.1, there exists $\epsilon_{0}>0$ s.t. for $\epsilon\in(-\epsilon_{0},\epsilon_{0})$, (4.12) has a hyperbolic (resp., linearly stable) periodic orbit $\overline{\mathcal{O}}_{\epsilon}$ close to $\overline{\mathcal{O}}$ with period $T_{\epsilon}$ close to $T$. It remains to show that $\mathcal{O}^{\prime}_{\epsilon}:=\overline{\mathcal{O}}_{\epsilon}\times\{0\}$ is hyperbolic (resp., linearly stable) for (4.11) for all sufficiently small $\epsilon>0$. This will imply immediately that $\mathcal{O}_{\epsilon}:=h(\mathcal{O}^{\prime}_{\epsilon})$ is an NPPO (resp., SPPO) of $\mathcal{R}^{\prime}$.

For each fixed $\epsilon\in(0,\epsilon_{0})$, choose $\psi_{\epsilon}$ to be some solution of (4.12) with initial condition on $\overline{\mathcal{O}}_{\epsilon}$. The variational equation of (4.12) about $\psi_{\epsilon}$ is:

[TABLE]

The fundamental matrix solution $\hat{Z}_{\epsilon}(t)$ of (4.13) with $\hat{Z}_{\epsilon}(0)=I_{r}$ can be written $\hat{Z}_{\epsilon}(t)=A_{\epsilon}(t)e^{tB_{\epsilon}}$ where $A_{\epsilon}(t)$ is a nonsingular periodic matrix of period $T_{\epsilon}>0$ and $B_{\epsilon}$ is a constant matrix. Hyperbolicity (resp., linear stability) of $\overline{\mathcal{O}}_{\epsilon}$ for (4.12) means that $\hat{Z}_{\epsilon}(T_{\epsilon})=e^{T_{\epsilon}B_{\epsilon}}$ has exactly one eigenvalue equal to $1$ with the remaining $r-1$ eigenvalues disjoint from (resp., inside) the unit circle.

For each $\psi_{\epsilon}$ chosen as above, $(\psi_{\epsilon},0)$ is clearly a periodic solution of (4.11), with image $\mathcal{O}^{\prime}_{\epsilon}$. The full variational equation of (4.11) about $(\psi_{\epsilon},0)$ is:

[TABLE]

Our goal is to compute $Z_{\epsilon}(t)$, the fundamental matrix solution of (4.14) satisfying $Z_{\epsilon}(0)=I$. Solving the second equation of (4.14) gives $\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\zeta}}}(t)=e^{-\epsilon t}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\zeta}}}(0)$. Substituting into the first equation of (4.14) gives

[TABLE]

Setting $\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\zeta}}}(0)=0$ gives back (4.13). The above calculations give:

[TABLE]

Here $A(t)$ is some matrix which can be determined by integration but which does not affect the subsequent argument. The characteristic multipliers of $\mathcal{O}^{\prime}_{\epsilon}$ are precisely the eigenvalues of $\hat{Z}_{\epsilon}(T_{\epsilon})$, namely the characteristic multipliers of $\overline{\mathcal{O}}_{\epsilon}$ for (4.12), and the single value $e^{-\epsilon T_{\epsilon}}$ occurring with multiplicity $n-r$. As $\overline{\mathcal{O}}_{\epsilon}$ is a hyperbolic (resp., linearly stable) periodic orbit of (4.12), and $e^{-\epsilon T_{\epsilon}}$ lies inside the unit circle for any $\epsilon>0$, $T_{\epsilon}>0$, $\mathcal{O}^{\prime}_{\epsilon}$ is a hyperbolic (resp., linearly stable) periodic orbit of (4.11), and consequently $\mathcal{O}_{\epsilon}=h(\mathcal{O}^{\prime}_{\epsilon})$ is an NPPO (resp., SPPO) of $\mathcal{R}^{\prime}$. $\square$

Remark 4.4* (Geometric interpretation of Theorem 2, the role of mass action kinetics).*

Inflows and outflows were chosen to guarantee that $S_{\mathcal{O}}$ remained invariant for $\mathcal{R}^{\prime}$: this necessitated mass action kinetics for the flow reactions. The construction ensured that for sufficiently small $\epsilon>0$ $S_{\mathcal{O}}$ is exponentially attracting and the vector field of $\mathcal{R}^{\prime}$ restricted to $S_{\mathcal{O}}$ is $\epsilon$-close to that of $\mathcal{R}$ restricted to $S_{\mathcal{O}}$, ensuring the existence on $S_{\mathcal{O}}$ of a hyperbolic (resp., linearly stable) periodic orbit $\mathcal{O}_{\epsilon}$ close to $\mathcal{O}$.

Remark 4.5* (Theorem 2 and fully open extensions).*

Suppose that $(\mathcal{R},\mathcal{K})$ is any CRN with $C^{1}$ kinetics such that if $v$ is an allowed rate for some reaction $X_{i}\rightarrow 0$ of $\mathcal{R}$, then so is $v+\epsilon x_{i}$ for all sufficiently small $\epsilon>0$, and if $v$ is an allowed rate for some reaction $0\rightarrow X_{i}$ of $\mathcal{R}$, then so is $v+\epsilon$ for all sufficiently small $\epsilon>0$. Then the condition that $\mathcal{R}$ excludes flow reactions can clearly be dropped in Theorem 2. In particular, if $(\mathcal{R},\mathcal{K})$ is any mass action CRN admitting an NPPO (resp., SPPO) then Theorem 2 tells us that its fully open extension with mass action kinetics admits an NPPO (resp., SPPO). The same holds for CRNs with positive general kinetics. However, we cannot arrive at this conclusion for CRNs with arbitrary fixed physical power-law kinetics.

Theorem 3** (Adding a trivial species).**

Let $(\mathcal{R},\mathcal{K})$ be a CRN with $C^{1}$ kinetics admitting an NPPO (resp., SPPO). Let $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ be a species-extension of $(\mathcal{R},\mathcal{K})$ created by adding into some reactions of $\mathcal{R}$ a new species $Y$ with concentration $y$, which occurs with the same stoichiometry on both sides of each reaction in which it participates. Then $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ admits an NPPO (resp., SPPO).

Proof 6*.*

Fix the rate function $v\in\mathcal{K}$ such that $\mathcal{R}$ has an NPPO (resp., SPPO) $\mathcal{O}$. Fix $w\in\mathcal{K}^{\prime}$ such that $w(x,1)=v(x)$, possible by assumption. With this rate function, the evolution of $\mathcal{R}^{\prime}$ is governed by

[TABLE]

(4.15) leaves $\mathbb{R}^{n}_{\gg 0}\times\{y\}$ locally invariant for each $y>0$. Since $w(x,1)=v(x)$, $\mathcal{O}^{\prime}:=\mathcal{O}\times\{1\}$ is a periodic orbit of $\mathcal{R}^{\prime}$. Let $S^{\prime}=S\times\{0\}$ so that $S^{\prime}_{\mathcal{O^{\prime}}}:=(\mathcal{O}^{\prime}+S^{\prime})\cap(\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}_{>0})=S_{\mathcal{O}}\times\{1\}$ is the positive stoichiometry class of $\mathcal{O}^{\prime}$ for $\mathcal{R}^{\prime}$. Define $h\colon\mathbb{R}^{r}\to x_{0}+S$ by $h(z)=x_{0}+\Gamma_{0}z$. $h$ is an affine bijection between $W:=h^{-1}(S_{\mathcal{O}})\subseteq\mathbb{R}^{r}$ and $S_{\mathcal{O}}$, and defines local coordinates on $S_{\mathcal{O}}$ which evolve according to

[TABLE]

By definition, as $\mathcal{O}$ is an NPPO (resp., SPPO), $h^{-1}(\mathcal{O})$ is nondegenerate (resp., linearly stable). Now define $h^{\prime}\colon\mathbb{R}^{r}\to(x_{0}+S)\times\{1\}$ by $h^{\prime}(z)=(x_{0}+\Gamma_{0}z,1)$ and note that $h^{\prime}$ is an affine bijection between $W$ and $S^{\prime}_{\mathcal{O^{\prime}}}$. Moreover $h^{\prime}$ gives rise to precisely the same evolution in local coordinates (since $w(x_{0}+\Gamma_{0}z,1)=v(x_{0}+\Gamma_{0}z)$) and $h^{\prime-1}(\mathcal{O}^{\prime})=h^{-1}(\mathcal{O})$. Thus, by definition, $\mathcal{O}^{\prime}$ is an NPPO (resp., SPPO) of $\mathcal{R}^{\prime}$. $\square$

Theorem 4** (Adding a new species with inflow and outflow).**

Let $(\mathcal{R},\mathcal{K})$ be a CRN with $C^{2}$ kinetics admitting an NPPO (resp., SPPO). Let $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ be a species-reaction-extension of $(\mathcal{R},\mathcal{K})$ created by

1. (i)

adding into the reactions of $\mathcal{R}$ a new species $Y$ with arbitrary stoichiometries; and 2. (ii)

adding the new reaction $0\rightleftharpoons Y$ with $C^{2}$ kinetics belonging to a scaling invariant subset of positive general kinetics.

Then $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ admits an NPPO (resp., SPPO).

Proof 7*.*

Fix the rate function $v\in\mathcal{K}$ such that $\mathcal{R}$ has an NPPO (resp., SPPO) $\mathcal{O}$. As in the proof of Theorem 1 define $h\colon\mathbb{R}^{r}\to x_{0}+S$ by $h(z)=x_{0}+\Gamma_{0}z$ and note that $h$ is an affine bijection between the open set $h^{-1}(S_{\mathcal{O}})\subseteq\mathbb{R}^{r}$ and $S_{\mathcal{O}}$. $h$ defines local coordinates on $S_{\mathcal{O}}$ via $x=h(z)$, and $z$ evolves according to

[TABLE]

(4.16) has a hyperbolic (resp., linearly stable) periodic orbit $\mathcal{O}^{\prime}=h^{-1}(\mathcal{O})$. The assumptions on the kinetics mean that:

The new rate function $w(x,y)$ of the existing reactions can be chosen to satisfy $w(x,1)=v(x)$. 2. 2.

There exists a $C^{2}$ function $f\colon\mathbb{R}_{>0}\to\mathbb{R}_{>0}$ satisfying $f(1)=1$ and $f^{\prime}(y)>0$ for all $y>0$ and such that we may choose the rate of $0\rightleftharpoons Y$ to be $\frac{1}{\epsilon}(1-f(y))$ where $\epsilon>0$ is a parameter to be controlled.

With these choices, $\mathcal{R}^{\prime}$ gives rise to the following singularly perturbed system:

[TABLE]

Here $s_{i}$ is the net change in the stoichometry of $Y$ in the $i$th reaction of $\mathcal{R}^{\prime}$, and $s:=(s_{1},\ldots,s_{m})^{\mathrm{t}}$. For any fixed $\epsilon>0$, rescaling time in the “slow time system” (4.17) gives the “fast time system”:

[TABLE]

Define $h^{\prime}\colon\mathbb{R}^{r}\times\mathbb{R}\to(x_{0}+S)\times\mathbb{R}$ by $h^{\prime}(z,y):=(h(z),y)=(x_{0}+\Gamma_{0}z,y)$. Note that $h^{\prime}$ is an affine bijection between $h^{-1}(S_{\mathcal{O}})\times\mathbb{R}_{>0}$ and $S_{\mathcal{O}}\times\mathbb{R}_{>0}$ and defines local coordinates on $S_{\mathcal{O}}\times\mathbb{R}_{>0}$ via $(x,y)=(h(z),y)$. In $(z,y)$ coordinates the slow time system (4.17) becomes:

[TABLE]

while the fast time system (4.18) becomes:

[TABLE]

(4.19) is the decoupled differential-algebraic system

[TABLE]

(as $f(y)=1$ if and only if $y=1$, and $w(x_{0}+\Gamma_{0}z,1)=v(x_{0}+\Gamma_{0}z)$). We observe that

The vector field $H(z)$ has a hyperbolic (resp., linearly stable) periodic orbit $\mathcal{O}^{\prime}$ by assumption, and hence (4.19) has a periodic orbit $\overline{\mathcal{O}}:=\mathcal{O}^{\prime}\times\{1\}$. 2. 2.

$y=1$ is a linearly stable equilibrium of $\dot{y}=1-f(y)$ or, equivalently, the Jacobian matrix of (4.20) evaluated at $y=1,\epsilon=0$, namely,

[TABLE]

has a single nontrivial eigenvalue $-f^{\prime}(1)<0$.

By Theorems 13.1 and 13.2 in [36], observations (1) and (2) together tell us that there exists $\epsilon_{0}>0$ s.t. that for $\epsilon\in(0,\epsilon_{0})$, (4.19) has a hyperbolic (resp., linearly stable) periodic orbit $\overline{\mathcal{O}}_{\epsilon}$ close to $\overline{\mathcal{O}}$. Thus, for $\epsilon\in(0,\epsilon_{0})$, $\mathcal{R}^{\prime}$ has an NPPO (resp., SPPO) $\mathcal{O}_{\epsilon}:=h^{\prime}(\overline{\mathcal{O}}_{\epsilon})$. $\square$

Remark 4.6* (Geometrical interpretation of the proof of Theorem 4).*

The differential algebraic system (4.19) defines a local flow on the $r$-dimensional (smooth) manifold $\mathcal{Y}:=h^{-1}(\mathcal{S}_{\mathcal{O}})\times\{1\}$ which includes the periodic orbit $\overline{\mathcal{O}}$. Let $\mathcal{Y}_{0}$ be some compact subset of $\mathcal{Y}$ containing $\overline{\mathcal{O}}$. The theory developed by Fenichel [36] shows (roughly, and omitting a myriad of technical details) that for sufficiently small $\epsilon>0$ (4.19) has an $r$-dimensional locally invariant manifold $\mathcal{Y}_{\epsilon}$ close to $\mathcal{Y}_{0}$. The vector field of (4.19) restricted to $\mathcal{Y}_{\epsilon}$ is $\epsilon$-close to that of (4.19) on $\mathcal{Y}_{0}$ and consequently, by regular perturbation theory, for sufficiently small $\epsilon\neq 0$, (4.19) has a periodic orbit $\overline{\mathcal{O}}_{\epsilon}$ on $\mathcal{Y}_{\epsilon}$ close to $\overline{\mathcal{O}}$. The technical assumption that all vector fields involved are $C^{2}$ is to ensure that the family of vector fields on $\mathcal{Y}_{\epsilon}$ is $C^{1}$, allowing use of regular perturbation theory. The $r-1$ nontrivial Floquet multipliers of $\overline{\mathcal{O}}_{\epsilon}$ relative to $\mathcal{Y}_{\epsilon}$ are close to those of $\overline{\mathcal{O}}$ relative to $\mathcal{Y}$ which, by assumption, are disjoint from (resp., inside) the unit circle. Meanwhile, the single Floquet multiplier of $\overline{\mathcal{O}}_{\epsilon}$ transverse to $\mathcal{Y}_{\epsilon}$ lies inside the unit circle as a consequence of the fact that $-f^{\prime}(1)<0$.

Remark 4.7* (Kinetic assumptions in Theorem 4).*

The added flow reaction $0\rightleftharpoons Y$ in Theorem 4 may have, for example, mass action kinetics, positive general kinetics, physical power-law kinetics, or any fixed physical power-law kinetics (these all define scaling invariant subsets of positive general kinetics).

Theorem 4, combined with Theorem 1 allows us to deduce an important corollary:

Proposition 4.8** (Inheritance in fully open species-reaction extensions).**

Let $(\mathcal{R},\mathcal{K})$ be a fully open CRN with $C^{2}$ kinetics admitting an NPPO (resp., SPPO). Let $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ be a fully open CRN with kinetics, which is a species-reaction extension of $(\mathcal{R},\mathcal{K})$ (Definition 3.15), and such that for each new reaction $R$ in $\mathcal{R}^{\prime}$, $\mathcal{K}^{\prime(R)}$ is $C^{2}$, and belongs to a scaling invariant subset of positive general kinetics. Then $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ admits an NPPO (resp., SPPO).

Proof 8*.*

Let $\mathcal{R}$ have $n_{1}$ species and $m_{1}$ non-flow reactions (i.e., reactions not of the form $0\rightarrow X_{i}$ or $X_{i}\rightarrow 0$), and $\mathcal{R}^{\prime}$ have $n_{2}$ species and $m_{2}$ non-flow reactions. We can construct $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ from $(\mathcal{R},\mathcal{K})$ via a sequence of steps as follows:

1. (i)

Beginning with $\mathcal{R}$, for each absent species $X_{j}$ (if any) we add the species to all existing reactions and add $0\rightleftharpoons X_{j}$. The kinetic assumptions ensure that this corresponds to $n_{2}-n_{1}$ applications of Theorem 4, one for each absent species. Note that, as $\mathcal{R}$ is fully open, the new CRN created at each stage is fully open. 2. (ii)

We add each remaining absent reaction (if any). The kinetic assumptions ensure that this corresponds to $m_{2}-m_{1}$ applications of Theorem 1, one for each reaction added. Theorem 1 applies because a fully open CRN has stoichiometric subspace which is the whole state space, and hence any added reaction is a dependent reaction.

We can see the above procedure as constructing a sequence of intermediate (fully-open) CRNs with kinetics, beginning with $(\mathcal{R},\mathcal{K})$ and terminating with $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$:

[TABLE]

(Here $p=n_{2}+m_{2}-n_{1}-m_{1}$.) If $(\mathcal{R},\mathcal{K})$ admits an NPPO (resp., SPPO), then each step of the above procedure preserves this property, and consequently $(\mathcal{R}^{\prime},\mathcal{K}^{\prime})$ admits an NPPO (resp., SPPO). $\square$

Remark 4.9* (Kinetic assumptions in Proposition 4.8, proof of Proposition 1.1).*

The somewhat unwieldy kinetic assumptions in Proposition 4.8 are in order to maximise generality. They are satisfied if $\mathcal{R}\leq\mathcal{R}^{\prime}$ and, for example,

Both $\mathcal{R}$ and $\mathcal{R}^{\prime}$ have mass action kinetics. 2. 2.

Both $\mathcal{R}$ and $\mathcal{R}^{\prime}$ have physical power-law kinetics. 3. 3.

$\mathcal{R}$ has any fixed power-law kinetics and $\mathcal{R}^{\prime}$ has any power-law kinetics derived from that of $\mathcal{R}$ (see Definition 3.9). 4. 4.

Both $\mathcal{R}$ and $\mathcal{R}^{\prime}$ have $C^{2}$ positive general kinetics.

Thus, in particular Proposition 1.1 follows immediately from Proposition 4.8. Explorations in Section 5 are carried out using Proposition 4.8 with $\mathcal{R}$ and $\mathcal{R}^{\prime}$ both given mass action kinetics or both given physical power-law kinetics.

In the light of Proposition 4.8, and adapting the terminology of [4], the following definitions make sense.

Definition 4.10** (Atoms of oscillation, atoms of stable oscillation).**

A fully open mass action CRN which admits an NPPO (resp., SPPO), and which is minimal with respect to the induced subnetwork ordering amongst fully open mass action CRNs admitting NPPOs (resp., SPPOs), is referred to as a fully open mass action atom of oscillation (resp., stable oscillation). Atoms with respect to other classes of kinetics, such as physical power-law kinetics, are similarly defined.

Observe that Definition 4.10 is restricted to fully open CRNs as the presence of an oscillatory induced subnetwork in a general CRN does not necessarily imply oscillation; an example is provided in the concluding section (Example 6.1). Note also that, as in the case of multistationarity [8], a fully open mass action atom of oscillation with respect to the induced subnetwork ordering may not be minimal with respect to other, better partial orders. Note finally that a fully open mass action atom of oscillation may include an induced subnetwork admitting an NPPO but which is not fully open; thus if we do not restrict attention to fully open CRNs, fully open atoms of oscillation may not be minimal oscillatory CRNs even with respect to the induced subnetwork ordering.

5 The occurrence of stable oscillation in small, fully open, CRNs

A fully open CRN is taken to be “small” if it has few species, few non-flow reactions, and is at most bimolecular, namely the total stoichiometry of all species on each side of every reaction is no more than two. The goal of this section is to provide some lower bounds on the frequency with which small fully open CRNs admit SPPOs under the assumptions of (i) mass action kinetics and (ii) physical power-law kinetics. This is done via a mixture of basic analysis, numerical simulation, and application of the inheritance result in Proposition 4.8.

Define a $(k,l)$ CRN to be a fully open, at most bimolecular, CRN with $k\geq 1$ species and $l\geq 0$ irreversible non-flow reactions. It is easy to see that $(1,l)$ and $(k,0)$ CRNs can admit no nontrivial periodic orbits for any reasonable kinetics: if $k=1$ then regardless of the kinetics (3.8) is a one dimensional autonomous system which forbids nontrivial oscillation; if $l=0$ then, with positive general kinetics, (3.8) is a decoupled system of $k$ autonomous univariate ODEs which again forbids nontrivial oscillation.

We proceed as follows. We first treat the smallest nontrivial case, namely $(k,l)=(2,1)$, which is simple enough to be fully analysed using fairly basic ideas from dynamical systems. The results of this analysis are summarised in Propostion 5.1. We then proceed as follows, ensuring that $(k,l)$ CRNs are treated after $(k-1,l)$ and $(k,l-1)$ CRNs, and treating the cases of mass action, and of physical power-law kinetics separately.

Whenever an SPPO is found in a $(k,l)$ CRN $\mathcal{R}$, we use the powerful and widely available graph-isomorphism software NAUTY [43] to identify all $(k+1,l)$ CRNs and $(k,l+1)$ CRNs $\geq\mathcal{R}$ (i.e., which include $\mathcal{R}$ as an induced subnetwork). Proposition 4.8 then tells us that these must admit SPPOs. 2. 2.

We search numerically for SPPOs in $(k,l)$ CRNs, limiting the search to those CRNs not already found to inherit SPPOs via step (1) or believed to forbid oscillation by Conjecture 5.2 below.

Via this process we obtain a lower bound on the occurrence of stable oscillation in small fully open CRNs. The methodological details are in A.

Proposition 5.1**.**

There are 14 non-isomorphic $(2,1)$ CRNs. These are, upto isomorphism, the fully open extensions of:

[TABLE]

Let $\mathcal{R}_{(k)}$ refer to the fully open extension of reaction (k), namely the CRN consisting of reaction (k) along with $X\rightleftharpoons 0\rightleftharpoons Y$.

With mass action kinetics $\mathcal{R}_{(i)}$ to $\mathcal{R}_{(xiv)}$ forbid oscillation. All but $\mathcal{R}_{(xiii)}$ have a unique equilibrium which is locally asymptotically stable and attracts all of $\mathbb{R}^{2}_{\geq 0}$. $\mathcal{R}_{(xiii)}$ either has a unique locally asymptotically stable equilibrium which attracts all of $\mathbb{R}^{2}_{\geq 0}$, or all orbits are unbounded. 2. 2.

With positive general kinetics $\mathcal{R}_{(i)}$ to $\mathcal{R}_{(xiii)}$ forbid oscillation. 3. 3.

With physical power-law kinetics or general kinetics $\mathcal{R}_{(xiv)}$ admits an SPPO.

The proof of Proposition 5.1 is fairly straightforward, but somewhat lengthy, and is in B. In order to proceed more efficiently, we make the following conjecture.

Conjecture 5.2**.**

Let $(\mathcal{R},\mathcal{K})$ be a fully open CRN with kinetics, where $\mathcal{K}$ is any scaling invariant subset of positive general kinetics (for example, $\mathcal{K}$ may be given by mass action kinetics or physical power-law kinetics). Suppose that $\Gamma$ is the stoichiometric matrix of $\mathcal{R}$ so that $\mathcal{R}$ gives rise to the family of ODEs on $\mathbb{R}^{n}_{\gg 0}$

[TABLE]

If, for all $x\in\mathbb{R}^{n}_{\gg 0}$ and all $v\in\mathcal{K}$, the Jacobian matrix $\Gamma Dv(x)$ has no purely imaginary eigenvalues, then $\mathcal{R}$ does not admit a positive periodic orbit.

A theoretical justification for Conjecture 5.2 is not attempted here, but it is not hard to believe the rather stronger claim that such families of CRNs admit oscillation if and only if they admit Hopf bifurcation (a similar conjecture is made in Section 2.2 of [26]). If Conjecture 5.2 holds, it is possible to rule out oscillation by examining, with the help of computer algebra, certain polynomials associated with $\Gamma Dv(x)$ whose positivity is sufficient to forbid purely imaginary eigenvalues. This process, which will be described in forthcoming work, is computationally much less expensive than simulating the differential equations with tens of thousands of parameter choices. No counterexamples to Conjecture 5.2 were found during a large number of numerical simulations. Note also that as the claims such as those drawn from the data in Table 1 concern lower bounds on the frequency of oscillation in CRNs, they are not invalidated if Conjecture 5.2 is false.

The results of simulations and analysis for $k=2,\ldots,4$ and $l=1,\ldots,4$ are summarised in Table 1. The table suggests, assuming that Conjecture 5.2 is true, and that large numbers of oscillatory CRNs were not missed by the numerical simulations, that the great majority of CRNs admitting stable oscillation do so as a consequence of inheritance (this becomes even more evident as we increase the number of reactions in the CRNs). As a particular example, the motif

$\scriptstyle{2}$

corresponding to the single reaction $X+Y\rightarrow 2Y$ occurs in $22\%$ of all the CRNs in Table 1, which consequently admit SPPOs with physical power-law kinetics by Propositions 5.1 and 4.8. A total of about $33\%$ of the CRNs in Table 1 were found to admit SPPOs with physical power-law kinetics and thus this single motif is responsible for about two thirds of the oscillation found under the assumption of physical power-law kinetics. Additional investigation revealed that this motif occurs in a total of about $2.52\times 10^{7}$ ($75\%$) of all $3.36\times 10^{7}$ $(2,l)$ CRNs ($l$ ranges from $1$ to $26$ by the counting arguments in A). Thus identifying small atoms of oscillation is worthwhile from a practical viewpoint, as these appear to be the source of most oscillation in CRNs.

Table 1 also highlights the importance of kinetics, and in particular how much more frequently stable oscillation occurs in small CRNs with physical power-law kinetics as compared to those with mass action kinetics. Presumably the linear or quadratic nature of at most bimolecular mass action systems significantly restricts the allowed dynamics in many cases.

As in the case of physical power-law kinetics, small oscillatory motifs account for most of the oscillation in the table found in mass action CRNs. For example, at least one of the five $(3,2)$ (presumed) mass action atoms of stable oscillation found in simulations occurs in about $5\%$ of all the CRNs in Table 1, which consequently admit SPPOs with mass action kinetics; this accounts for almost $90\%$ of the oscillation in mass action CRNs detailed in Table 1. While numerical investigations in A indicate that the lower bounds in Table 1 can be improved with additional simulation, it remains true that inheritance results applied to a few small oscillatory motifs automatically give us large numbers of oscillatory CRNs.

Not visible in the table are relationships amongst the atoms of stable oscillation. For example, the five $(3,2)$ mass action atoms of stable oscillation are the fully open extensions of: (i) $X+Z\rightarrow 2Y\rightarrow Y+Z$, (ii) $X+Z\rightarrow 2Y$, $Y+Z\rightarrow 2Z$, (iii) $X+Z\rightarrow Y,Y+Z\rightarrow 2Z$, (iv) $X+Z\rightarrow 2Y\rightarrow 2Z$; (v) $X+Z\rightarrow 0,Y+Z\rightarrow 2Z$. These correspond to the following motifs:

(i)\scriptstyle{2}$$\scriptstyle{2}$${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\scriptstyle{1}}

(ii)\scriptstyle{2}$$\scriptstyle{2}

(iii)\scriptstyle{2}$${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\scriptstyle{2}}

(iv)\scriptstyle{2}$$\scriptstyle{2}$$\scriptstyle{2}

(v)\scriptstyle{2}$${\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\scriptstyle{2}}

Representing these motifs pictorially highlights the close relationships between them. Observe that there are various subnetwork relationships between the motifs. For example, (v) is a subnetwork of (iii), but not an induced subnetwork of (iii), and hence oscillation in the fully open extension of (iii) cannot be predicted from that in the fully open extension of (v) using the theorems in this paper. There remains the possibility that there exists an inheritance result rather different from those in this paper which predicts oscillation in the fully open extension of (iii) from that in the fully open extension of (v). More generally, it seems likely that there are interesting theorems to be discovered on sufficient conditions for stable oscillation in mass action CRNs which might explain something about the structures of oscillatory motifs.

6 Conclusions

Armed with the results in this paper one can predict the occurrence of oscillation in CRNs from its occurrence in smaller CRNs. Our main conclusion is:

Any CRN built from an oscillatory CRN via a sequence of modifications of the kind described in Theorems 1 to 4 is again oscillatory.

Here “oscillatory” may be taken to mean either “which admits an NPPO” or “which admits an SPPO”, and the conclusion is valid under mild assumptions on the kinetics and for general CRNs (not necessarily fully open). We emphasised the consequence that a fully open, mass action, CRN which includes a fully open oscillatory subnetwork is itself oscillatory, illustrating how certain motifs are associated with oscillation in fully open CRNs. It was mentioned, however, that this particular conclusion does not extend to CRNs which are not fully open: such a CRN may include an oscillatory subnetwork but fail to be oscillatory. The following is a typical example:

Example 6.1**.**

Consider the following CRNs $\mathcal{R}$, $\mathcal{R}^{\prime}$ and $\mathcal{R}^{\prime\prime}$ which satisfy $\mathcal{R}\leq_{S}\mathcal{R}^{\prime}\leq_{R}\mathcal{R}^{\prime\prime}$:

[TABLE]

$\mathcal{R}$* admits an SPPO with mass action kinetics as it is just the fully open extension of motif (i) above, with the reverse of some reactions added (see Remark 4.2). On the other hand $\mathcal{R}^{\prime}$ is a weakly reversible, deficiency zero, network and, consequently, with mass action kinetics, forbids oscillation by the deficiency zero theorem [12]. Finally, by Theorem 4 applied twice to $\mathcal{R}$, $\mathcal{R}^{\prime\prime}$ admits an SPPO with mass action kinetics.*

Example 6.1 illustrates that predicting oscillation in CRNs is rather subtle: enlarging a CRN in natural ways can both destroy and create oscillation. Moreover, $\mathcal{R}^{\prime}$ and $\mathcal{R}^{\prime\prime}$ involve the same set of species and have the same stoichiometric subspace (namely, all of $\mathbb{R}^{5}$); but adding the flow reactions $0\rightleftharpoons V,\,0\rightleftharpoons W$ to $\mathcal{R}^{\prime}$ gives rise to oscillation. This corresponds to adding constant and linear terms to the differential equations describing the evolution of $\mathcal{R}^{\prime}$ with mass action kinetics.

It is highly likely that further results of the kind presented in this paper hold: following Theorems 5 and 6 in [8] we expect that modifications such as including new reactions with new species, or inserting intermediate complexes involving new species into reactions should, with mild additional hypotheses, preserve the capacity for NPPOs or SPPOs. Some oscillatory CRNs, minimal w.r.t. to the modifications described in Theorems 1 to 4 of this paper, may cease to be minimal under the improved partial order such results would bring.

There are also interesting questions on the connections between inheritance approaches as described here, and known sufficient conditions for oscillation such as those in [26, 27, 28, 29]. The families of chemical oscillators described in these papers can provide a starting point for application of the inheritance results here. It is also possible that some of the theory on families of chemical oscillators or algorithmic conditions for oscillation might suggest further inheritance results not described here. These possibilities remain to be explored.

The investigation of small, fully open, CRNs in Section 5 highlights two important points:

• •

identifying small oscillatory motifs is a worthwhile pursuit as it automatically implies oscillation in the large number of CRNs which “inherit” these motifs; and

• •

stable oscillation is much more common given larger classes of kinetics such as physical power-law kinetics as compared to mass action kinetics.

Similar studies could also be carried out for general CRNs (not necessarily fully open), using Theorems 1 to 4. The difficulty of finding oscillation in mass action CRNs by numerical experiment is evidenced by additional data in A. This data suggests that often oscillation is confined to small parameter regions, and encourages the use of more systematic algorithmic approaches to the detection of oscillation such as those in [29].

Finally, the “enumerate and simulate” methodology which provided the data in Section 5 and is described in more detail in A may also prove useful for studying the frequency of other behaviours such as chaos in CRNs [44]. Some modification to the approach may be needed to explore sets of CRNs too large to be studied exhaustively. For example, there are more than $10^{8}$ nonisomorphic $(4,5)$ CRNs, and exploring the dynamics of such large numbers numerically becomes challenging; however, it should be possible either to restrict attention to certain interesting subsets of these CRNs, such as those which are weakly reversible for example, or to explore randomly chosen CRNs from such sets in order to draw some conclusions about how often various behaviours might occur.

Acknowledgements

I would like to thank Anne Shiu and the anonymous referees for a number of helpful comments on the manuscript.

Appendix A Methodological notes

The processes of generating CRNs, and of searching numerically for oscillation, are described briefly. Further detail can be found in [40].

Generating unlabelled CRNs. Note first that two fully open CRNs are isomorphic if and only if they are isomorphic after removal from both of the flow reactions $0\rightleftharpoons X_{i}$. All at most bimolecular unlabelled, fully open CRNs with $k$ species and $l$ non-flow reactions can be generated as follows:

All at most bimolecular complexes on $k$ species are listed: there are $n_{C}(k):={k+2\choose 2}$ such complexes. 2. 2.

Irreversible reactions can be viewed as ordered pairs of distinct complexes: and consequently there are a total of $n_{C}(k)(n_{C}(k)-1)$ distinct irreversible reactions involving these complexes. 3. 3.

Excluding the reactions $0\rightarrow X_{i}$ and $X_{i}\rightarrow 0$ leaves $n_{R}(k):=n_{C}(k)(n_{C}(k)-1)-2k$ distinct non-flow reactions from which to build the CRNs. 4. 4.

All sets of $l$ distinct non-flow reactions are chosen and represented (in digraph6 format) as two-layer vertex-coloured digraphs, as described in the section Isomorphism of edge-coloured graphs of the NAUTY documentation at http://users.cecs.anu.edu.au/~bdm/nauty/nug26.pdf. There are

[TABLE]

of these CRNs. This number corresponds to the total number of $(k,l)$ CRNs with labelled species, but unlabelled reactions. 5. 5.

The NAUTY program shortg is used to canonically label and remove isomorphs from this list of CRNs, respecting the species-reaction bipartition.

Details of the enumeration methodology can be found in [40] with data at https://reaction-networks.net/networks/.

Generating $(k,l)$ CRNs which inherit oscillation. Given lists of oscillatory $(k-1,l)$ CRNs and $(k,l-1)$ CRNs (with the reactions $0\rightleftharpoons X_{i}$ removed), the following procedure generates $(k,l)$ CRNs which are oscillatory by inheritance:

Each possible new non-flow reaction is added to each oscillatory $(k,l-1)$ CRN, giving a list of oscillatory $(k,l)$ CRNs. 2. 2.

The new species $X_{k}$ is added into the reactions of each oscillatory $(k-1,l)$ CRN in every possible way which preserves bimolecularity, giving a list of oscillatory $(k,l)$ CRNs. 3. 3.

The two lists obtained in the steps above are merged and shortg is used to canonically label and remove isomorphs from the combined list. The open extensions of CRNs in the merged list are the $(k,l)$ CRNs which inherit stable oscillation, by Proposition 4.8.

Numerical simulations. The following procedure was set up to search for oscillation in CRNs which neither inherited oscillation, nor were conjectured to forbid oscillation by Conjecture 5.2. For each such CRN:

The differential equations were constructed from a combinatorial description of the CRN, along with the assumption of mass action kinetics or physical power-law kinetics. 2. 2.

A minimum of 10,000 parameter-sets were chosen randomly using uniform distributions on each parameter. The parameters chosen were rate constants and initial conditions, and additionally exponents in the case of power-law kinetics. 3. 3.

Simulations of the differential equations were carried out using RADAU5 software [45]. 4. 4.

A script was written to analyse the outputs, searching for bounded but apparently nonconvergent trajectories. Where such behaviour was identified, plots of the trajectories were output and later examined by eye to confirm that indeed oscillation had been found numerically.

Several choices were necessarily somewhat arbitrary, particularly the number of simulations, the length of simulation, and the upper and lower limits on the magnitudes of parameters. It was also assumed throughout that what appeared in a plot as a periodic orbit was indeed an SPPO.

To explore the likelihood of finding oscillation in numerical simulations, from 100 to 100,000 simulations were carried out on each of the 444 $(3,3)$ CRNs known to admit SPPOs with mass action kinetics (see Table 1). The results, shown in Table 2, reflect the fact that oscillation often occurs only in small regions of parameter space, and so can be hard to find by brute-force approaches. In larger CRNs this problem becomes even more acute, highlighting the importance of theoretical approaches, including the inheritance results in this paper.

Appendix B Proof of Proposition 5.1

*Proof of Proposition 5.1. * The 14 non-isomorphic CRNs listed in the proposition are easily enumerated with NAUTY [43] as described in A or even by eye. Under the assumption of positive general kinetics each of the CRNs $\mathcal{R}_{(i)}$ to $\mathcal{R}_{(xiv)}$ gives rise to an ODE system of the form

[TABLE]

where $a$ and $c$ are positive constants; $g_{1}$ and $g_{2}$ are positive-valued $C^{1}$ functions on $\mathbb{R}_{>0}$ with positive derivative; and $f(x,y)$, the rate of the non-flow reaction, is a positive-valued $C^{1}$ function on $\mathbb{R}^{2}_{\gg 0}$ satisfying $\frac{\partial f}{\partial x}>0$ (resp., $\frac{\partial f}{\partial y}>0$) on $\mathbb{R}^{2}_{\gg 0}$ if $X$ (resp., $Y$) occurs on the left of the non-flow reaction. $n_{X}$ and $n_{Y}$ are the net production of $X$ and $Y$ respectively in the non-flow reaction. (B.21) defines a local flow $\phi$ on $\mathbb{R}^{2}_{\gg 0}$.

(B.21) has Jacobian matrix

[TABLE]

on $\mathbb{R}^{2}_{\gg 0}$, and so

[TABLE]

The assumption of positive general kinetics implies that $g_{1}^{\prime}(x)>0$ and $g_{2}^{\prime}(y)>0$ for all positive $x$ and $y$, and so $\mathrm{det}\,J(x,y)>0$ and $\mathrm{Tr}\,J(x,y)<0$ provided $n_{X}f_{x}(x,y)\leq 0$ and $n_{Y}f_{y}(x,y)\leq 0$.

$\mathcal{R}_{(i)}$** to $\mathcal{R}_{(xii)}$.** In these cases, it is easily seen that the assumption of positive general kinetics ensures that $n_{X}f_{x}(x,y)\leq 0$ and $n_{Y}f_{y}(x,y)\leq 0$. Consequently, $\mathrm{Tr}\,J(x,y)<0$, and so each CRN forbids nontrivial periodic orbits in $\mathbb{R}^{2}_{\gg 0}$ by the Bendixson criterion (Theorem 4.1.1 of [35]). Notice that since $J(x,y)$ is everywhere Hurwitz stable on $\mathbb{R}^{2}_{\gg 0}$, all positive equilibria are locally asymptotically stable and, by the 2D Markus-Yamabe Theorem [46], existence of a positive equilibrium guarantees that it attracts all of $\mathbb{R}^{2}_{\gg 0}$. (To apply the Theorem as it is usually stated, we may first pass to logarithmic coordinates via $(x,y)\mapsto(\ln x,\ln y)$, a smooth diffeomorphism between $\mathbb{R}^{2}_{\gg 0}$ and $\mathbb{R}^{2}$.)

Under the assumption of mass action kinetics we can go further. We have $g_{1}(x)=bx$ and $g_{2}(y)=dy$ for some $b>0$ and $d>0$ and $f(x,y)=\gamma x^{\alpha}y^{\beta}$ where $\gamma>0$ and $\alpha$ and $\beta$ refer to the stoichiometries of $X$ and $Y$ on the left of the non-flow reaction. It is an easy exercise to prove that with mass action kinetics (and indeed more general assumptions), each of $\mathcal{R}_{(i)}$ to $\mathcal{R}_{(xii)}$ has a unique positive equilibrium which must, by the arguments above, attract all of $\mathbb{R}^{2}_{\gg 0}$. In fact, since each forward trajectory originating on $\partial\mathbb{R}^{2}_{\geq 0}$ immediately enters $\mathbb{R}^{2}_{\gg 0}$, the unique positive equilibrium attracts all of $\mathbb{R}^{2}_{\geq 0}$. As it is locally asymptotically stable, it is in fact globally asymptotically stable.

$\mathcal{R}_{(xiii)}$. In this case, for positive general kinetics, (B.21) is a decoupled system of two ODEs which clearly forbids periodic orbits. With mass action kinetics, the system can be written

[TABLE]

For $b>\gamma$, (B.22) has a unique, globally attracting, positive equilibrium at $(a/(b-\gamma),c/d)$, while for $b\leq\gamma$ it has only unbounded trajectories.

What remains is the only nontrivial case, namely $\mathcal{R}_{(xiv)}$:

[TABLE]

Existence of an SPPO for $\mathcal{R}_{(xiv)}$ for general kinetics and physical power-law kinetics. The two cases are dealt with at once by proving the existence of an SPPO for reaction rates which belong simultaneously to both classes of kinetics. Fixing mass-action kinetics for $Y\rightleftharpoons 0\rightleftharpoons X$, but allowing physical power-law kinetics for the non-flow reaction, the system takes the form:

[TABLE]

where $\alpha>0$ and $\beta>0$. We now prove that (B.23) has an attracting periodic orbit for appropriate choices of the parameters. We choose

[TABLE]

$k$ is a bifurcation parameter which affects two rate constants, and the choices ensure that the kinetics is polynomial and belongs both to the class of general kinetics and of physical power-law kinetics. With these parameters (B.23) becomes

[TABLE]

Note that $F(1,1;k)=0$ for all $k$ and $DF(1,1;k)$ has eigenvalues

[TABLE]

We confirm that a Hopf bifurcation occurs at $k=0$. Let $\omega=\frac{\sqrt{3}}{2}$ so that $\lambda_{\pm}(0)=\pm i\omega$. A quick calculation gives that

[TABLE]

namely $\lambda_{\pm}(k)$ cross with nonzero speed from left half to right half plane as $k$ increases through [math]. To confirm that an attracting periodic orbit indeed exists for sufficiently small $k>0$ we bring (B.24) with $k=0$ into a standard form by defining new variables $(u,v)$ via

[TABLE]

to get

[TABLE]

Here (B.25) has an equilibrium at $(0,0)$ and $f^{1}(u,v)$, $f^{2}(u,v)$ consist of the nonlinear terms in $u,v$, i.e., $f^{1}(0,0)=f^{2}(0,0)=f^{1}_{u}(0,0)=f^{1}_{v}(0,0)=f^{2}_{u}(0,0)=f^{2}_{v}(0,0)=0$. We can calculate the relevant combination of higher partial derivatives of $f^{1}$ and $f^{2}$ at $(0,0)$ to get:

[TABLE]

This computation, carried out with the help of MAXIMA [47], implies that as $k$ increases through zero we have a supercritical Hopf bifurcation and the creation of an asymptotically orbitally stable periodic orbit near to the origin in $u\mathchar 45\relax v$ space (see Section 2 of Chapter 20 in [35]). In practice, stable periodic orbits are easily found in numerical simulations of (B.23) for many choices of parameters.

Nonexistence of a periodic orbit for $\mathcal{R}_{(xiv)}$ with mass action kinetics. We now show that with mass action kinetics, $\mathcal{R}_{(xiv)}$ admits no periodic orbits, and in fact has a unique positive equilibrium which is globally asymptotically stable. Specialising to mass action kinetics, (B.23) becomes

[TABLE]

Rescaling $x\mapsto\gamma x$, $y\mapsto\gamma y$, $a\mapsto\gamma a$ and $c\mapsto\gamma c$ allows us to eliminate $\gamma$ from the system to get:

[TABLE]

(B.26) defines a local semiflow $\phi$ on $\mathbb{R}^{2}_{\geq 0}$. Let $G(x,y):=\frac{1}{y}F(x,y)$. Then

[TABLE]

As this holds everywhere on $\mathbb{R}^{2}_{\gg 0}$, (B.26) has no nontrivial periodic orbits in $\mathbb{R}^{2}_{\gg 0}$ by the Bendixson-Dulac criterion (Theorem 4.1.2 of [35]).

Next, we show that (B.26) has a unique linearly stable equilibrium which attracts all of $\mathbb{R}^{2}_{\geq 0}$. As each point of $\partial\mathbb{R}^{2}_{\geq 0}$ is a start point of $\phi$, $\phi$ has no $\omega$-limit points on $\partial\mathbb{R}^{2}_{\geq 0}$ (Theorem 5.6 in [41] for example). Note that the triangle $T$ defined by

[TABLE]

is a global attractor for (B.26) as

[TABLE]

Thus all orbits of $\phi$ are bounded and $\phi$ is in fact a semiflow on $\mathbb{R}^{2}_{\geq 0}$. Moreover each point of $\partial T$ is a start point of $\left.\phi\right|_{T}$, and so $\left.\phi\right|_{T}$ (and hence $\phi$) has no $\omega$-limit points on $\partial T$.

There is a unique positive solution $(x_{0},y_{0})$ to $F(x,y)=0$ given by

[TABLE]

As $\theta>a+c-bd$, we see that $x_{0}<d$. The Jacobian matrix of the system is

[TABLE]

At $(x_{0},y_{0})$, as $x_{0}<d$,

[TABLE]

and so $(x_{0},y_{0})$ is linearly stable. Certainly, $(x_{0},y_{0})$ has no homoclinic orbits. Since all trajectories of (B.26) enter $\mathrm{int}\,T$ and $\phi$ has no nontrivial periodic or homoclinic orbits in $\mathrm{int}\,T$, by the Poincaré-Bendixson theorem (Theorem 9.0.6 in [35] for example), $(x_{0},y_{0})$ is a global attractor of (B.26) and hence is globally asymptotically stable. $\square$