Building oscillatory chemical reaction networks by adding reversible reactions

TL;DR

This paper demonstrates that adding reversible reactions to a chemical reaction network that admits oscillation can preserve that oscillatory behavior under certain conditions, extending understanding of network inheritance properties.

Contribution

It establishes conditions under which oscillation is inherited when reversible reactions are added to a CRN, advancing the theory of oscillatory network construction.

Findings

Oscillation inheritance is possible with reversible reactions under mild conditions.

The results extend previous inheritance theorems for multistationarity.

Provides criteria for designing larger oscillatory networks from smaller ones.

Abstract

We show that if a chemical reaction network (CRN) admits nondegenerate (resp., linearly stable) oscillation, and we add new reversible reactions involving new species to this CRN, then the new CRN so created also admits nondegenerate (resp., linearly stable) oscillation provided certain mild and easily checked conditions are met. This claim that the larger CRN "inherits" oscillation from the smaller one, provided it is built from the smaller CRN in an appropriate way, follows an analogous result involving multistationarity. It also adds to a number of prior results on the inheritance of oscillation; these collectively often allow us to determine the capacity of a given network for oscillation based on an analysis of its subnetworks.

Equations

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

$$a_i\cdot X+b_i\cdot Y\rightleftharpoons a_i^{\prime }\cdot X+b_i^{\prime }\cdot Y,(i=1,\dots ,m).$$

$$\left(\begin{array}[]{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}[]{cc}\Gamma&\alpha\\ 0&\beta\end{array}\right)\left(\begin{array}[]{c}v(x)\\ q(x,y)\end{array}\right)\,$$

$${\mathcal{R}}_1$$

$${\mathcal{R}}_2$$

$$\left(\begin{array}[]{rr}1&-1\\ 1&2\\ 0&1\end{array}\right)$$

$${\mathbb{R}}_{\gg 0}^n:=\{x\in {\mathbb{R}}^n:x_i>0(i=1,\dots ,n)\}.$$

$${\mathbb{R}}_{\ge 0}^n:=\{x\in {\mathbb{R}}^n:x_i\ge 0(i=1,\dots ,n)\}.$$

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

$$D{\Psi }_t(h(x_0))=Dh({\Phi }_t(x_0))D{\Phi }_t(x_0)[D{h}^{-1}(h(x_0))].$$

$$D{\Psi }_T(y_0)=[Dh(x_0)]D{\Phi }_T(x_0)[Dh(x_0){]}^{-1}.$$

$$\{(x,y)\in U\times V:F(x,y)=0\}=\{(x,\varphi (x)):x\in U\},$$

$$\{(x,y)\in U\times B_t:F(x,y)=0\}=\{(x,\varphi (x)):x\in U\}.$$

$$\{(x,y)\in B_{s_x}(x)\times B_{t_x}:F(x,y)=0\}=\{(x,{\varphi }_x(x)):x\in B_{s_x}(x)\}.$$

$$\hat{\theta}(z,\eta):=\left\{\begin{array}[]{rcl}\theta(z,\eta)/\eta&(\eta\neq 0)\\ D_{\eta}\theta(z,0)&(\eta=0)\end{array}\right.$$

$$|\hat{\theta }(z,\eta )-\hat{\theta }(z_0,0)|=|\hat{\theta }(z,0)-\hat{\theta }(z_0,0)|=|D_{\eta }\theta (z,0)-D_{\eta }\theta (z_0,0)|$$

$$|\hat{\theta }(z,\eta )-\hat{\theta }(z_0,0)|$$

$$|\hat{\theta }(z_0,\eta )-\hat{\theta }(z_0,0)|=|\theta (z_0,\eta )/\eta -D_{\eta }\theta (z_0,0)|=|R(z_0,\eta )|.$$

$$|\theta (z,\eta )-\theta (z_0,\eta )|=\left|{\int }_0^{\eta }(D_{\eta }\theta (z,\sigma )-D_{\eta }\theta (z_0,\sigma ))\mathrm{d}\sigma \right|\le {\int }_0^{\eta }\left|D_{\eta }\theta (z,\sigma )-D_{\eta }\theta (z_0,\sigma )\right|\mathrm{d}\sigma .$$

$$\frac{|\theta (z,\eta )-\theta (z_0,\eta )|}{\eta }\le \frac{1}{\eta }{\int }_0^{\eta }K|z-z_0|\mathrm{d}\sigma \le K|z-z_0|.$$

$$|\hat{\theta }(z,\eta )-\hat{\theta }(z_0,0)|\le K|z-z_0|+|R(z_0,\eta )|.$$

$$\Gamma^{\prime}:=\left(\begin{array}[]{cc}\Gamma&\alpha\\ 0&\beta\end{array}\right)\,.$$

$$\phi_{0}\left(\begin{array}[]{c}x\\ y\end{array}\right)=\left(\begin{array}[]{c}x-\alpha\hat{\beta}^{-1}\hat{y}\\ y\end{array}\right),\quad\mbox{so that}\quad\phi_{0}^{-1}\left(\begin{array}[]{c}z\\ y\end{array}\right)=\left(\begin{array}[]{c}z+\alpha\hat{\beta}^{-1}\hat{y}\\ y\end{array}\right)\,.$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{im}}\,\left(\begin{array}[]{cc}\Gamma&\alpha\\ 0&\beta\end{array}\right)\quad\mbox{in $(x,y)$-space to cosets of }{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{im}}\,\left(\begin{array}[]{cc}\Gamma&0\\ 0&\beta\end{array}\right)\quad\mbox{in $(z,y)$-space.}$$

$$\begin{array}[]{rcl}\dot{z}&=&\Gamma v(z+\alpha\hat{\beta}^{-1}\hat{y})\,,\\ \dot{y}&=&\beta q(z+\alpha\hat{\beta}^{-1}\hat{y},y)\,,\end{array}$$

$$\mathcal{H}:=\{(z,y)\in {\mathbb{R}}^n\times {\mathbb{R}}^{m+k}:{\delta }^{\mathrm{t}}\hat{y}+\bar{y}=1\}.$$

$${\varphi }_0^{-1}\circ {\varphi }_1^{-1}(S)=S^{\ast }.$$

$$W_{+}:=\{(z,\hat{y})\in {\mathbb{R}}^n\times {\mathbb{R}}^m:{\varphi }_0^{-1}\circ {\varphi }_1^{-1}(z,\hat{y})\in {\mathbb{R}}_{\gg 0}^n\times {\mathbb{R}}_{\gg 0}^{m+k}\},$$

$$S_{+}:=\{(z,\hat{y})\in S:{\varphi }_0^{-1}\circ {\varphi }_1^{-1}(z,\hat{y})\in {\mathbb{R}}_{\gg 0}^n\times {\mathbb{R}}_{\gg 0}^{m+k}\},$$

$$\begin{array}[]{rcl}\dot{z}&=&\Gamma v(z+\alpha\hat{\beta}^{-1}\hat{y})\,,\\ \dot{\hat{y}}&=&\hat{\beta}q_{*}(z,\hat{y})\,.\end{array}$$

Full text

1 Introduction and statement of the main result

A mathematically interesting and practically important question is when we can infer some dynamical behaviour in a network model based on knowledge that this behaviour occurs in some model of a subnetwork. Results in this area focussed on chemical reaction networks (CRNs), for example in [1, 2, 3, 4, 5, 6], have illustrated that this is a rather subtle question. Some intuitively plausible claims turn out to be hard to prove, or to be false. This paper contributes to this literature. The dynamical behaviour of interest here is oscillation, and the main result to be proved (Theorem 1 below) was conjectured to hold in the concluding sections of [6]; however the proof turned out to be somewhat harder than expected.

The work here can be motivated by the following natural questions about sequestration or inhibition. To make the arguments concrete, consider any standard biochemical model exhibiting stable oscillations, such as the classical two-pool calcium model of Goldbeter et al [7] or models of glycolytic oscillations [8]. Now suppose we introduce a reversible “sequestration” process where one of the existing chemical species can bind reversibly to some species, either already in the model or new, to create a new inactive complex. Equivalently (from a mathematical point of view), we introduce an allosteric inhibitor which binds reversibly to some existing species converting it to an inactive form. Could such a process, with certainty, destroy stable oscillation? The answer is no – maintaining the kinetics of the original reactions, the model will still stably oscillate if we choose mass action kinetics and appropriate rate constants for the inhibition/sequestration process. This claim is intuitively plausible. The idea is that provided (i) the binding and unbinding rates are fast so the new reaction is trying hard to equilibriate, and (ii) the ratio of rates is chosen so the bound, inactive, form has low equilibrium concentration, then the dynamics of the enlarged model projected onto the original species space is “close” to that of the original model, and stable oscillation should survive by perturbation arguments.

Theorem 1 of this paper is much more general than this motivating discussion suggests, but is guided by the desire to find under what circumstances similar arguments work. The Theorem states that we can build a new oscillatory netork by adding into an existing oscillatory network new reversible reactions, but with a caveat: some new species must figure nondegenerately in the new reactions. This condition is made precise later, but can easily be illustrated in the special case of a single added reaction, when it becomes: “there must be a net change in at least one new species in the added reaction”. Example 1.1 provides a simple illustration of the result in this special case. Meanings of the terms, and assumptions about reaction kinetics, will follow later.

Example 1.1**.**

Suppose that we have a CRN $\mathcal{R}$ on chemical species $X_{1},\ldots,X_{n}$ admitting linearly stable oscillation. We now build a larger CRN $\mathcal{R}^{\prime}$ by adding to $\mathcal{R}$ a reaction $\mathcal{R}_{0}$ involving some new species. Then, for example:

• (i)

If $\mathcal{R}_{0}$ is $X_{1}+X_{n+1}\rightleftharpoons 2X_{n+1}$ then $\mathcal{R}^{\prime}$ also admits linearly stable oscillation: there is a net change in the new species $X_{n+1}$ in the added reaction.

• (ii)

If $\mathcal{R}_{0}$ is $X_{1}+X_{n+1}\rightleftharpoons X_{n+1}+X_{n+2}$ then $\mathcal{R}^{\prime}$ admits linearly stable oscillation: there is a net change in the new species $X_{n+2}$ in the added reaction.

• (iii)

If $\mathcal{R}_{0}$ is $X_{1}+X_{n+1}\rightleftharpoons X_{n+1}$ then we cannot conclude from Theorem 1 that $\mathcal{R}^{\prime}$ admits oscillation: there is a new species involved, but the added reaction does not cause any net change in this new species.

It is little surprise that perturbation theory (both regular and singular) forms the backbone of the proof of Theorem 1. The challenge which takes up the majority of our effort here is to recast the basic problem in a form amenable to geometric singular perturbation theory approaches. This requires effort since, firstly, adding reactions into a system alters the dynamics of both the new and the original model species and, secondly, various mathematical necessary conditions for the application of geometric singular perturbation theory are easily violated.

Returning to the broader context of the work, this paper contributes to the literature on oscillation in CRNs. This literature has a considerable history and includes theoretical, numerical, and algorithmic work, focussed on both ruling out oscillation, and finding oscillation or bifurcations which give rise to oscillation. It would be hard to compile a complete list of papers about, or with important implications for, oscillations in CRNs. Instead, the following is a small sample, illustrating both numerical and applied work, and some key strands of classical and modern theory: [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. Inheritance results of the kind here provide an important theoretical tool for guaranteeing the occurrence of oscillation in CRNs without resorting to numerical investigation.

Presenting a detailed introduction to the mathematical theory of CRNs can often take up a considerable chunk of a paper, and a minimal approach is adopted here: we intersperse key definitions into the text without much discussion. The reader is referred to some of the papers referenced above and to [25] for a more thorough background using notation and conventions close to those adopted here.

We now turn to statement of the main result. Consider a chemical reaction network $\mathcal{R}$ involving $n$ species $X_{1},\ldots,X_{n}$ collectively termed $X$. Let species $X_{i}$ have concentration $x_{i}$ ($i=1,\ldots,n$). We are interested in positive concentrations, namely $x:=(x_{1},\ldots,x_{n})^{\mathrm{t}}\in\mathbb{R}^{n}_{\gg 0}:=\{x\in\mathbb{R}^{n}\colon x_{i}>0\,\,\mbox{for all}\,\,i\}$. Suppose that we have $r_{0}$ chemical reactions involving $X$, and that the $i$th reaction has reaction vector $\Gamma_{i}$ and reaction rate $v_{i}(x)$. We assume only that $v_{i}\colon\mathbb{R}^{n}_{\gg 0}\to\mathbb{R}$ is $C^{2}$; other than this condition, the kinetics of the reactions is arbitrary. $v(x):=(v_{1}(x),\ldots,v_{r_{0}}(x))^{\mathrm{t}}$ is termed the reaction rate vector for the system, and $\Gamma:=[\Gamma_{1}|\cdots|\Gamma_{r_{0}}]$ is termed the stoichiometric matrix of the system. Then the following system of ODEs on $\mathbb{R}^{n}_{\gg 0}$ describes the evolution of the concentration vector $x$.

[TABLE]

Note that the RHS of (1) belongs to $\mathrm{im}\,\Gamma$, a linear subspace of $\mathbb{R}^{n}$ termed the stoichiometric subspace of the system, and consequently cosets of $\mathrm{im}\,\Gamma$ are invariant under the local flow defined by (1) on $\mathbb{R}^{n}_{\gg 0}$. The intersection of these cosets of $\mathrm{im}\,\Gamma$ with $\mathbb{R}^{n}_{\gg 0}$ are termed the positive stoichiometry classes of the system.

Now let $m\geq 1$ and $k\geq 0$ be integers and suppose that we add to the system $m$ new reversible reactions involving $m+k$ new species $Y_{1},\ldots,Y_{m+k}$, collectively termed $Y$. The new CRN obtained from $\mathcal{R}$ by adding in the new reversible reactions will be termed $\mathcal{R}^{\prime}$. In order to state a nondegeneracy condition on the added reactions we need to describe these reactions, and for this we introduce some notation.

Given a list of species, say $X_{1},\ldots,X_{n}$, and a nonnegative integer vector of the same length, say $c=(c_{1},\ldots,c_{n})^{\mathrm{t}}$, we write $c\cdot X$ for the formal sum (i.e., complex in CRN terminology) $c_{1}X_{1}+c_{2}X_{2}+\cdots+c_{n}X_{n}$. We simply write “[math]” for the zero complex $0X_{1}+\cdots+0X_{n}$. Using this notation, let the new reactions be:

[TABLE]

Here $a_{i},a_{i}^{\prime},b_{i}$ and $b_{i}^{\prime}$ are nonnegative integer vectors of length $n$, $n$, $m+k$ and $m+k$ respectively. Define $a=(a_{1}|a_{2}|\cdots|a_{m})\in\mathbb{R}^{n\times m}$, with $a^{\prime}\in\mathbb{R}^{n\times m}$, $b\in\mathbb{R}^{(m+k)\times m}$ and $b^{\prime}\in\mathbb{R}^{(m+k)\times m}$ defined similarly. Define $\alpha=a^{\prime}-a\in\mathbb{R}^{n\times m}$ and $\beta=b^{\prime}-b\in\mathbb{R}^{(m+k)\times m}$. $\alpha$ records the net stoichiometric changes of the old species $X$ in the added reactions. $\beta$ records the net stoichiometric changes of the new species $Y$ in the added reactions and occurs in a nondegeneracy condition in Theorem 1 below. Let $y_{i}$ denote the concentration of $Y_{i}$ ($i=1,\ldots,m+k$), and define $y:=(y_{1},\ldots,y_{m+k})^{\mathrm{t}}$. If the new reactions have reaction rate vector $q\colon\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0}\to\mathbb{R}^{m}$, then $\mathcal{R}^{\prime}$ evolves according to:

[TABLE]

on $\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0}$. We are now ready to state our main result, although some of the definitions to make it precise will follow.

Theorem 1**.**

*Suppose the CRN $\mathcal{R}$ with evolution described by (1) has a nondegenerate (resp., linearly stable) positive periodic orbit. Let $\mathcal{R}^{\prime}$ be the CRN with evolution described by (3), obtained by adding in the reactions of (2) to $\mathcal{R}$. Suppose (i) $\beta$ has rank equal to $m$, its number of columns, and (ii) the added reactions are given mass action kinetics. Then rate constants can be chosen for the added reactions in such a way that $\mathcal{R}^{\prime}$ has a nondegenerate (resp., linearly stable) positive periodic orbit. *

By a positive periodic orbit, we mean a periodic orbit that lies in the (strictly) positive orthant. By a nondegenerate periodic orbit, we mean one that is hyperbolic relative to its stoichiometry class. Linearly stable is also taken to mean linearly stable relative to its stoichiometry class. Precise statement of these latter conditions is deferred to Section 3.

Theorem 1 is exactly analogous (including the condition that $\beta$ has rank $m$) to Theorem 5 in [5], replacing “multiple positive nondegenerate (resp., linearly stable) equilibria” with “a nondegenerate (resp., linearly stable) positive periodic orbit”. The proof draws heavily both on techniques developed in [5], and on singular perturbation theory approaches which formed the basis for the proof of Theorem 4 in [6]. It is hoped that the proof of Theorem 1 provides a template for, or at least insight towards, the proof of further inheritance results on CRNs.

2 An example

The proof of Theorem 1 is constructive: it not only tells us about inheritance of oscillation, but also gives information about parameter regions at which oscillation occurs. Before the proof, we present an example illustrating the result, including how to choose parameter values at which we can observe inherited oscillation. We will re-examine in detail this example in Section 5 after presenting the proof of Theorem 1 and use it to help elucidate various ideas and quantities in the proof. Here we use it only to outline the conclusions.

In [6], the following was presented as an example of a CRN which admits stable oscillation with mass action kinetics:

[TABLE]

This is an example of a so-called fully open CRN on three species $X,Y$ and $Z$, as it includes the inflow-outflow reactions $0\rightleftharpoons X$, $0\rightleftharpoons Y$ and $0\rightleftharpoons Z$. The reactions are labelled with their rate constants. In numerical simulations we easily find parameter regions where the CRN admits a stable periodic orbit. For example, setting $k_{1}=4,\,k_{2}=3,\,k_{3}=0.2,\,k_{4}=2,\,k_{5}=0.3,\,k_{6}=2.5,\,k_{7}=2.5$ and $k_{8}=0.2$, and choosing initial conditions $X_{0}=Y_{0}=Z_{0}=1$ we find the system settles, after initial transients, onto the periodic orbit shown in Figure 1. We assume that the system does indeed have a positive, linearly stable periodic orbit at these parameter values.

Now suppose that we add in two further reversible reactions $Y\rightleftharpoons U+V$ and $U+X\rightleftharpoons 2V+W$ involving three new species $U,V$ and $W$ to obtain the system

[TABLE]

The matrix $\beta$ representing the net stoichiometric changes of the new species in the added reactions is, in this case

[TABLE]

which clearly has rank $2$. Consequently, the nondegeneracy condition in Theorem 1 is satisfied, and the theorem tells us that $\mathcal{R}_{2}$ admits stable oscillation with mass action kinetics.

The proof of Theorem 1 also tells us how to set rate constants for the added reactions to obtain oscillation in $\mathcal{R}_{2}$. We define two parameters $\epsilon$ and $\eta$, set $k_{9}=\epsilon^{-1}$, $k_{10}=\epsilon^{-1}\eta^{-2}$, $k_{11}=\epsilon^{-1}\eta^{-1}$ and $k_{12}=\epsilon^{-1}\eta^{-2}$, and leave all other rate constants as before; then the proof of Theorem 1 tells us that by choosing and fixing $\eta>0$ sufficiently small, and subsequently choosing and fixing $\epsilon>0$ sufficiently small, we can ensure that $\mathcal{R}_{2}$ has a positive periodic orbit which is linearly stable relative to its stoichiometry class. Moreover, with these choices, variation in the values of $U$, $V$ and $W$ on the periodic orbit will be small; the values of $U$ and $V$ on the periodic orbit will be small; the values of $W$ on the periodic orbit can be controlled by the choice of initial data; and the values of $X$, $Y$ and $Z$ on the periodic orbit will be close to their original values in the absence of the added reactions. Some plots of the periodic orbit (omitting transient behaviour) in the case $\epsilon=\eta=0.2$, and with initial values of the new variables $U_{0}=V_{0}=0$, $W_{0}=1$ are shown in Figure 2. Note that $\mathcal{R}_{2}$ now has a conserved quantity $3W+U-V$.

3 Technical preliminaries

Before presenting the proof of Theorem 1 we need some notational and mathematical preliminaries from analysis and the theory of differential equations.

3.1 Basic notation, conventions, and definitions

We draw heavily on [5] and [6] here.

Definition 3.1** (Positive subsets of $\mathbb{R}^{n}$).**

We refer to a subset of $\mathbb{R}^{n}$ as positive if it is a subset of

[TABLE]

Given $x\in\mathbb{R}^{n}$, we write $x\gg 0$ to mean $x\in\mathbb{R}^{n}_{\gg 0}$. We also define

[TABLE]

Notation 3.2** (Vector of ones, identity matrix).**

$\mathbf{1}$* denotes a vector of ones whose length is inferred from the context. If $\eta$ is any real constant, then $\bm{\eta}$ denotes a vector whose entries are all $\eta$ and whose length is inferred from the context. $I_{n}$ refers to the $n\times n$ identity matrix. *

Definition 3.3** (Empty vectors and matrices).**

*In order to simplify some arguments, we formally allow vectors and matrices to be empty. An empty matrix is one with zero rows, zero columns, or both: when we define a matrix to be $n\times m$, one or both of $n$ or $m$ is allowed to be zero. An empty vector, for example, is a $0\times 1$ matrix. Empty vectors and matrices obey the following natural rules. (i) If $A$ is an $n\times m$ matrix and $B$ is an $m\times k$ matrix, then $AB$ is defined and is an $n\times k$ matrix, even if some of $n,m$ or $k$ are zero. If $m=0$, but $n$ and $k$ are nonzero, then $AB$ is defined to be the $n\times k$ zero matrix. (ii) Any equation, inequality, or claim involving empty vectors or matrices is vacuously satisfied (provided that it makes sense, dimensionally). (iii) Given an empty vector $y$ and a $k\times 0$ matrix $A$, $y^{A}$ is defined to be $\mathbf{1}$, a vector of ones of length $k$. *

Notation 3.4** (Sum of a point and a set).**

*Given a point $x_{0}\in\mathbb{R}^{n}$ and a set $A\subseteq\mathbb{R}^{n}$, $x_{0}+A$ means, naturally, the following subset of $\mathbb{R}^{n}$: $\{x_{0}+y\colon y\in A\}$. *

Notation 3.5** (Open ball in $\mathbb{R}^{n}$).**

*For any $x\in\mathbb{R}^{n}$ and $r>0$, let $B_{r}(x)$, be the open ball in $\mathbb{R}^{n}$ of radius $r$ with centre $x$, namely $B_{r}(x):=\{y\in\mathbb{R}^{n}\colon|y-x|<r\}$. If the argument $x$ is omitted, it is taken to be zero. The dimension $n$ is to be inferred from the context. *

Notation 3.6** (Hausdorff distance).**

*Given two nonempty sets $A$ and $B$ in $\mathbb{R}^{n}$ with the Euclidean metric, $d_{\mathrm{H}}(A,B)$ denotes the Hausdorff distance between $A$ and $B$. *

Notation 3.7** (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 3.8** (Entrywise product and entrywise functions).**

*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}$. When we apply functions such as $\ln(\cdot)$ and $\exp(\cdot)$ with a vector or matrix argument, we mean entrywise application. Similarly, if $x=(x_{1},\ldots,x_{n})^{\mathrm{t}}$ and $y=(y_{1},\ldots,y_{n})^{\mathrm{t}}$, then $x/y$ means $(x_{1}/y_{1},x_{2}/y_{2},\ldots,x_{n}/y_{n})^{\mathrm{t}}$. *

Definition 3.9** (Mass action kinetics, rate constants).**

*A chemical reaction $a\cdot X\rightarrow b\cdot X$ is said to have mass action kinetics if the rate of reaction is $kx^{a^{\mathrm{t}}}$ for some positive constant $k$ termed the rate constant of the chemical reaction. *

Notation 3.10** (Preimages of sets).**

*Given a function $f\colon A\to B$, and any $\mathcal{B}\subseteq B$, $f^{-1}(\mathcal{B})$ refers, naturally, to $\{a\in A\colon f(a)\in\mathcal{B}\}$. *

Remark 3.11** (Differentiability of functions).**

When we refer to a function $f$ as being $C^{r}$ on some set $U\subseteq\mathbb{R}^{n}$, not necessarily open, we mean that there exists a function $\hat{f}$ defined and $C^{r}$ on an open set $V\subseteq\mathbb{R}^{n}$ containing $U$ and such that $\hat{f}$ coincides with $f$ on $U$.

Notation 3.12** (Derivatives of functions).**

*Given a differentiable function $f\colon U\subseteq\mathbb{R}^{n}\to\mathbb{R}^{m}$, $Df$ refers both to the derivative of $f$ and also its matrix representation where the bases on $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$ are the standard bases or are to be inferred from the context. Given a set of positive integers $n_{1},\ldots,n_{m},k$, variables $x_{i}\in\mathbb{R}^{n_{i}}$ ($i=1,\ldots,m$) and a differentiable function $f\colon\mathbb{R}^{n_{1}}\times\cdots\times\mathbb{R}^{n_{m}}\to\mathbb{R}^{k}$, $D_{x_{i}}f$ refers to the derivative of $f$ w.r.t. the variable $x_{i}$ and also its matrix representation. We may also write $D_{i}f$ for the derivative of a function $f$ w.r.t. its $i$th argument, or the matrix representation of this derivative. *

The following three examples, reproduced or adapted from [6], demonstrate how entrywise and monomial notation greatly abbreviate otherwise lengthy calculations.

Example 3.13** (Rules of exponentiation).**

*Let $x,y\in\mathbb{R}^{m}_{\gg 0}$, $A,B\in\mathbb{R}^{n\times m}$ and $C\in\mathbb{R}^{k\times n}$. Let $O$ refer to the $n\times m$ matrix of zeros. Then (i) $x^{O}=\mathbf{1}$, (ii) $x^{A+B}=x^{A}\circ x^{B}$, (iii) $x^{A}\circ y^{A}=(x\circ y)^{A}$ and (iv) $(x^{A})^{C}=x^{CA}$. *

Example 3.14** (Logarithm of monomials).**

*Suppose $x\in\mathbb{R}^{m}_{\gg 0}$, $y_{i}\in\mathbb{R}^{n_{i}}_{\gg 0}$ ($i=1,\ldots,k$) and $A_{i}\in\mathbb{R}^{m\times n_{i}}$ ($i=1,\ldots,k$). If $w=x\circ y_{1}^{A_{1}}\circ\cdots\circ y_{k}^{A_{k}}$, then $\ln w=\ln x+A_{1}\ln y_{1}+\cdots+A_{k}\ln y_{k}$. *

Example 3.15** (Differentiation of monomials).**

*Suppose $k\in\mathbb{R}^{m}_{\gg 0}$, $x\in\mathbb{R}^{n}_{\gg 0}$, $A\in\mathbb{R}^{m\times n}$. Let $w\colon\mathbb{R}^{n}_{\gg 0}\to\mathbb{R}^{m}_{\gg 0}$ be defined by $w(x):=k\circ x^{A}$. Then $Dw=\mathrm{diag}(w)\,A\,\mathrm{diag}(\mathbf{1}/x)$. *

3.2 Periodic orbits

We need a number of standard results on periodic orbits and Floquet theory largely as described in Section 2 of [6]. We summarise these here, but the reader is referred to [6] and the original sources ([26] for example) for more detail.

Consider some system of ODEs $\dot{x}=F(x)$ on an open set $U\subseteq\mathbb{R}^{n}$, satisfying conditions for existence and uniqueness of solutions and hence defining a local flow $\Phi_{t}\colon U\to U$. For each $x\in U$, $t$ belongs to an open interval including [math] which in general depends on $x$, and $\Phi_{t}(x)$ is the point that initial condition $x$ “reaches” at time $t$. Given some $T>0$ the orbit of a nontrivial $T$-periodic solution of the ODE system is termed a periodic orbit of the system. Associated with any such periodic orbit are its Floquet multipliers (or characteristic multipliers), namely eigenvalues of $D\Phi_{T}(x_{0})$ where $x_{0}$ is any point on the periodic orbit, and $D\Phi_{T}(x_{0})$ refers to the derivative of $\Phi_{T}$ w.r.t. $x$ evaluated at $x_{0}$. Here $D\Phi_{t}(x_{0})$ can be regarded as the fundamental matrix solution of the $T$-periodic variational equation $\dot{y}=DF(\Phi_{t}(x_{0}))y$ satisfying $y_{0}=I_{n}$. The choice of $x_{0}$ does not affect the Floquet multipliers.

Any periodic orbit always has one Floquet multiplier of $1$ corresponding to the direction tangential to the periodic orbit; the remaining Floquet multipliers are termed the nontrivial Floquet multipliers of the periodic orbit. If none of the nontrivial Floquet multipliers lie on the unit circle in the complex plane, then the periodic orbit is hyperbolic, and, in our terminology here, nondegenerate. If, further, all of the nontrivial Floquet multipliers lie inside the unit circle in the complex plane, then the periodic orbit is linearly stable and attracts a neighbourhood of itself.

Given any nondegenerate (resp., linearly stable) periodic orbit, by regular perturbation theory arguments involving, for example, the construction of a Poincaré map, a nearby nondegenerate (resp., linearly stable) periodic orbit (with nearby period) exists for all ODEs on $U$ $C^{1}$-close to $\dot{x}=F(x)$. More precisely, we have the following result, which appears as Lemma 2.1 in [6]. A proof can be found in Section IV of [27].

Lemma 3.16**.**

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

[TABLE]

*Suppose that $\dot{x}=F(x,0)$ has a nontrivial hyperbolic (resp., linearly stable) $T$-periodic orbit $\mathcal{O}\subseteq U$. Then there exists $\epsilon_{0}\in(0,\epsilon^{\prime}]$ s.t. for $\epsilon\in(-\epsilon_{0},\epsilon_{0})$ (4) 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$. *

An important basic observation that we will frequently need is that the Floquet multipliers of a periodic orbit are invariant under $C^{1}$-diffeomorphisms, and hence the notions of “nondegeneracy” and “linear stability” of a periodic orbit are invariant under $C^{1}$-diffeomorphisms. To be more precise:

Lemma 3.17** (Invariance of Floquet multipliers).**

*Suppose we have a $C^{1}$ local flow $\Phi_{t}$ on an open set $U\subseteq\mathbb{R}^{n}$ with a periodic orbit $\mathcal{O}\subseteq U$. Let $h\colon U\to h(U):=V\subseteq\mathbb{R}^{n}$ be a $C^{1}$ diffeomorphism, so that we have a new local flow $\Psi_{t}=h\circ\Phi_{t}\circ h^{-1}$ on $V$, and let $\mathcal{O}^{\prime}:=h(\mathcal{O})$ be the corresponding periodic orbit of $\Psi_{t}$. Then $\mathcal{O}$ and $\mathcal{O}^{\prime}$ have the same Floquet multipliers. *

Proof.

Given any point $x_{0}\in\mathcal{O}$ and the corresponding point $y_{0}=h(x_{0})\in\mathcal{O}^{\prime}$, $D\Phi_{t}(x_{0})$ and $D\Psi_{t}(y_{0})$ are linear operators from $T_{x_{0}}U\cong\mathbb{R}^{n}$ to $T_{\Phi_{t}(x_{0})}U\cong\mathbb{R}^{n}$, and $T_{y_{0}}V\cong\mathbb{R}^{n}$ to $T_{\Psi_{t}(y_{0})}V\cong\mathbb{R}^{n}$ respectively. They are defined for all $t\in\mathbb{R}$ (since $\Phi_{t}(x_{0})$ is periodic in $t$), and clearly satisfy $D\Phi_{0}=\mathrm{id}$ and $D\Psi_{0}=\mathrm{id}$. Applying the chain rule to $\Psi_{t}=h\circ\Phi_{t}\circ h^{-1}$ gives

[TABLE]

Setting $t=T$, and observing that $h(x_{0})=y$, $\Phi_{T}(x_{0})=x_{0}$ and $Dh^{-1}(h(x_{0}))=[Dh(x_{0})]^{-1}$ gives:

[TABLE]

Now the Floquet multipliers of $\mathcal{O}$ are the eigenvalues of $D\Phi_{T}(x_{0})$ and the Floquet multipliers of $\mathcal{O}^{\prime}$ are the eigenvalues of $D\Psi_{T}(y_{0})$. Since the final equation shows that $D\Phi_{T}(x_{0})$ and $D\Psi_{T}(y_{0})$ are similar, the two sets of Floquet multipliers are equal.

Remark 3.18** (Floquet multipliers relative to a given set).**

Let $\Phi_{t}$ be a local flow on an open set $U^{\prime}\subseteq\mathbb{R}^{n}$, let $U\subseteq U^{\prime}$ be locally invariant under $\Phi_{t}$, and let $\mathcal{O}\subseteq U$ be a periodic orbit of $\Phi_{t}$. Let $V\subseteq\mathbb{R}^{m}$ be open and suppose that $h\colon U\to V$ is a $C^{1}$-diffeomorphism. In the light of Lemma 3.17 it makes sense to refer to the Floquet multipliers of $\mathcal{O}$ relative to $U$. We mean the Floquet multipliers of $h(\mathcal{O})$ for the derived flow $\Psi_{t}=h\circ\Phi_{t}\circ h^{-1}$ on $V$, which, by Lemma 3.17, do not depend on $V$ or $h$.

Nondegenerate/linearly stable periodic orbits for a CRN. Suppose now that we have a chemical reaction network with stoichiometric matrix $\Gamma$ defining a system of ODEs $\dot{x}=\Gamma v(x)$ as in (1). Since cosets of $\mathrm{im}\,\Gamma$ are invariant under the local flow defined by such a system, any periodic orbit must belong to one of these sets. If $\Gamma$ has rank $r$, less than its number of rows $n$, then it is easily seen that no periodic orbit can be nondegenerate or linearly stable in the sense described above since any periodic orbit has $n-r$ nontrivial Floquet multipliers with value $1$ corresponding to directions transverse to the coset of $\mathrm{im}\,\Gamma$ on which it lies. In this situation, we follow [6] and overload the terms nondegenerate and linearly stable as follows. We say that a periodic orbit $\mathcal{O}$ is nondegenerate (resp., linearly stable) if it has $r-1$ Floquet multipliers which are disjoint from (resp., inside) the unit circle, or equivalently if it is nondegenerate (resp., linearly stable) relative to the coset of $\mathrm{im}\,\Gamma$ on which it lies in the sense of Remark 3.18. This abuse of terminology should cause no confusion.

3.3 Some results from analysis

We need the following form of the implicit function theorem (IFT):

Lemma 3.19** (Implicit Function Theorem).**

Let $W\subseteq\mathbb{R}^{n}\times\mathbb{R}^{m}$ be open and $F\colon W\to\mathbb{R}^{m}$ be $C^{r}$ ($r\geq 1$). Suppose $F(a,b)=0$ for some $(a,b)\in W$ and the Jacobian matrix $D_{2}F(a,b)$ (namely with respect to the second variables) is nonsingular. Then there exist $U\subseteq\mathbb{R}^{n}$, $V\subseteq\mathbb{R}^{m}$ both open, with $(a,b)\in U\times V\subseteq W$, and a $C^{r}$ function $\phi\colon U\to V$ satisfying $\phi(a)=b$, and such that

[TABLE]

*namely, the zero set of $F$ in $U\times V$ is precisely the graph of $\phi$. *

Proof.

See, for example, Chapter 5 of [28].

The reader may easily verify that the sets $U$ and $V$ in the statement of the IFT may, without loss of generality, be chosen to be open balls with centres $a$ and $b$ in $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$ respectively.

We need the following consequence of the IFT.

Lemma 3.20** (IFT extended to a compact set).**

Let $W\subseteq\mathbb{R}^{n}\times\mathbb{R}^{m}$ be open and $F\colon W\to\mathbb{R}^{m}$ be $C^{r}$ ($r\geq 1$). Let $X$ be a compact set in $\mathbb{R}^{n}$ such that $X\times\{0\}\subseteq W$. Suppose that, for all $x\in X$, $F(x,0)=0$ and the Jacobian matrix $D_{2}F(x,0)$ is nonsingular for each $x\in X$. Then there exist an open set $U\subseteq\mathbb{R}^{n}$ containing $X$, $t>0$ such that $U\times B_{t}\subseteq W$, and a $C^{r}$ function $\phi\colon U\to B_{t}$ whose graph is precisely the zero-set of $F$ in $U\times B_{t}$, namely,

[TABLE]

Proof.

We apply the IFT at $(x,0)$ for each $x\in X$. For each $x\in X$, there exist $s_{x}>0$ and $t_{x}>0$ such that $B_{s_{x}}(x)\times B_{t_{x}}\subseteq W$, and a $C^{r}$ function $\phi_{x}\colon B_{s_{x}}(x)\to B_{t_{x}}$ such that the zero set of $F$ in $B_{s_{x}}(x)\times B_{t_{x}}$ is precisely the graph of $\phi_{x}$, namely,

[TABLE]

We now choose a finite set $\{x_{i}\}\subseteq X$ such that $U^{\prime}:=\cup B_{s_{x_{i}}}(x_{i})$ forms an open cover of $X$. We define the function $\hat{\phi}\colon U^{\prime}\to\mathbb{R}^{m}$ via $\hat{\phi}(x)=\phi_{x_{i}}(x)$ where $x_{i}$ is chosen as any element such that $x\in B_{s_{x_{i}}}(x_{i})$. $\hat{\phi}$ is a well defined function since if $x\in B_{s_{x_{i}}}(x_{i})\cap B_{s_{x_{j}}}(x_{j})$, then $\phi_{x_{i}}(x)=\phi_{x_{j}}(x)$ (as $\phi_{x_{i}}(x)$ and $\phi_{x_{j}}(x)$ must certainly both lie in one of $B_{t_{x_{i}}}$ or $B_{t_{x_{j}}}$). It is also clear that $\hat{\phi}$ is $C^{r}$ since it coincides with the $C^{r}$ functions $\{\phi_{x_{i}}\}$. Let $t:=\min\{t_{x_{i}}\}$. Since $\hat{\phi}$ is continuous, $X$ is compact, and $\hat{\phi}(x)=0$ for $x\in X$, there exists an open neighbourhood $U\subseteq U^{\prime}$ of $X$ such that $x\in U$ implies $|\hat{\phi}(x)|<t$. Define $\phi:=\left.\hat{\phi}\right|_{U}$. Clearly $\phi$ satisfies the claims of the lemma, and in particular $\{(x,y)\in U\times B_{t}\,\colon\,F(x,y)=0\}=\{(x,\phi(x))\,|\,x\in U\}\,.$

We will need the following technical lemma in order to make uniform estimates on compact sets. Notation is fixed to be consistent with those proofs where Lemma 3.21 is used.

Lemma 3.21**.**

Let $\mathcal{Z}\subseteq\mathbb{R}^{r}$ be compact, $m$ a positive integer, and $\eta^{\prime}>0$ a positive constant. Let $\theta\colon\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]\to\mathbb{R}^{m}$ satisfy

• •

$\theta$* is $C^{1}$ with Lipschitz continuous derivative on its domain of definition (see Remark 3.11). For example, if $\theta$ is $C^{2}$, then this condition certainly holds.*

• •

$\theta(z,0)=0$* for all $z\in\mathcal{Z}$.*

Then $\theta/\eta\colon\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]\backslash\{0\}\to\mathbb{R}^{m}$ has a continuous extension to $\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]$. Explicitly, $\hat{\theta}\colon\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]\to\mathbb{R}^{m}$ defined by

[TABLE]

*is continuous on $\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]$. Consequently $R(z,\eta):=\hat{\theta}(z,\eta)-D_{\eta}\theta(z,0)$ is continuous on $\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]$. *

Proof.

It is trivial that $\hat{\theta}$ is continuous (in fact, $C^{1}$) at points in its domain where $\eta\neq 0$. So we need to show that it is continuous at an arbitrary point of $\mathcal{Z}\times\{0\}$. Define $R(z,\eta):=\hat{\theta}(z,\eta)-D_{\eta}\theta(z,0)$ on $\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]$. Then $R(z,0)=0$ and $\theta(z,\eta)=\eta[D_{\eta}\theta(z,0)+R(z,\eta)]$ is an identity; on the other hand, Taylor’s theorem tells us that, for any fixed $z$, $\lim_{\eta\to 0}R(z,\eta)=0$. We would like to show that given any $z_{0}\in\mathcal{Z}$, $|\hat{\theta}(z,\eta)-\hat{\theta}(z_{0},0)|\to 0$ as $(z,\eta)\to(z_{0},0)$. If $\eta=0$, then

[TABLE]

which, by continuity of $D_{\eta}\theta$ at $(z_{0},0)$ can be made arbitrarily small by choosing $|z-z_{0}|$ sufficiently small. So now consider the case $\eta\neq 0$. We assume that $\eta>0$; the case $\eta<0$ requires minor modifications below. By the triangle inequality:

[TABLE]

The final term is simply the magnitude of the remainder $R(z_{0},\eta)$ in the Taylor expansion:

[TABLE]

On the other hand, using the fundamental theorem of calculus, we get

[TABLE]

Lipschitz continuity of $D_{\eta}\theta$ means that $\left|(D_{\eta}\theta(z,\sigma)-D_{\eta}\theta(z_{0},\sigma))\right|\leq K|z-z_{0}|$, where $K$ is the Lipschitz constant of $D_{\eta}\theta$ on $\mathcal{Z}\times[-\eta^{\prime},\eta^{\prime}]$. Thus

[TABLE]

We thus have from (5), (6) and (7):

[TABLE]

The first term on the RHS can be made small by choosing $|z-z_{0}|$ sufficiently small; the second term can be made small by choosing $\eta$ sufficiently small. This completes the proof that $\hat{\theta}$ is continuous on $\mathcal{Z}\times\{0\}$ and hence on its entire domain. As $\hat{\theta}$ is continuous and $D_{\eta}\theta(z,0)$ is continuous (as $\theta$ is $C^{1}$), $R$ is continuous as the difference of two continuous functions.

4 Proof of Theorem 1

The proof is presented with some surrounding explanation and, for readability, several subclaims are separated from the main proof into “subproofs”. This is to allow the reader to follow the main argument without necessarily digressing into the details of each technical claim. We break the proof into numbered parts which can be referred back to.

The basic set-up. We suppose that the original CRN $\mathcal{R}$ described by (1), namely, $\dot{x}=\Gamma v(x)$, admits a nondegenerate (resp., linearly stable) positive periodic orbit $\mathcal{O}$. Recall that $\Gamma$ is an $n\times r_{0}$ matrix, and that $v\colon\mathbb{R}^{n}_{\gg 0}\to\mathbb{R}^{r_{0}}$ is assumed to be $C^{2}$. Define $S_{\mathcal{O}}$ to be the coset of $\mathrm{im}\,\Gamma$ which contains $\mathcal{O}$, namely $S_{\mathcal{O}}=x_{0}+\mathrm{im}\,\Gamma$ for some $x_{0}\in\mathcal{O}$. Let $\mathcal{Z}^{o}$ be some connected subset of $S_{\mathcal{O}}$ containing $\mathcal{O}$, relatively open w.r.t. $S_{\mathcal{O}}$, and whose closure, $\mathcal{Z}$, is compact and lies in $\mathbb{R}^{n}_{\gg 0}$. We do not need to introduce local coordinates on $S_{\mathcal{O}}$ explicitly, although this can be done (see Remark 4.1). The geometry is illustrated schematically in Figure 3.

Recall the definitions of the matrices $a$, $a^{\prime}$, $b$, $b^{\prime}$, $\alpha$ and $\beta$. Given the assumption that the $(m+k)\times m$ matrix $\beta$ has rank $m$, we can assume without loss of generality (i.e., by reordering the added species $Y$ if necessary) that \beta=\left(\begin{array}[]{c}\hat{\beta}\\ \bar{\beta}\end{array}\right), where $\hat{\beta}$ is a nonsingular $m\times m$ matrix, and $\bar{\beta}$ is a $k\times m$ matrix. If $k=0$, then $\bar{\beta}$ is empty. $\hat{y}\in\mathbb{R}^{m}$, $\bar{y}\in\mathbb{R}^{k}$, $\hat{b}\in\mathbb{R}^{m\times m}$, $\hat{b}^{\prime}\in\mathbb{R}^{m\times m}$, $\bar{b}\in\mathbb{R}^{k\times m}$ and $\bar{b}^{\prime}\in\mathbb{R}^{k\times m}$ are defined in the natural way ($\bar{y}$, $\bar{b}$, and $\bar{b}^{\prime}$ are empty if $k=0$).

The stoichiometric matrix of $\mathcal{R}^{\prime}$ is

[TABLE]

Let $S^{*}$ denote the coset of $\mathrm{im}\,\Gamma^{\prime}$ which includes the point $(x,\hat{y},\bar{y})=(x_{0},0,\mathbf{1})$. (In the case $k=0$, this just means that $(x,\hat{y})=(x_{0},0)$.) Let $S^{*}_{+}$ denote the positive part of $S^{*}$, namely $S^{*}_{+}:=S^{*}\cap(\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0})$. $S^{*}_{+}$ is the positive stoichiometry class of interest to us: our goal is to show that we can choose rates for the added reactions from the class of mass action kinetics such that $\mathcal{R}^{\prime}$ admits a periodic orbit on $S^{*}_{+}$ which is nondegenerate (resp., linearly stable) relative to $S^{*}$. 2. 2.

A coordinate transformation to simplify the geometry. With the ultimate aim of setting up a singular perturbation problem amenable to the techniques of geometric singular perturbation theory, we now carry out a number of steps. The first is to define a new variable $z=x-\alpha\hat{\beta}^{-1}\hat{y}$. More precisely, we define the linear bijection $\phi_{0}\colon\mathbb{R}^{n}\times\mathbb{R}^{m+k}\to\mathbb{R}^{n}\times\mathbb{R}^{m+k}$ by

[TABLE]

We refer to the domain of $\phi_{0}$ as $(x,y)$-space and its codomain as $(z,y)$-space. It is easily shown that $\phi_{0}$ takes cosets of

[TABLE]

(See Subproof 4.2.) This transformation has permitted us to write the initial stoichiometric matrix in block diagonal form at the cost of somewhat complicating the subsequent rate functions. We will see shortly that this also allows us to formulate a singular perturbation problem. In $(z,y)$-space, (3) becomes

[TABLE]

(8) defines a local flow on $\phi_{0}(\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0})=\{(z,y)\in\mathbb{R}^{n}\times\mathbb{R}^{m+k}\colon z+\alpha\hat{\beta}^{-1}\hat{y}\gg 0,\,\,y\gg 0\}$. 3. 3.

Removing the variables $\bar{y}$. We now carry out a further transformation to eliminate the variable $\bar{y}$ from (8). The goal is to remove complications arising if the new reactions introduce new conservation laws. Conserved quantities involving only the new species can immediately be fixed at arbitrary positive values, eliminating some added species and simplifying later arguments. A similar strategy was adopted in the proof of Theorem 5 of [5]. Define $\delta:=-(\bar{\beta}\hat{\beta}^{-1})^{\mathrm{t}}$. Then $\delta^{\mathrm{t}}\hat{y}+\bar{y}$ is constant along trajectories of (8). (See Subproof 4.3.) We fix the value of $\delta^{\mathrm{t}}\hat{y}+\bar{y}$ as $\mathbf{1}$ (this choice is arbitrary: any fixed vector in $\mathbb{R}^{k}_{\gg 0}$ would do) and define the hyperplane

[TABLE]

Note that according to our conventions, $\mathcal{H}:=\mathbb{R}^{n}\times\mathbb{R}^{m}$ if $k=0$. Trajectories of (8) beginning on $\mathcal{H}$ remain on $\mathcal{H}$, as $\mathcal{H}$ is a union of cosets of $\mathrm{im}\,\Gamma\times\mathrm{im}\beta$ in $\mathbb{R}^{n}\times\mathbb{R}^{m+k}$. Define the affine bijection $\phi_{1}\colon\mathcal{H}\to\mathbb{R}^{n}\times\mathbb{R}^{m}$ by $\phi_{1}(z,(\hat{y},\mathbf{1}-\delta^{\mathrm{t}}\hat{y}))=(z,\hat{y})$. $\phi_{1}$ is the projection $(z,\hat{y},\bar{y})\mapsto(z,\hat{y})$ restricted to $\mathcal{H}$ and is just the identity on $\mathbb{R}^{n}\times\mathbb{R}^{m}$ if $k=0$. Define $S:=S_{\mathcal{O}}\times\mathbb{R}^{m}$, and observe that

[TABLE]

(See Subproof 4.4.) Also define $\mathcal{H}^{\prime}=\phi_{0}^{-1}(\mathcal{H})$, namely $\mathcal{H}^{\prime}$ is the subset of $(x,y)$-space defined by the equation $\delta^{\mathrm{t}}\hat{y}+\bar{y}=\mathbf{1}$. The action of the transformations $\phi_{0}$ and $\phi_{1}$ is summarised in Figure 4.

The claim that $\phi_{0}^{-1}\circ\phi_{1}^{-1}(S)=S^{*}$ is consistent with the following notational convention which we now adopt: given any set $X$ in $(z,\hat{y})$-space, $X^{*}$ refers to the corresponding set in $(x,y)$-space, namely $X^{*}=\phi_{0}^{-1}\circ\phi_{1}^{-1}(X)$. Notice that if $X\subseteq S$, then $X^{*}\subseteq S^{*}$. In the other direction, given a function $F$ on $(x,y)$-space, $F_{*}$ refers to the corresponding function on $(z,\hat{y})$-space, namely $F_{*}=F\circ\phi_{0}^{-1}\circ\phi_{1}^{-1}$. 4. 4.

The system in $(z,\hat{y})$-space.

Define the following subsets of $(z,\hat{y})$-space.

[TABLE]

or explicitly $W_{+}=\{(z,\hat{y})\in\mathbb{R}^{n}\times\mathbb{R}^{m}\colon z+\alpha\hat{\beta}^{-1}\hat{y}\gg 0,\,\,\hat{y}\gg 0,\,\,\mathbf{1}-\delta^{\mathrm{t}}\hat{y}\gg 0\}$. Similarly, define

[TABLE]

namely, $S_{+}:=\{(z,\hat{y})\in S\colon z+\alpha\hat{\beta}^{-1}\hat{y}\gg 0,\,\,\hat{y}\gg 0,\,\,\mathbf{1}-\delta^{\mathrm{t}}\hat{y}\gg 0\}$. Note that $S_{+}^{*}=\phi_{0}^{-1}\circ\phi_{1}^{-1}(S_{+})$ as expected from our notational conventions.

$W_{+}$ is the image in $(z,\hat{y})$-space of the positive part of $\mathcal{H}^{\prime}$, namely, $W_{+}$ is the image under $\phi_{1}\circ\phi_{0}$ of a union of positive stoichiometry classes; and $S_{+}$ is the image under $\phi_{1}\circ\phi_{0}$ in $(z,\hat{y})$-space of the positive part of $S^{*}$, the positive stoichiometry class of interest to us.

Define $q_{*}=q\circ\phi_{0}^{-1}\circ\phi_{1}^{-1}$, namely, $q_{*}(z,\hat{y})=q(z+\alpha\hat{\beta}^{-1}\hat{y},(\hat{y},\mathbf{1}-\delta^{\mathrm{t}}\hat{y}))$. As $q$ is defined on $\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0}$, $q_{*}$ is defined on $W_{+}$. The restriction of (8) to $\mathcal{H}$, followed by projection by $\phi_{1}$, gives the following system on $W_{+}$:

[TABLE]

We want to restrict attention to a region of $S_{+}$ on which the values of $\hat{y}$ are small and positive. To this end, we observe that there exists $y_{pos}>0$ such that

[TABLE]

lies in $S_{+}$, namely,

[TABLE]

(See Subproof 4.5.) Since $\mathcal{Z}_{+}\subseteq S_{+}$, it follows that $\mathcal{Z}_{+}^{*}\subseteq S^{*}_{+}$. Thus (10) defines a local flow on $\mathcal{Z}_{+}^{o}$, the relative interior of $\mathcal{Z}_{+}$ in $S$, and applying $\phi_{0}^{-1}\circ\phi_{1}^{-1}$ we get a corresponding local flow on a relatively open subset of $S^{*}_{+}$. If $\hat{y}\in\mathbb{R}$ for example, then $\mathcal{Z}_{+}=\mathcal{Z}\times(0,y_{pos}]$, and $\mathcal{Z}_{+}^{o}=\mathcal{Z}^{o}\times(0,y_{pos})$. The eventual goal of our constructions is to prove the existence of a periodic orbit of (10) in $\mathcal{Z}_{+}^{o}$ for appropriate choices of $q_{*}$. 5. 5.

Choosing the rates for the new reactions. The rate function $q$ of the new reactions has so far been left undetermined. $q$ will be chosen from some class of rate functions depending on two real parameters $\epsilon$ and $\eta$, and will take the form

[TABLE]

$f$ will be chosen shortly, but for the moment we assume that $f$ is defined and $C^{2}$ on $\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0}\times\mathbb{R}_{>0}$. The reason for introducing the two parameters is that, roughly speaking, we need to be able to make both rates of each reversible reaction arbitrarily large (via $\epsilon$) while independently controlling the ratio of forward and backward rates (via $\eta$). The intuition is that (i) we want to be able to ensure that the added reactions are close to equilibrium via $\epsilon$ (this will be our singular perturbation parameter); and (ii) we want to ensure, via $\eta$, that the values of $\hat{y}$ are small when the added reactions are at equilibrium.

In arguments below we sometimes extend $\phi_{0}$ and $\phi_{1}$ to $(x,y,\eta)$-space and $(z,y,\eta)$-space respectively, namely we let $\phi_{0}$ and $\phi_{1}$ refer to $\phi_{0}\times\mathrm{id}$ and $\phi_{1}\times\mathrm{id}$. This should cause no confusion.

Define $f_{*}:=f\circ\phi_{0}^{-1}\circ\phi_{1}^{-1}$, or more explicitly $f_{*}(z,\hat{y},\eta)=f(z+\alpha\hat{\beta}^{-1}\hat{y},(\hat{y},\mathbf{1}-\delta^{\mathrm{t}}\hat{y}),\eta)$. By construction $f_{*}$ is defined and $C^{2}$ on $W_{+}\times\mathbb{R}_{>0}$ which includes $\mathcal{Z}_{+}\times\mathbb{R}_{>0}$. (10) now takes the form of a typical singular perturbation problem:

[TABLE]

When we study (12) we restrict attention to $\mathcal{Z}_{+}^{o}$.

We have left fixing of the kinetics of the added reactions to this late stage in order to facilitate generalisation of the results. But now we assume mass action kinetics for the added reactions (2), and set

[TABLE] 6. 6.

The positive zero-set of $f$. We will be interested in the positive zero-set of $f$, corresponding to the added reactions being at equilibrium. Let $\gamma:=-(\alpha\,\hat{\beta}^{-1})^{\mathrm{t}}$ and recall that we defined $\delta=-(\bar{\beta}\,\hat{\beta}^{-1})^{\mathrm{t}}$. With some manipulation we find that for each fixed $\eta>0$, solutions to $f(x,y,\eta)=0$ on $\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0}$ are precisely solutions to

[TABLE]

on $\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m+k}_{\gg 0}$ (see Subproof 4.6.) Note that, for any $\gamma$ and $\delta$, $g$ is defined on a larger domain than $f$, namely, $\mathbb{R}^{n}_{\gg 0}\times(\mathbb{R}^{m}\times\mathbb{R}^{k}_{\gg 0})\times\mathbb{R}$ (this is to be interpreted as $\mathbb{R}^{n}_{\gg 0}\times\mathbb{R}^{m}\times\mathbb{R}$ if $k=0$). Define $g_{*}:=g\circ\phi_{0}^{-1}\circ\phi_{1}^{-1}$, namely,

[TABLE]

$g_{*}$ is defined and $C^{2}$ provided $z+\alpha\hat{\beta}^{-1}\hat{y}\gg 0$ and $\mathbf{1}-\delta^{\mathrm{t}}\hat{y}\gg 0$, and so the domain of $g_{*}$ includes an open neighbourhood of $\mathcal{Z}\times B_{y_{pos}}\times\mathbb{R}$ (recall Part 4). It is important in constructions below involving the IFT that the domain of $g_{*}$ includes $\mathcal{Z}\times\{0\}\times\{0\}$. 7. 7.

An upper bound on $\eta$. We need to put a number of upper bounds on $\eta$. The first of these ensures (via the IFT) that a portion of the zero set of $g_{*}$, and hence $f_{*}$, is the graph of a function. Let $\mathcal{E}$ be the zero set of $g_{*}$ on its domain, namely

[TABLE]

We now claim that there exists $y_{max}\in(0,y_{pos}]$ and $\eta_{1}>0$ such that the zero-set of $g_{*}$ in $\mathcal{Z}\times B_{y_{max}}\times(-\eta_{1},\eta_{1})$ is the graph of a $C^{2}$ function $\theta\colon\mathcal{Z}\times(-\eta_{1},\eta_{1})\to\mathbb{R}^{m}$, namely,

[TABLE]

This is shown by applying the extended IFT in Lemma 3.20 – (see Subproof 4.7). It is also clear that $\theta$ must satisfy $\theta(z,0)=0$. For each fixed $\eta\in(-\eta_{1},\eta_{1})$ we then have

[TABLE]

The geometry of the situation is illustrated schematically in Figure 5.

A second upper bound on $\eta$. This is needed to ensure that the function $\theta$ just constructed is strictly positive for positive $\eta$. In other words, we claim that there exists $\eta_{2}\in(0,\eta_{1}]$, such that $(z,\eta)\in\mathcal{Z}\times(0,\eta_{2})$ implies $\theta(z,\eta)\gg 0$ (see Subproof 4.8.) Recall, additionally, (Part 7 above) that $(z,\eta)\in\mathcal{Z}\times[0,\eta_{1})$ implies that $|\theta(z,\eta)|<y_{max}$. Thus $(z,\eta)\in\mathcal{Z}\times(0,\eta_{2})$ implies that $(z,\theta(z,\eta))\in\mathcal{Z}_{+}$. As a result, provided $\eta\in(0,\eta_{2})$, $\mathcal{E}_{\eta}\subseteq\mathcal{Z}_{+}\subseteq S_{+}$ (see Part 4), and consequently, defining $\mathcal{E}_{\eta}^{*}:=\phi_{0}^{-1}\circ\phi_{1}^{-1}(\mathcal{E}_{\eta})$, we have,

[TABLE] 9. 9.

A third upper bound on $\eta$. This is needed to ensure that the differential algebraic system obtained in a singular limit has a nondegenerate (resp., linearly stable) periodic orbit. Fix any $\eta\in(0,\eta_{2})$ and consider the following system

[TABLE]

on $\mathcal{Z}^{o}$. The reader may glance ahead to Part 11 of the proof and note that the vector field $F(z,\eta):=\Gamma v(z+\alpha\hat{\beta}^{-1}\theta(z,\eta))$ occurring on the right of (15) can be regarded as the so-called “reduced vector field” associated with (12). Note that $F$ is $C^{2}$ on $\mathcal{Z}\times(0,\eta_{2})$ and, by assumption, $\dot{z}=F(z,0)=\Gamma v(z)$ has a periodic orbit $\mathcal{O}$ in $\mathcal{Z}^{o}$, nondegenerate relative to $S_{\mathcal{O}}$. By regular perturbation theory arguments (Lemma 3.16) there exists $\eta_{3}\in(0,\eta_{2}]$ such that provided $\eta\in(0,\eta_{3})$ (15) has a periodic orbit $\mathcal{O}^{\prime}_{\eta}$ in $\mathcal{Z}^{o}$ close to $\mathcal{O}$, and such that the number of Floquet multipliers of $\mathcal{O}^{\prime}_{\eta}$ relative to $S_{\mathcal{O}}$ inside and outside the unit circle is the same as that of $\mathcal{O}$. Consequently if $\mathcal{O}$ is nondegenerate (resp., linearly stable) relative to $S_{\mathcal{O}}$, then so is $\mathcal{O}^{\prime}_{\eta}$. 10. 10.

A fourth upper bound on $\eta$. We need one more upper bound on $\eta$, connected with ensuring that normal hyperbolicity conditions needed to apply results in [27] hold. For each $(z,\eta)\in\mathcal{Z}\times(0,\eta_{3})$ define $W(z,\eta):=\hat{\beta}D_{\hat{y}}f_{*}(z,\theta(z,\eta),\eta)$, namely $W(z,\eta)$ is the Jacobian matrix of $\hat{\beta}f_{*}(z,\hat{y},\eta)$ w.r.t. $\hat{y}$, evaluated at $\mathcal{E}_{\eta}$. Since $f_{*}$ is defined and $C^{2}$ on $\mathcal{Z}_{+}\times\mathbb{R}_{>0}$, $D_{\hat{y}}f_{*}$ is defined and $C^{1}$ on $\mathcal{Z}_{+}\times\mathbb{R}_{>0}$. Since, additionally, $\theta$ is $C^{2}$ on $\mathcal{Z}\times(0,\eta_{3})$, $W$ is defined and $C^{1}$ on $\mathcal{Z}\times(0,\eta_{3})$. We claim that there exists $\eta_{4}\in(0,\eta_{3}]$, such that $(z,\eta)\in\mathcal{Z}\times(0,\eta_{4})$ implies that $W(z,\eta)$ is Hurwitz stable, namely the eigenvalues of $W(z,\eta)$ lie in the open left half of the complex plane. The calculations are fairly lengthy and are presented in Subproof 4.9. 11. 11.

Singular perturbation theory: completing the argument. We have done the preliminary work and are ready to apply perturbation theory results of Fenichel [27]. We fix some $\eta\in(0,\eta_{4})$, and return to system (12) on $\mathcal{Z}_{+}^{o}$, and the equivalent “fast time” system obtained by rescaling time. These are

[TABLE]

with their respective limiting systems in the limit $\epsilon\to 0+$:

[TABLE]

At this point, it may help to visit Remark 4.1: we could, if we chose, identify $S_{\mathcal{O}}$ with $\mathbb{R}^{r}$ and hence $S=S_{\mathcal{O}}\times\mathbb{R}^{m}$ with $\mathbb{R}^{r+m}$ via a linear change of coordinates, in which case $\mathcal{Z}^{o}$ would be identified with an open subset of $\mathbb{R}^{r}$ and $\mathcal{Z}_{+}^{o}$ would be identified with an open subset of $\mathbb{R}^{r+m}$. However, instead, we choose to keep our current coordinate system and bear in mind that we are restricting attention to the locally invariant set $\mathcal{Z}_{+}^{o}$. Claims about “nondegeneracy” and “linear stability” are relative to $S_{\mathcal{O}}$ or $S$ as appropriate.

(i) Note that $\mathcal{E}_{\eta}^{o}:=\{(z,\hat{y})\colon z\in\mathcal{Z}^{o},\,\,\hat{y}=\theta(z,\eta)\}$ is an invariant manifold of ($\mathrm{B_{0}}$) consisting entirely of equilibria. Our computation showing that $W(z,\eta)=\hat{\beta}D_{\hat{y}}f_{*}(z,\theta(z,\eta),\eta)$ is Hurwitz (Part 10 above) means that the eigenvalues associated with these equilibria corresponding to directions within $S$ but transverse to $\mathcal{E}_{\eta}^{o}$ all have negative real parts. The situation is illustrated schematically in Figure 6.

(ii) ($\mathrm{A_{0}}$) arises from the differential-algebraic system $\dot{z}=\Gamma v(z+\alpha\hat{\beta}^{-1}\hat{y}),\,\,0=\hat{\beta}f_{*}(z,\hat{y},\eta)$ bearing in mind that $\hat{\beta}$ is nonsingular and that solutions to $\hat{y}=\theta(z,\eta)$ satisfying $(z,\eta)\in\mathcal{Z}\times(0,\eta_{4})$ form a portion of the zero set of $f_{*}(z,\hat{y},\eta)$ (parts 6 and 7 above). For each fixed $\eta\in(0,\eta_{4})$, ($\mathrm{A_{0}}$) defines a local flow on $\mathcal{Z}^{o}$ which includes the periodic orbit $\mathcal{O}_{\eta}^{\prime}$ which is nondegenerate (resp., linearly stable) relative to $S_{\mathcal{O}}$ (recall Part 9 above). The projection of the vector field of ($\mathrm{A_{0}}$), namely $\Gamma v(z+\alpha\hat{\beta}^{-1}\theta(z,\eta))$, onto the tangent space of $\mathcal{E}_{\eta}^{o}$ is the reduced vector field associated with this system in the terminology of [27]. The reduced vector field has the periodic orbit

[TABLE]

on $\mathcal{E}_{\eta}^{o}$. $\mathcal{O}_{\eta}$ is nondegenerate (resp., linearly stable) relative to $\mathcal{E}_{\eta}^{o}$ because $\mathcal{O}^{\prime}_{\eta}$ is nondegenerate (resp., linearly stable) relative to $S_{\mathcal{O}}$ and $z\mapsto(z,\theta(z,\eta))$ is a $C^{1}$ diffeomorphism taking $\mathcal{Z}^{o}\subseteq S_{\mathcal{O}}$ to $\mathcal{E}_{\eta}^{o}$ (see Remark 3.18). The situation is illustrated in Figure 7.

According to Theorem 13.1 in [27], (i) and (ii) together tell us that the periodic orbit $\mathcal{O}_{\eta}$ “survives” perturbation namely, given any $\zeta>0$, we can choose $\epsilon_{0}>0$ such that for $\epsilon\in[0,\epsilon_{0})$, ($\mathrm{A}_{\epsilon}$) has a periodic orbit $\mathcal{O}_{\eta,\epsilon}$ satisfying $d_{H}(\mathcal{O}_{\eta},\mathcal{O}_{\eta,\epsilon})<\zeta$. In particular, we choose $\epsilon_{0}>0$ such that $\epsilon\in[0,\epsilon_{0})$ implies that $d_{H}(\mathcal{O}_{\eta},\mathcal{O}_{\eta,\epsilon})<d_{H}(\mathcal{O}_{\eta},\partial\mathcal{Z}_{+})$, and hence $\mathcal{O}_{\eta,\epsilon}\subseteq\mathcal{Z}_{+}^{o}$.

Moreover, according to Theorem 13.2 in [27], as $\mathcal{O}_{\eta}$ is nondegenerate (resp., linearly stable) relative to $S$, $\epsilon_{1}\in(0,\epsilon_{0}]$ can be chosen such that for all $\epsilon\in[0,\epsilon_{1})$, $\mathcal{O}_{\eta,\epsilon}$ is nondegenerate (resp., linearly stable) relative to $S$ as a periodic orbit of ($\mathrm{A}_{\epsilon}$). The claim about linear stability follows because the nontrivial eigenvalues of ($\mathrm{B_{0}}$) relative to $S$ at points of $\mathcal{E}_{\eta}$ all have negative real parts (Part 10 above).

We now fix $\epsilon\in(0,\epsilon_{1})$ and return to $(x,y)$-space. Let $\mathcal{O}_{\eta,\epsilon}^{*}:=\phi_{0}^{-1}\circ\phi_{1}^{-1}(\mathcal{O}_{\eta,\epsilon})$. Explicitly,

[TABLE]

Since $\mathcal{O}_{\eta,\epsilon}\subseteq\mathcal{Z}_{+}$, it follows that $\mathcal{O}_{\eta,\epsilon}^{*}\subseteq\mathcal{Z}_{+}^{*}\subseteq S_{+}^{*}$ (recall (14)), namely, $\mathcal{O}_{\eta,\epsilon}^{*}$ is a positive periodic orbit on $S^{*}$ for the enlarged CRN $\mathcal{R}^{\prime}$ governed by (3).

On the other hand, nondegeneracy (resp., linear stability) of $\mathcal{O}_{\eta,\epsilon}$ for ($\mathrm{A}_{\epsilon}$) relative to $S$, is equivalent (since $\phi_{0}^{-1}\circ\phi_{1}^{-1}$ is a $C^{1}$ diffeomorphism taking $S$ to $S^{*}$ – recall (9) and Lemma 3.17) to nondegeneracy (resp., linear stability) of $\mathcal{O}_{\eta,\epsilon}^{*}$ relative to $S^{*}$.

We have thus constructed a family of nondegenerate (resp., linearly stable) periodic orbits of $\mathcal{R}^{\prime}$ as desired. This completes the proof.

4.1 Remarks and subproofs

Remark 4.1**.**

In previous work [6], we introduced local coordinates on $S_{\mathcal{O}}$ as follows. Let $r$ be the rank of $\Gamma$, let $\Gamma_{0}$ be a matrix whose columns form a basis of $\mathrm{im}\,\Gamma$, and define the $r\times r_{0}$ matrix $Q$ with rank $r$ by $\Gamma=\Gamma_{0}Q$. Choose an arbitrary point $x_{0}\in S_{\mathcal{O}}$ and define $h\colon\mathbb{R}^{r}\to S_{\mathcal{O}}$ by $h(w)=x_{0}+\Gamma_{0}w$. Then $w\in\mathbb{R}^{r}$ defines a local coordinate on $S_{\mathcal{O}}$ which evolves according to the differential equation

[TABLE]

We could work with (16) rather than (1) as our starting point. However, we avoid explicitly introducing local coordinates as it is not strictly necessary here and introduces an extra transformation for us to track, obscuring the fundamental geometrical meaning of calculations. Instead, the reader may find it helpful to bear in mind that we can identify $S_{\mathcal{O}}$ with $\mathbb{R}^{r}$, $\mathcal{Z}$ with a compact subset of $\mathbb{R}^{r}$, and $\mathcal{Z}^{o}$ with an open subset of $\mathbb{R}^{r}$, all via some linear bijection such as $h$.

Subproof 4.2**.**

Given any $z_{0}\in\mathbb{R}^{n}$, $y_{0}\in\mathbb{R}^{m+k}$, $t\in\mathbb{R}^{r_{0}}$ and $s\in\mathbb{R}^{m}$,

[TABLE]

$\square$

Subproof 4.3**.**

The claim that $\delta^{\mathrm{t}}\hat{y}+\bar{y}$ is constant along trajectories of (8) is vacuously true in the case $k=0$. So, suppose that $k>0$. Let $P:=[-\bar{\beta}\hat{\beta}^{-1}|I_{k}]\in\mathbb{R}^{k\times(m+k)}$. By a quick calculation, $P\beta=0$, and since $P$ has rank $k$ and $\mathrm{dim}(\mathrm{ker}\,\beta^{\mathrm{t}})=k$ by the rank-nullity theorem, the rows of $P$ must form a basis of $\mathrm{ker}\,\beta^{\mathrm{t}}$. Multiplying $\dot{y}=\beta q(x,y)$ by $P$ on the left gives $P\dot{y}=0$, namely $Py=-\bar{\beta}\hat{\beta}^{-1}\hat{y}+\bar{y}$ is constant along trajectories. Thus, the value of $\bar{y}$ along a trajectory at any point is specified by the value of $\hat{y}$ on the trajectory at that point, and an additional parameter, $\delta^{\mathrm{t}}\hat{y}+\bar{y}\in\mathbb{R}^{k}$. $\square$

Subproof 4.4**.**

Let $x_{0}$ be some arbitrary point on $S_{\mathcal{O}}$. Note that

[TABLE]

Thus $(x,y)\in\phi_{0}^{-1}\circ\phi_{1}^{-1}(S)$, implies that $x-\alpha\hat{\beta}^{-1}\hat{y}=x_{0}+\Gamma t$ for some $t\in\mathbb{R}^{r_{0}}$, namely $x=x_{0}+\Gamma t+\alpha\hat{\beta}^{-1}\hat{y}$, and $\bar{y}=\mathbf{1}-\delta^{\mathrm{t}}\hat{y}$, so that

[TABLE]

Thus $(x,y)\in\phi_{0}^{-1}\circ\phi_{1}^{-1}(S)\Rightarrow(x,y)\in S^{*}$. On the other hand, $(x,y)\in S^{*}$ means, by definition, that

[TABLE]

for some $t\in\mathbb{R}^{r_{0}}$ and $s\in\mathbb{R}^{m}$, from which we see that $x-\alpha\hat{\beta}^{-1}\hat{y}=x_{0}+\Gamma t$ (namely, $x-\alpha\hat{\beta}^{-1}\hat{y}\in S_{\mathcal{O}}$), and $\delta^{\mathrm{t}}\hat{y}+\bar{y}=\mathbf{1}$. Thus $(x,y)\in S^{*}\Rightarrow(x,y)\in\phi_{0}^{-1}\circ\phi_{1}^{-1}(S)$, confirming the claim. $\square$

Subproof 4.5**.**

Let $D=d_{H}(\mathcal{Z},\partial\mathbb{R}^{n}_{\geq 0})$, and let $\|\cdot\|$ refer to the matrix norm induced by the Euclidean norm. If $k=0$, then define $y_{pos}:=D/(2\|\alpha\hat{\beta}^{-1}\|)\,.$ Otherwise, set

[TABLE]

Then $(z,\hat{y})\in\mathcal{Z}_{+}$ implies that (i) $z+\alpha\hat{\beta}^{-1}\hat{y}\gg 0$ and (ii) $\delta^{\mathrm{t}}\hat{y}\ll\mathbf{1}$ and, consequently, $\mathbf{1}-\delta^{\mathrm{t}}\hat{y}\gg 0\,.$ Consequently,

[TABLE]

$\square$

Subproof 4.6**.**

Solving $f(x,y,\eta)=0$ with the assumption that $x\gg 0$ and $y\gg 0$ gives $y^{\beta^{\mathrm{t}}}=\bm{\eta}^{\hat{\beta}^{\mathrm{t}}}\circ x^{-\alpha^{\mathrm{t}}}$. This can be written $\hat{y}^{\hat{\beta}^{\mathrm{t}}}=\bm{\eta}^{\hat{\beta}^{\mathrm{t}}}\circ x^{-\alpha^{\mathrm{t}}}\circ\bar{y}^{-\bar{\beta}^{\mathrm{t}}}$. Taking logs gives:

[TABLE]

Multiplying through by $(\hat{\beta}^{\mathrm{t}})^{-1}$ and exponentiating again gives the result. $\square$

Subproof 4.7**.**

First, we calculate that $g_{*}(z,0,0)=0$ for any $z$ in the domain of $g_{*}$, and that

[TABLE]

We now apply the extension of the IFT in Lemma 3.20. Note that given any neighbourhood $V$ of $\mathcal{Z}\times\{0\}$ in $\mathbb{R}^{n}\times\mathbb{R}$, there exists $\eta_{1}>0$ such that $\mathcal{Z}\times(-\eta_{1},\eta_{1})\subseteq V$. Consequently, since $g_{*}$ is $C^{2}$, according to Lemma 3.20 there exists $\eta_{1}>0$, $y_{max}\in(0,y_{pos}]$ and a $C^{2}$ function $\theta\colon\mathcal{Z}\times(-\eta_{1},\eta_{1})\to B_{y_{max}}$ whose graph is precisely the zero-set of $g_{*}$ in $\mathcal{Z}\times B_{y_{max}}\times(-\eta_{1},\eta_{1})$, namely,

[TABLE]

$\square$

Subproof 4.8**.**

We find $\eta_{2}\in(0,\eta_{1}]$, such that $\theta(z,\eta)\gg 0$ for $(z,\eta)\in\mathcal{Z}\times(0,\eta_{2})$. Consider the function $\hat{\theta}\colon\mathcal{Z}\times(-\eta_{1},\eta_{1})$ defined by

[TABLE]

Since $\theta$ is $C^{2}$ on its domain and $\theta(z,0)=0$ for all $z\in\mathcal{Z}$ we have, by Lemma 3.21, that $\hat{\theta}$ is continuous on its domain. On the other hand $\theta(\mathcal{Z},0)\subseteq\mathbb{R}^{m}_{\gg 0}$: implicitly differentiating $g_{*}(z,\theta(z,\eta),\eta)=0$ w.r.t. $\eta$, evaluating at $(z,0)$ for any $z\in\mathcal{Z}$, recalling that $D_{\hat{y}}g_{*}(z,0,0)=\mathrm{id}$, and using the fact that $\theta(z,0)=0$ gives

[TABLE]

If $\hat{\theta}(z,\eta)\gg 0$ on $\mathcal{Z}\times(-\eta_{1},\eta_{1})$, then $\theta(z,\eta)\gg 0$ on $\mathcal{Z}\times(0,\eta_{1})$ and we define $\eta_{2}=\eta_{1}$. Otherwise, we define

[TABLE]

The last inequality follows by continuity of $\hat{\theta}$: the pre-image of the closed set $\partial\mathbb{R}^{m}_{\geq 0}$ under $\hat{\theta}$ is closed, and would otherwise have to include a point in $\mathcal{Z}\times\{0\}$, contradicting $\hat{\theta}(\mathcal{Z},0)\subseteq\mathbb{R}^{m}_{\gg 0}$. $\square$

Subproof 4.9**.**

We can rewrite $f$ as a function of four variables $x$, $\hat{y},\bar{y}$ and $\eta$, namely

[TABLE]

Noting that $f_{*}(z,\hat{y},\eta)=\overline{f}(z+\alpha\hat{\beta}^{-1}\hat{y},\hat{y},\mathbf{1}-\delta^{\mathrm{t}}\hat{y},\eta)$, we get

[TABLE]

$W\colon\mathcal{Z}\times(0,\eta_{3})\to\mathbb{R}^{m\times m}$ is defined and continuous. We want to show that $\eta W$ has a continuous extension to $\mathcal{Z}\times[0,\eta_{3})$, namely, there exists a continuous function $\overline{W}\colon\mathcal{Z}\times[0,\eta_{3})\to\mathbb{R}^{m\times m}$ which coincides with $\eta W$ on $\mathcal{Z}\times(0,\eta_{3})$. Further, $\overline{W}(z,0)$ has real, negative eigenvalues for each fixed $z\in\mathcal{Z}$. By continuity of $\overline{W}$ and compactness of $\mathcal{Z}$, it will follow that there exists $\eta_{4}\in(0,\eta_{3}]$ such that for $(z,\eta)\in\mathcal{Z}\times(0,\eta_{4})$, eigenvalues of $\overline{W}(z,\eta)$ lie in the open left half of the complex plane. The same must hold for $W(z,\eta)$ whose eigenvalues are positive multiples of those of $\overline{W}(z,\eta)$.

Note first that $f_{*}(z,\theta(z,\eta),\eta)=0$ for $(z,\eta)\in\mathcal{Z}\times(0,\eta_{3})$, implying that

[TABLE]

$T^{\prime}\colon\mathcal{Z}\times(0,\eta_{3})\to\mathbb{R}^{m}$ has a continuous extension to $\mathcal{Z}\times[0,\eta_{3})$. To see this note, via Lemma 3.21 and Subproof 4.8, that

[TABLE]

where $R(z,\eta)$ is continuous on $\mathcal{Z}\times[0,\eta_{3})$ and satisfies $R(z,0)=0$. Moreover $\hat{\theta}(z,\eta):=z^{\gamma}+R(z,\eta)$ is continuous and positive on $\mathcal{Z}\times[0,\eta_{3})$ (see Subproof 4.8). So

[TABLE]

is defined and continuous on $\mathcal{Z}\times[0,\eta_{3})$ and clearly extends $T^{\prime}(z,\eta)$. Simple evaluation gives

[TABLE]

We are now ready to compute $\overline{W}(z,\eta)$. Differentiating $\overline{f}$ w.r.t. its first argument (see Example 3.15) gives

[TABLE]

where $T_{1}$ and $T_{2}$ are abbreviations for the first and second terms in $\overline{f}$ respectively. But when $\hat{y}=\theta(z,\eta)$ and $\bar{y}=\mathbf{1}-\delta^{\mathrm{t}}\theta(z,\eta)$, then $\overline{f}=0$, and $T_{1}=T_{2}=T(z,\eta)$. Define $D_{z,\eta}:=\mathrm{diag}(T(z,\eta))$ on $\mathcal{Z}\times[0,\eta_{3})$. We calculate:

[TABLE]

This quantity is defined and continuous on $\mathcal{Z}\times[0,\eta_{3})$ since $(z,\eta)\in\mathcal{Z}\times[0,\eta_{3})$ implies that $(z,\theta(z,\eta))\in\mathcal{Z}_{+}$ and consequently $z+\alpha\hat{\beta}^{-1}\theta(z,\eta)\gg 0$ (see Subproof 4.5). Similarly,

[TABLE]

This is again defined and continuous on $\mathcal{Z}\times[0,\eta_{3})$ (we have already observed in Subproof 4.8 that $\hat{\theta}$ is continuous and positive on $\mathcal{Z}\times[0,\eta_{3})$.) Finally,

[TABLE]

(Note that if $k=0$, then according to our conventions $D_{3}\overline{f}$ is a zero matrix.) Again, this quantity is defined and continuous on $\mathcal{Z}\times[0,\eta_{3})$ since $(z,\eta)\in\mathcal{Z}\times[0,\eta_{3})$ implies that $(z,\theta(z,\eta))\in\mathcal{Z}_{+}$ and consequently $\mathbf{1}-\delta^{\mathrm{t}}\theta(z,\eta)$ is positive (see Subproof 4.5). We thus have that

[TABLE]

is defined and continous on $\mathcal{Z}\times[0,\eta_{3})$ and coincides with $\eta W(z,\eta)$ on $(0,\eta_{3})$. Moreover, as calculated in Subproof 4.8, $\hat{\theta}(z,0)=z^{\gamma}\gg 0$, and so

[TABLE]

Defining $Y_{z}:=\mathrm{diag}(\mathbf{1}/z^{\gamma})$, we have

[TABLE]

where $M_{z}:=Y_{z}^{1/2}\hat{\beta}(D_{z,0})^{1/2}$ is nonsingular as each of its factors is nonsingular. From this we see that $\overline{W}(z,0)$ is similar to a negative definite matrix, and so its eigenvalues are real and negative as claimed. $\square$

5 Elucidating the proof via an example

We now return to the example in Section 2 and use it to elucidate the proof of Theorem 1. We derive for this example many of the matrices, sets and functions which occur in the proof above. Recall the CRNs

[TABLE]

and

[TABLE]

We fix the order $[X,Y,Z,U,V,W]$ on the chemical species and let the corresponding lower case letters refer to chemical concentrations. The stoichiometric matrices (Part 1 above) of $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are, respectively

[TABLE]

We also have following matrices defined in Part 1 of the proof:

[TABLE]

and

[TABLE]

From these we obtain (Parts 3 and 6 above):

[TABLE]

The added reactions give rise to the conservation relation (Part 3)

[TABLE]

$v$, the rate vector of $\mathcal{R}_{1}$, is given by

[TABLE]

Thus, the ODE system for $\mathcal{R}_{1}$ is

[TABLE]

As $\Gamma$ has rank $3$, the stoichiometric subspace is the whole of $(x,y,z)$-space, and consequently, the unique positive stoichiometry class of the system is the positive orthant in $(x,y,z)$-space. We know that for some choice of rate constants this system has a linearly stable periodic orbit.

Turning to $\mathcal{R}_{2}$, the rate vector $q$ of the new reactions, set according to Part 5, is

[TABLE]

Thus the ODE system for $\mathcal{R}_{2}$ is

[TABLE]

Theorem 1 tells us that for sufficiently small $\epsilon$ and $\eta$ (21) has a positive periodic orbit on the set defined by $u/3-v/3+w=1$, and that this periodic orbit is linearly stable relative to its stoichiometry class.

Note at this stage that (21) does not have the form of a singular perturbation problem if we treat $\epsilon$ as a small parameter. In order to remedy this we follow Part 2 of the proof and introduce new variables $(x^{\prime},y^{\prime},z^{\prime})$ defined by

[TABLE]

Conversely

[TABLE]

The dynamical system in the new variables becomes:

[TABLE]

Note the block diagonal structure of the new “stoichiometric matrix”. The system now has the form of a singular perturbation problem with $\epsilon$ as the perturbation parameter. $w$ can be eliminated by fixing the conserved quantity and setting $w=1-u/3-v/3$ to give:

[TABLE]

This singularly perturbed set of ODEs is, for this example, the system referred to as ($\mathrm{A}_{\epsilon}$) in Part 11 of the proof. The function $f_{*}$ referred to in Part 5 is, in this case, the following function of $x^{\prime},y^{\prime},z^{\prime},u,v$ and $\eta$ which occurs in the second part of (22).

[TABLE]

In the limit $\epsilon\to 0$, (22) reduces to the differential-algebraic system

[TABLE]

while the associated fast-time system gives the ODEs:

[TABLE]

(24) is the system referred to as ($\mathrm{B}_{0}$) in Part 11 of the proof.

When we observe (23), our instinct is to “solve” the algebraic equations for $u$ and $v$ in terms of $x^{\prime},y^{\prime}$ and $z^{\prime}$, and substitute into the differential equations. Several steps of the proof are focussed on showing that for fixed small $\eta$ this is indeed possible. Note, however, that even for this simple example explicit solution is difficult, and the sensible approach is to appeal to the implicit function theorem to solve the equations locally.

The function $g$ defined in Part 6 whose positive zeros are precisely the positive equilibria of the new reactions takes the form

[TABLE]

Substituting in $w=1-u/3+v/3$ from the new conservation relation and using the new variables $x^{\prime},y^{\prime},z^{\prime}$ we obtain for this example the function referred to as $g_{*}$ in Part 6 of the proof:

[TABLE]

Note that $g_{*}$ is defined provided $x^{\prime}-v/3+u/3>0$, $y^{\prime}-v/3-2u/3>0$ and $1-u/3+v/3>0$. We are interested in the zeros of $g_{*}$ satisfying, additionally, $u>0$, $v>0$ and $z^{\prime}>0$, corresponding to the original variables $u,v,w,x,y,z$ all being positive.

$g_{*}$ (unlike $f_{*}$) is defined when $\eta=0$ and $g_{*}(x^{\prime},y^{\prime},z^{\prime},0,0,0)=0$. In Part 7 we find an upper bound $\eta_{1}$ on $\eta$ ensuring, via the IFT, that for $x^{\prime},y^{\prime}$ and $z^{\prime}$ lying in the compact positive set termed $\mathcal{Z}$ in the proof, and $0<\eta<\eta_{1}$, we can solve $g_{*}=0$ for $u$ and $v$ in terms of $x^{\prime},y^{\prime},z^{\prime}$ and $\eta$ (although $z^{\prime}$ does not explicitly occur here). In other words, we can locally describe the zero set of $g_{*}$ as the graph of a function, termed $\theta$ in the proof, with input $(x^{\prime},y^{\prime},z^{\prime},\eta)$ and output $(u,v)$. Moreover, in the region where this description holds $x=x^{\prime}-v/3+u/3>0$, $y=y^{\prime}-v/3-2u/3>0$, $z=z^{\prime}>0$ and $w=1-u/3+v/3>0$. The fact that we solve $g_{*}=0$ not just in the neighbourhood of a single point, but in a neighbourhood of a larger set necessitates use of the slight extension to the IFT in Lemma 3.20.

The construction in Part 8 of the proof gives another upper bound $\eta_{2}$ on $\eta$ which ensures that $\theta$ is a positive function. In other words, for $0<\eta<\eta_{2}$ the solutions to $g_{*}=0$ obtained in the previous part additionally satisfy $u>0$ and $v>0$. Thus for $x^{\prime},y^{\prime}$ and $z^{\prime}$ lying in the compact positive set $\mathcal{Z}$, and $\eta$ sufficiently small and positive, the graph of $\theta$ corresponds to a positive subset of the original $(x,y,z,u,v,w)$-space.

To understand the third upper bound on $\eta$ in Part 9, consider the “reduced” vector field in a portion of $(x^{\prime},y^{\prime},z^{\prime})$-space defined by solving the two algebraic equations of (23) for $u$ and $v$ in terms of $x^{\prime},y^{\prime},z^{\prime}$ (possible by the previous arguments for each small enough, fixed, positive $\eta$) and substituting into the differential equations of (23). Notice that since $g_{*}(x^{\prime},y^{\prime},z^{\prime},0,0,0)=0$, the values of $u$ and $v$ on the graph of $\theta$ are small for small $\eta$, and consequently the reduced vector field is close to the original vector field in (20). The third upper bound on $\eta$ ensures that the two are close enough to ensure, by regular perturbation theory arguments, that the periodic orbit of (20) survives in the reduced system and retains the same number of multipliers inside and outside the unit circle. Thus the reduced system has a periodic orbit of the same stability type as the original vector field (20) for sufficiently small values of $\eta$.

The fourth and final upper bound on $\eta$ in Part 10 is, roughly speaking, a normal hyperbolicity condition. To understand it, fix any sufficiently small $\eta>0$ and consider (24), namely the fast-time system in the limit $\epsilon\to 0$. Note that setting $\dot{u}=\dot{v}=0$ in (24) defines the set of equilibria of (24), which includes a portion of the graph of the function $\theta$ constructed above. Each equilibrium on this graph has a set of eigenvalues [math] corresponding to directions tangential to the graph and a set of “nontrivial” eigenvalues corresponding to directions transverse to the graph. This final upper bound on $\eta$ ensures that all of the nontrivial eigenvalues have negative real part, and so this portion of the equilibrium manifold of (24) is locally, exponentially attracting. This allows use of the perturbation theory results of [27]: it ensures, roughly speaking, that the manifold defined by the graph of $\theta$ survives for sufficiently small $\epsilon>0$ in (22) and has sufficiently smooth dependence on $\epsilon$. Consequently, by regular perturbation theory, the periodic orbit of (23) on this graph survives for (22) with sufficiently small $\epsilon>0$.

Since the original periodic orbit of (20) was linearly stable, so is the perturbed one for (22) – the previous two bounds on $\eta$ took care of its Floquet multipliers tangential to, and transverse to, the invariant set on which it lies respectively. Thus for sufficiently small, positive $\eta$ and $\epsilon$, (22) has a linearly stable periodic orbit. Moreover, for small enough $\epsilon>0$, this periodic orbit lies in a region of $(x^{\prime},y^{\prime},z^{\prime},u,v)$-space corresponding to positive $(x,y,z,u,v,w)$-space, since it can be made arbitrarily close to the periodic orbit of the limiting system (23). We have thus obtained, as desired, a periodic orbit of (21) which is positive and linearly stable relative to its stoichiometry class.

6 Final remarks

We remark, first, that giving the added reactions mass action kinetics was convenient and simplified many calculations, but was not fundamental to the techniques of proof of Theorem 1. The motivated reader could, with some effort, reprove the result with the added reactions having kinetics from other classes than mass action. Key to the proof is scalability of the reaction rates, and the characterisation of the equilibrium set of the added reactions as a graph over $\mathcal{Z}$.

The techniques of the proof of Theorem 1 also provide an alternative proof of Theorem 5 in [5]. The set-up requires only minor and formal modifications: $\mathcal{Z}$ is now a compact subset of some stoichiometry class of $\mathcal{R}$ containing two nondegenerate (resp., linearly stable) equilibria (rather than a periodic orbit), and we need to apply Theorems 12.1 and 12.2 of [27] (rather than Theorems 13.1 and 13.2). Using the same approach and the very general and powerful Theorem 9.1 in [27], essentially any compact limit set of (1), hyperbolic relative to its stoichiometry class, survives for (3) under the assumptions of Theorem 1.

The condition that $\beta$ has rank $m$ can be rephrased, roughly, as “the new species feature nondegenerately in the new reactions”. Although this condition is essential to our proof (at several points we invert $\hat{\beta}$), it is possible that a weaker condition might suffice here. The difficulty in exploring this question arises partly because techniques for proving the nonexistence of periodic orbits are limited.

In Theorem 6 of [5] we proved that nondegenerate (resp., linearly stable) multistationarity is preserved by “splitting” reactions and adding in intermediate complexes involving new species, provided this is done in a way satisfying a nondegeneracy condition very similar to the condition on the rank of $\beta$ in Theorem 1. We believe an analogous result for periodic orbits to hold, and certain special cases are easily proved. However, the result in full generality cannot be proved using the arguments in the proof of Theorem 1 here. The set up in Theorem 6 of [5] does not afford sufficient freedom to control reaction rates to apply the techniques of proof used in Theorem 1 here: the construction involving two independently controlled parameters, $\eta$ and $\epsilon$ cannot be simply reused, and an alternative approach needs to be found. This task will be undertaken in future work.

Acknowledgements

I would like to thank Casian Pantea, Amlan Banaji, and the anonymous reviewers of this paper for helpful discussions and useful comments on drafts of this paper.