††thanks: BB was supported by the Austrian Science Fund (FWF), project P32532.

# Some minimal bimolecular mass-action
systems with limit cycles

Balázs Boros  Department of Mathematics, University of Vienna
Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria<EMAIL_ADDRESS>Josef Hofbauer  Department of Mathematics, University of Vienna
Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria<EMAIL_ADDRESS>
###### Abstract.

We discuss three examples of bimolecular mass-action systems with three
species, due to Feinberg, Berner, Heinrich, and Wilhelm. Each system has a
unique positive equilibrium which is unstable for certain rate constants and
then exhibits stable limit cycles, but no chaotic behaviour. For some rate
constants in the Feinberg–Berner system, a stable equilibrium, an unstabe
limit cycle, and a stable limit cycle coexist. All three networks are minimal
in some sense.

By way of homogenising the three known examples, we construct bimolecular
mass-conserving mass-action systems with four species that admit a stable
limit cycle. The homogenised Feinberg–Berner system and the homogenised
Wilhelm–Heinrich system admit the coexistence of a stable equilibrium, an
unstable limit cycle, and a stable limit cycle.

###### Key words and phrases:

bimolecular network, competitive system, permanence, limit cycle,
Andronov–Hopf bifurcation, focal value

###### 1991 Mathematics Subject Classification:

34C25, 34C12, 34C23

##  1\. Introduction

We proved recently that rank-two bimolecular mass-action systems do not admit
limit cycles [6, Section 4]. Thus, to construct simple mass-action systems
with limit cycles, one could focus either on rank-two trimolecular (or
tetramolecular) networks or rank-three bimolecular networks. The former
direction is followed in [6, Section 3], while the latter one is the topic of
the present paper. We analyse three known three-species bimolecular mass-
action systems that are all minimal in some sense:

  1. (i)

The underlying reaction network of the _Feinberg–Berner oscillator_ (studied
in Section 3) is reversible with $5$ complexes and $4$ reversible reactions.
As we will discuss, any reversible bimolecular mass-action system with a limit
cycle has at least $5$ complexes and at least $4$ reversible reactions.

  2. (ii)

The r.h.s. of the differential equation of the _Wilhelm–Heinrich oscillator_
(studied in Section 4) has only one nonlinear term. No linear differential
equation admits a limit cycle.

  3. (iii)

The underlying reaction network of the _Wilhelm oscillator_ (studied in
Section 5) has only $4$ reactions. As we will discuss, any bimolecular mass-
action system with a limit cycle has at least $4$ reactions.

Having available the inheritance results summarised in [2, Theorem 2], it is
apparent that finding small networks with limit cycles is a useful tool for
studying the capacity for oscillations in more complex networks.

Since physically realistic systems are often mass-conserving, it is also
interesting to find small mass-conserving bimolecular networks whose
associated mass-action system admits limit cycle oscillation. Following [3,
Section 4], we homogenise each of the three three-species oscillators
mentioned above, and thereby construct three four-species mass-conserving
bimolecular mass-action systems that all exhibit limit cycles.

After collecting the necessary preliminaries in Section 2, we study the
(homogenised) Feinberg–Berner oscillator, the (homogenised) Wilhelm–Heinrich
oscillator, and the (homogenised) Wilhelm oscillator in Sections 3, 4 and 5,
respectively. Appendix A is devoted to recall the computation of the first
focal value at an equilibrium of a $3$-dimensional differential equation.

##  2\. Preliminaries

In this section we briefly summarise the necessary background and terminology.
For more details on reaction network theory, consult e.g. [29] and the
references therein. The symbols $\mathbb{R}_{+}$, $\mathbb{R}_{\geq 0}$, and
$\mathbb{Z}_{\geq 0}$ denote the set of positive reals, nonnegative reals, and
nonnegative integers, respectively.

###  2.1. Reaction networks and mass-action systems

A _reaction network_ (or _network_ for short) consists of species, complexes,
and reactions. For given _species_ $\mathsf{X}_{1},\ldots,\mathsf{X}_{n}$,
_complexes_ are formal linear combinations
$a_{1}\mathsf{X}_{1}+\cdots+a_{n}\mathsf{X}_{n}$ with nonnegative integers
$a_{1},\ldots,a_{n}$. A _reaction_ is an ordered pair of a reactant complex
and a product complex. With the symbols $V\subseteq\mathbb{Z}^{n}_{\geq 0}$
and $E\subseteq V\times V$ denoting the set of complexes and the set of
reactions, respectively, $(V,E)$ is a directed graph. The _reaction vector_
associated to a reaction $a_{1}\mathsf{X}_{1}+\cdots+a_{n}\mathsf{X}_{n}\to
b_{1}\mathsf{X}_{1}+\cdots+b_{n}\mathsf{X}_{n}$ is the vector
$(b_{1}-a_{1},\ldots,b_{n}-a_{n})^{\top}$. The _stoichiometric matrix_ is
denoted by $\Gamma$, its columns are the reaction vectors. By definition, the
_rank_ of a network is $\operatorname{rank}\Gamma$.

The vector $x(\tau)\in\mathbb{R}^{n}_{\geq 0}$ encodes the concentrations of
the $n$ species at time $\tau$. Assuming _mass-action kinetics_ , the _rate
function_ associated to a reaction
$a_{1}\mathsf{X}_{1}+\cdots+a_{n}\mathsf{X}_{n}\to
b_{1}\mathsf{X}_{1}+\cdots+b_{n}\mathsf{X}_{n}$ is $\kappa x_{1}^{a_{1}}\cdots
x_{n}^{a_{n}}$, where $\kappa>0$ is called the _rate constant_ of the
reaction. The rate functions of the individual reactions are collected in the
vector $v(x)$. The species concentrations then evolve according to the mass-
action differential equation $\dot{x}(\tau)=\Gamma v(x(\tau))$ with state
space $\mathbb{R}^{n}_{\geq 0}$. Notice that its solutions are confined to the
linear manifolds $(p+\operatorname{im}\Gamma)\cap\mathbb{R}^{n}_{\geq 0}$ with
$p\in\mathbb{R}^{n}_{\geq 0}$, termed _stoichiometric classes_. It is known
that $(p+\operatorname{im}\Gamma)\cap\mathbb{R}^{n}_{+}$ with
$p\in\mathbb{R}^{n}_{+}$ is forward invariant. A reaction network is _mass-
conserving_ if the stoichiometric classes are bounded (or equivalently, there
exists a $d\in\mathbb{R}^{n}_{+}$ such that $d^{\top}\Gamma=0$).

###  2.2. Deficiency

We now recall some special cases of the classical Deficiency-Zero Theorem [9,
Theorem 3.1] and the classical Deficiency-One Theorem [9, Theorem 3.2]. The
_deficiency_ of a network $(V,E)$ is the nonnegative integer
$|V|-\ell-\operatorname{rank}\Gamma$, where $\ell$ is the number of connected
components of the directed graph $(V,E)$. A directed graph is _strongly
connected_ if for any two vertices $i$ and $j$ there exists a directed path
from $i$ to $j$.

######  Theorem 1.

Deficiency-zero mass-action systems do not admit periodic solutions in the
positive orthant $\mathbb{R}^{n}_{+}$.

######  Theorem 2.

If the deficiency-one network $(V,E)$ is strongly connected and
$\operatorname{rank}\Gamma=n$ holds then its associated mass-action system has
exactly one positive equilibrium.

###  2.3. Permanence

The mass-action differential equation $\dot{x}(\tau)=\Gamma v(x(\tau))$ is
said to be _permanent_ if for every stoichiometric class
$\mathcal{P}=(p+\operatorname{im}\Gamma)\cap\mathbb{R}^{n}_{\geq 0}$ with
$p\in\mathbb{R}^{n}_{+}$ there exists a compact and forward invariant set
$K\subseteq\mathcal{P}\cap\mathbb{R}^{n}_{+}$ with the property that for each
solution $\tau\mapsto x(\tau)$ with
$x(0)\in\mathcal{P}\cap\mathbb{R}^{n}_{+}$, there exists a $\tau_{0}\geq 0$
such that for all $\tau\geq\tau_{0}$ we have $x(\tau)\in K$. Due to the
following theorem, a mass-action system is permanent provided that its
underlying reaction network is strongly connected (see [17, Theorem 1.3], [1,
Theorem 5.5], or [5, Theorem 4.2]).

######  Theorem 3.

If the network $(V,E)$ is strongly connected then its associated mass-action
system is permanent.

We also use special cases of [20, Theorems 4.1 and 4.2] to prove permanence.
Though the results in [20] are phrased for discrete-time systems, they also
hold in the continuous-time case with almost identical proofs. Let $\Delta$ be
the simplex $\\{x\in\mathbb{R}^{n}_{\geq 0}\colon x_{1}+\cdots+x_{n}=c\\}$ for
some $c>0$, and $F=(F^{t})_{t\geq 0}$ a semi-flow on $\Delta$ such that the
relative interior $\Delta^{\circ}$ of $\Delta$ is forward invariant. Let $M$
be the maximal invariant subset (i.e., the union of all complete orbits) of
$\partial\Delta=\Delta\setminus\Delta^{\circ}$ under $F$. The _stable set_ of
a compact invariant set $K$ is defined as
$W^{s}(K)=\\{x\in\Delta\colon\omega(x)\subseteq K\\}$. We say that a compact
invariant set $K$ is _isolated_ in $\Delta$ if there is a neighborhood $U$ of
$K$ in $\Delta$ such that $K$ is the maximal invariant set in $U$.

######  Theorem 4.

If $M$ is isolated in $\Delta$ and $W^{s}(M)\subseteq\partial\Delta$ then $F$
is permanent.

######  Theorem 5.

Let $M_{1},M_{2}\subseteq M$ be invariant subsets of $M$, such that
$\alpha(x)\subseteq M_{1}$ and $\omega(x)\subseteq M_{2}$ hold for all $x\in
M\setminus(M_{1}\cup M_{2})$. If both $M_{1}$ and $M_{2}$ are isolated in
$\Delta$ and their stable sets satisfy $W^{s}(M_{i})\subseteq\partial\Delta$
then $F$ is permanent.

###  2.4. Molecularity

The _molecularity_ of a complex
$a_{1}\mathsf{X}_{1}+\cdots+a_{n}\mathsf{X}_{n}$ is the nonnegative integer
$a_{1}+\cdots+a_{n}$. We say that a reaction network is _bimolecular_ if every
complex in the network has molecularity at most two. Due to the following
theorem, rank-two bimolecular networks do not admit limit cycles, see [6,
Section 4].

######  Theorem 6.

Rank-two bimolecular mass-action systems do not admit limit cycles.

###  2.5. Homogenisation

We now recall a construction, _homogenisation_ , from [3, Section 4]. Start
with a mass-action system that admits $r\geq 1$ hyperbolic limit cycles and
add a new species $\mathsf{X}_{n+1}$ to each reaction in the network in a way
that the reactant and the product complex have the same molecularity in the
new network (notice that the new network as a directed graph might not be
isomorphic to the original one). Then, by [3, Theorem 1], the new mass-action
system admits at least $r$ hyperbolic limit cycles. Limit cycles that are born
via Andronov–Hopf bifurcation are hyperbolic (at least for parameters that are
close enough to the critical value).

###  2.6. Competitive systems

Finally, we recall some basic facts about competitive systems [18]. With
$(a_{1},b_{1})$, $(a_{2},b_{2})$, $(a_{3},b_{3})$ being (possibly unbounded)
open intervals in $\mathbb{R}$ let
$D=(a_{1},b_{1})\times(a_{2},b_{2})\times(a_{3},b_{3})$. Further, let $f\colon
D\to\mathbb{R}^{3}$ be analytic, and consider the differential equation

$\displaystyle\dot{x}=f(x).$ (1)

Let $J(x)$ denote the Jacobian matrix of the r.h.s. of (1) at $x\in D$. The
system (1) is called _competitive_ in $D$ if all the off-diagonal entries of
$J$ are nonpositive everywhere in $D$. The system (1) is called _irreducible_
in $D$ if $J$ is irreducible everywhere in $D$.

The most important result about three-dimensional competitive systems is the
Poincaré–Bendixson Theorem, due to Morris Hirsch and Hal Smith, see [18,
Theorem 3.23]. In particular, it excludes chaotic behaviour.

######  Theorem 7.

Assume that the differential equation (1) is competitive in $D$. Then every
nonempty bounded $\omega$-limit set either contains an equilibrium or is a
periodic orbit.

The following theorem is due to Hsiu-Rong Zhu and Hal Smith [30, Theorem 1.2].

######  Theorem 8.

Assume that the differential equation (1) has the following three properties.

  1. (A)

It is permanent in $D$.

  2. (B)

It is competitive and irreducible in $D$.

  3. (C)

The domain $D$ contains a unique equilibrium $x^{*}$, the determinant of
$J(x^{*})$ is negative, and $x^{*}$ is unstable.

Then there is at least one but no more than finitely many periodic orbits and
at least one of them is orbitally asymptotically stable.

Throughout the paper, whenever we write _stable limit cycle_ , we mean it is
orbitally asymptotically stable.

##  3\. Feinberg–Berner oscillator

In 1979, Martin Feinberg introduced the mass-action system

in [9, Example 3.D.3] as a perturbation of the Edelstein mass-action system
[8]. The Edelstein network is obtained from (3) by deleting the reversible
reaction $\mathsf{X}\rightleftarrows\mathsf{Z}$ (which amounts to setting
$\kappa_{3}=\kappa_{4}=0$ in the differential equation). By the Deficiency-One
Theorem (Theorem 2 in Section 2), the differential equation (3) has exactly
one positive equilibrium for all choices of rate constants. This holds despite
the fact that the Edelstein mass-action system does admit multiple positive
equilibria (for some rate constants in some stoichiometric class).

The additional purpose of Feinberg for presenting example (3) is to
demonstrate that not all of the conclusions of the Deficiency-Zero Theorem
hold true for networks under the scope of the Deficiency-One Theorem. In
short, though the existence and the uniqueness of the positive equilibrium
hold true, its asymptotic stability in general does not follow. According to
[9, Remark 3.5], Paul Berner identified rate constants for which the unique
positive equilibrium of (3) is unstable and numerical solutions suggest the
existence of a stable limit cycle. The same remark is made in [10, page 108]
and [11, page 66]. The mass-action system (3) appears further in [12, Remark
6.2.B] and [14, (4.12)]. To our knowledge, the only place, where specific rate
constants are given for which the unique positive equilibrium is claimed to be
unstable is [16, (8.5)]. However, unfortunately, that claim is not correct:
for the given rate constants all the three eigenvalues of the Jacobian matrix
at the unique positive equilibrium are real and negative, and thus, the
equilibrium is linearly stable.

Below we prove that there exist rate constants such that the unique positive
equilibrium of the mass-action system (3) is unstable and is surrounded by a
stable limit cycle. Furthermore, we also prove that there exist rate constants
such that the positive equilibrium is asymptotically stable and is surrounded
by two limit cycles (an unstable and a stable one).

Before we turn to Theorem 10, the analysis of the mass-action system (3), we
discuss in what sense the network is minimal. A reaction network is called
_reversible_ if along with any reaction $i\to j$, the reverse reaction $j\to
i$ is also present in the network. The network (3) is reversible and it has
$5$ complexes and $4$ reversible reactions. The following proposition states
that neither the number of complexes nor the number of reactions could be less
for a reversible bimolecular network whose associated mass-action system
admits a limit cycle.

######  Proposition 9.

Any reversible bimolecular network whose associated mass-action system admits
a limit cycle has

  1. (i)

at least $5$ complexes and

  2. (ii)

at least $4$ reversible reactions.

###### Proof.

By Theorem 6, the rank of a bimolecular mass-action system with a limit cycle
is at least $3$. For the rank to be at least $3$, it must have at least $4$
complexes. However, with only $4$ complexes, the deficiency is $0$, and that
precludes the possibility of a limit cycle by Theorem 1. This concludes the
proof of (i). If the network has only $3$ reversible reactions, and the rank
of the network is $3$ then it is not hard to see that, again, the deficiency
is zero, which, again, precludes the possibility of a limit cycle. This
concludes the proof of (ii). ∎

######  Theorem 10.

For the mass-action system (3) the following statements hold.

  1. (i)

The system is permanent.

  2. (ii)

The system is competitive.

  3. (iii)

If the unique positive equilibrium is unstable then there exists a stable
limit cycle and there are only finitely many periodic orbits.

  4. (iv)

The unique positive equilibrium $(x^{*},y^{*},z^{*})$ equals $(1,1,1)$ if and
only if

$\displaystyle\kappa_{3}=\kappa_{4}\text{ and
}\kappa_{1}-\kappa_{2}=\kappa_{5}-\kappa_{6}=\kappa_{7}-\kappa_{8}.$ (2)

  5. (v)

If (2) and $\kappa_{2}+2\kappa_{4}\geq\kappa_{1}$ hold then $(1,1,1)$ is
linearly stable.

  6. (vi)

Set $\kappa_{1}=1$, $\kappa_{2}=\frac{1}{5}$, $\kappa_{4}=\frac{1}{5}$ and
eliminate $\kappa_{3}$, $\kappa_{5}$, $\kappa_{7}$ by (2) (thus, the only free
parameters are $\kappa_{6}>0$ and $\kappa_{8}>0$).

    1. (a)

Then the positive equilibrium $(1,1,1)$ is asymptotically stable on
$\mathcal{H}_{s}$, undergoes an Andronov–Hopf bifurcation at
$\mathcal{H}_{c}$, and is unstable on $\mathcal{H}_{u}$, where

$\displaystyle\mathcal{H}_{s}$
$\displaystyle=\\{(\kappa_{6},\kappa_{8})\in\mathbb{R}^{2}_{+}\colon
h(\kappa_{6},\kappa_{8})>0\\},$ $\displaystyle\mathcal{H}_{c}$
$\displaystyle=\\{(\kappa_{6},\kappa_{8})\in\mathbb{R}^{2}_{+}\colon
h(\kappa_{6},\kappa_{8})=0\\},$ $\displaystyle\mathcal{H}_{u}$
$\displaystyle=\\{(\kappa_{6},\kappa_{8})\in\mathbb{R}^{2}_{+}\colon
h(\kappa_{6},\kappa_{8})<0\\}$

with
$h(\kappa_{6},\kappa_{8})=50\kappa_{6}^{2}\kappa_{8}+100\kappa_{6}\kappa_{8}^{2}+55\kappa_{6}^{2}+260\kappa_{6}\kappa_{8}+50\kappa_{8}^{2}+128\kappa_{6}+40\kappa_{8}-26$.

    2. (b)

Along the curve $\mathcal{H}_{c}$, the first focal value is negative for
$\kappa_{6}$ close to $0$, while it is positive for $\kappa_{8}$ close to $0$.
Thus, the Andronov–Hopf bifurcation could be supercritical, subcritical, and
degenerate.

  7. (vii)

There exist rate constants such that the unique positive equilibrium is
unstable and is surrounded by a stable limit cycle.

  8. (viii)

There exist rate constants such that the unique positive equilibrium is
asymptotically stable and is surrounded by two limit cycles (an unstable and a
stable one).

###### Proof.

Since the network is strongly connected, it is permanent by Theorem 3, thus,
(i) follows.

With $\widetilde{z}=-z$, the differential equation (3) in the new coordinates
$(x,y,\widetilde{z})$ take the form

$\displaystyle\begin{split}\dot{x}&=-\kappa_{1}xy-\kappa_{2}\widetilde{z}-\kappa_{3}\widetilde{z}-\kappa_{4}x+\kappa_{5}x-\kappa_{6}x^{2},\\\
\dot{y}&=-\kappa_{1}xy-\kappa_{2}\widetilde{z}-\kappa_{7}\widetilde{z}-\kappa_{8}y,\\\
\dot{\widetilde{z}}&=-\kappa_{1}xy-\kappa_{2}\widetilde{z}-\kappa_{3}\widetilde{z}-\kappa_{4}x-\kappa_{7}\widetilde{z}-\kappa_{8}y\end{split}$
(3)

with state space $\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}_{-}$. The
Jacobian matrix of the r.h.s. equals

$\displaystyle\begin{bmatrix}-\kappa_{1}y-\kappa_{4}+\kappa_{5}-2\kappa_{6}x&-\kappa_{1}x&-\kappa_{2}-\kappa_{3}\\\
-\kappa_{1}y&-\kappa_{1}x-\kappa_{8}&-\kappa_{2}-\kappa_{7}\\\
-\kappa_{1}y-\kappa_{4}&-\kappa_{1}x-\kappa_{8}&-\kappa_{2}-\kappa_{3}-\kappa_{7}\\\
\end{bmatrix},$ (4)

a matrix whose off-diagonal entries are negative on
$\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}_{-}$, the system is thus
competitive, proving (ii).

Statement (iii) follows by the application of Theorem 8. The determinant of
the Jacobian matrix at the unique positive equilibrium is negative by [6,
Theorem 5(iii)(b)].

Statement (iv) follows by a short calculation.

Next, we prove (v). We eliminate $\kappa_{3}$, $\kappa_{5}$, $\kappa_{7}$ by
(2). Thus, the remaining parameters are $\kappa_{1}$, $\kappa_{2}$,
$\kappa_{4}$, $\kappa_{6}$, $\kappa_{8}$ with all of them being positive and
additionally both of $\kappa_{6}+\kappa_{1}-\kappa_{2}$ and
$\kappa_{8}+\kappa_{1}-\kappa_{2}$ are also positive. The Jacobian matrix at
the unique positive equilibrium $(1,1,1)$ equals

$\displaystyle\begin{bmatrix}-\kappa_{2}-\kappa_{4}-\kappa_{6}&-\kappa_{1}&\kappa_{2}+\kappa_{4}\\\
-\kappa_{1}&-\kappa_{1}-\kappa_{8}&\kappa_{1}+\kappa_{8}\\\
\kappa_{1}+\kappa_{4}&\kappa_{1}+\kappa_{8}&-\kappa_{1}-\kappa_{4}-\kappa_{8}\end{bmatrix}.$

The characteristic polynomial of this matrix equals
$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$, where

$\displaystyle a_{2}$
$\displaystyle=2\kappa_{1}+\kappa_{2}+2\kappa_{4}+\kappa_{6}+2\kappa_{8},$
$\displaystyle a_{1}$
$\displaystyle=\kappa_{1}(\kappa_{2}+2\kappa_{4}-\kappa_{1})+2\kappa_{2}\kappa_{8}+\kappa_{4}\kappa_{6}+3\kappa_{4}\kappa_{8}+2\kappa_{1}\kappa_{6}+2\kappa_{6}\kappa_{8},$
$\displaystyle a_{0}$
$\displaystyle=\kappa_{4}(\kappa_{1}\kappa_{6}+\kappa_{1}\kappa_{8}+\kappa_{6}\kappa_{8}).$

By the Routh-Hurwitz criterion, all the three eigenvalues have negative real
part if and only if $a_{0}$, $a_{2}$, $a_{2}a_{1}-a_{0}$ are all positive.
Clearly, $a_{0}$ and $a_{2}$ are positive. Under
$\kappa_{2}+2\kappa_{4}-\kappa_{1}\geq 0$, the product $a_{2}a_{1}$ is
positive as well. The positivity of $a_{2}a_{1}-a_{0}$ then follows, because
the negative terms $-\kappa_{1}\kappa_{4}\kappa_{6}$,
$-\kappa_{1}\kappa_{4}\kappa_{8}$, $-\kappa_{4}\kappa_{6}\kappa_{8}$ coming
from $-a_{0}$ are compensated by the positive terms
$6\kappa_{1}\kappa_{4}\kappa_{6}$, $6\kappa_{1}\kappa_{4}\kappa_{8}$,
$9\kappa_{4}\kappa_{6}\kappa_{8}$ coming from $a_{2}a_{1}$. This concludes the
proof of (v).

To prove (vi), notice that in the special case $\kappa_{1}=1$,
$\kappa_{2}=\frac{1}{5}$, $\kappa_{4}=\frac{1}{5}$, the expression
$a_{2}a_{1}-a_{0}$ simplifies to $h(\kappa_{6},\kappa_{8})$ (up to a positive
multiplier). Part (a) in (vi) then follows immediately by the Routh-Hurwitz
criterion. For proving part (b), one computes the first focal value, $L_{1}$,
along the curve $\mathcal{H}_{c}$. This can be performed e.g. by applying
formula (5.39) in [22] using a software that is capable for symbolic
computations, like Mathematica. Unfortunately, in this particular case, the
obtained expression for $L_{1}$ is complicated. One can nevertheless analyse
it by Mathematica and find that $L_{1}<0$ when $\kappa_{6}$ is close to $0$,
while $L_{1}>0$ when $\kappa_{8}$ is close to $0$. We have depicted in the
left panel of Figure 1 the sign of $h(\kappa_{6},\kappa_{8})$ and the sign of
the first focal value. This concludes the proof of (vi).

To prove (vii), it suffices to note that a supercritical Andronov–Hopf
bifurcation is possible by (vi).

To prove (viii), fix rate constants such that the unique positive equilibrium
has a negative real eigenvalue, a pair of purely imaginary eigenvalues, and
the first focal value is positive. Such rate constants exist by (vi), see the
orange curve in the left panel of Figure 1. Since the equilibrium is unstable,
there exists a stable limit cycle by (iii). Now perturb the rate constants
slightly and create an unstable limit cycle on the center manifold via a
subcritical Andronov–Hopf bifurcation. The equilibrium then becomes
asymptotically stable. We illustrated the two limit cycles in the right panel
of Figure 1. This concludes the proof of (viii).

|
---|---
Figure 1. Left panel: the bifurcation diagram of the mass-action system (3)
for $\kappa_{1}=1$, $\kappa_{2}=\frac{1}{5}$, $\kappa_{3}=\frac{1}{5}$,
$\kappa_{4}=\frac{1}{5}$, $\kappa_{5}=\kappa_{6}+\frac{4}{5}$,
$\kappa_{7}=\kappa_{8}+\frac{4}{5}$, while $\kappa_{6}$ and $\kappa_{8}$ are
kept as parameters. Right panel: illustrated are the two limit cycles of the
mass-action system (3) with $\kappa_{6}=0.187$ and $\kappa_{8}=0.0052$ (notice
that these parameters fall slightly above the orange curve of the left panel).

We now show an alternative way to prove (viii) that avoids the application of
(iii) and furthermore guarantees that the two limit cycles are hyperbolic. Fix
rate constants such that the unique positive equilibrium has a negative real
eigenvalue, a pair of purely imaginary eigenvalues, the first focal value
vanishes, and the second focal value is negative. One finds such rate
constants numerically e.g. by MATCONT [7]: for instance, along the curve
$\mathcal{H}_{c}$ in (vi) there is one such point (it is shown in blue on the
left panel in Figure 1). Then perturb the rate constants $\kappa_{6}$ and
$\kappa_{8}$ slightly along $\mathcal{H}_{c}$ to make the first focal value
positive. Then a stable limit cycle is born via a degenerate Andronov–Hopf
bifurcation (or Bautin bifurcation) and the equilibrium becomes unstable.
Finally, perturb the rate constants $\kappa_{6}$ and $\kappa_{8}$ away from
$\mathcal{H}_{c}$ to $\mathcal{H}_{s}$. With this, an unstable limit cycle is
born via a subcritical Andronov–Hopf bifurcation and the equilibrium becomes
asymptotically stable. ∎

We remark that the choice $\kappa_{1}=1$, $\kappa_{2}=\frac{1}{5}$,
$\kappa_{4}=\frac{1}{5}$ in part (vi) of Theorem 10 does not play a crucial
role in the sense that for any $\kappa_{1}$, $\kappa_{2}$, $\kappa_{4}$ with
$\kappa_{2}+2\kappa_{4}<\kappa_{1}$ we have the same qualitative picture:
there is a curve in the positive quadrant of the
$(\kappa_{6},\kappa_{8})$-plane that connects the two axes, and where the
unique positive equilibrium undergoes an Andronov–Hopf bifurcation, which is
supercritical close to the $\kappa_{8}$-axis and subcritical close to the
$\kappa_{6}$-axis.

Before we turn to the homogenisation of the Feinberg–Berner oscillator (3), we
briefly discuss detailed balance and complex balance [21, 29]. For a
reversible mass-action system, a positive vector $x$ is _detailed balanced_ if
for each reversible reaction the rate of the forward and the backward reaction
at $x$ equals. For a mass-action system, a positive vector $x$ is _complex
balanced_ if at every complex $i$ the total rate at $x$ over the reactions
whose reactant complex is $i$ equals the total rate at $x$ over the reactions
whose product complex is $i$. Clearly, a detailed balanced vector is also
complex balanced. Further, a detailed balanced or a complex balanced vector is
an equilibrium. Finally, if a mass-action system admits a detailed balanced
(respectively, complex balanced) equilibrium then all equilibria are detailed
balanced (respectively, complex balanced). The following proposition tells us
about detailed balance and complex balance in the mass-action system (3).

######  Proposition 11.

For the mass-action system (3) the following are equivalent.

  1. (a)

The system is complex balanced,

  2. (b)

The system is detailed balanced.

  3. (c)

The relation
$\frac{\kappa_{1}}{\kappa_{2}}\frac{\kappa_{5}}{\kappa_{6}}\frac{\kappa_{7}}{\kappa_{8}}=1$
holds.

###### Proof.

Note that if we replace each reversible reaction with an undirected edge in
(3) then we get a tree. For such reversible networks, complex balance and
detailed balance are equivalent. By definition, the mass-action system (3) is
detailed balanced if there exists $(x,y,z)\in\mathbb{R}^{3}_{+}$ such that

$\displaystyle\begin{split}\kappa_{1}xy&=\kappa_{2}z,\\\
\kappa_{3}z&=\kappa_{4}x,\\\ \kappa_{5}x&=\kappa_{6}x^{2},\\\
\kappa_{7}z&=\kappa_{8}y.\\\ \end{split}$ (5)

Solving the last three equations in (5) for $x$, $y$, $z$ yields
$x=\frac{\kappa_{5}}{\kappa_{6}}$,
$y=\frac{\kappa_{4}}{\kappa_{3}}\frac{\kappa_{5}}{\kappa_{6}}$,
$z=\frac{\kappa_{4}}{\kappa_{3}}\frac{\kappa_{5}}{\kappa_{6}}\frac{\kappa_{7}}{\kappa_{8}}$.
It is compatible with the first equation in (5) if and only if
$\frac{\kappa_{1}}{\kappa_{2}}\frac{\kappa_{5}}{\kappa_{6}}\frac{\kappa_{7}}{\kappa_{8}}=1$
holds. This concludes the proof. ∎

We now homogenise the network in (3) and obtain the mass-conserving mass-
action system Since the network is reversible, there exists a positive
equilibrium in every stoichiometric class $x+y+z+w=c$ with $c>0$ [4, Theorem
1]. The uniqueness of the positive equilibrium in each of these stoichiometric
classes follows by application of the Deficiency-One Algorithm [13, 15]. In
fact, since it is possible to parametrise the set of positive equilibria, one
sees directly the existence and uniqueness. In the following proposition we
collected a few properties of the mass-action system (LABEL:eq:F_homog_ode).

######  Theorem 12.

For the mass-action system (LABEL:eq:F_homog_ode) the following statements
hold.

  1. (i)

The set of positive equilibria is
$\\{(x^{*}(t),y^{*}(t),z^{*}(t),w^{*}(t))\colon t>0\\}$, where

$\displaystyle x^{*}(t)$ $\displaystyle=t,$ $\displaystyle y^{*}(t)$
$\displaystyle=\frac{\kappa_{2}\kappa_{3}\kappa_{6}t+\kappa_{7}(\kappa_{2}\kappa_{4}+\kappa_{3}\kappa_{5})}{\kappa_{1}\kappa_{3}\kappa_{5}t+\kappa_{8}(\kappa_{2}\kappa_{4}+\kappa_{3}\kappa_{5})}\frac{\kappa_{4}}{\kappa_{3}}t,$
$\displaystyle z^{*}(t)$ $\displaystyle=\frac{\kappa_{4}}{\kappa_{3}}t,$
$\displaystyle w^{*}(t)$
$\displaystyle=\frac{\kappa_{1}\kappa_{3}\kappa_{6}t+(\kappa_{1}\kappa_{4}\kappa_{7}+\kappa_{3}\kappa_{6}\kappa_{8})}{\kappa_{1}\kappa_{3}\kappa_{5}t+\kappa_{8}(\kappa_{2}\kappa_{4}+\kappa_{3}\kappa_{5})}t.$

  2. (ii)

There exists a unique positive equilibrium in every stoichiometric class
$x+y+z+w=c$ with $c>0$.

  3. (iii)

The following are equivalent.

    1. (a)

The system is complex balanced.

    2. (b)

The system is detailed balanced.

    3. (c)

The relation
$\frac{\kappa_{1}}{\kappa_{2}}\frac{\kappa_{5}}{\kappa_{6}}\frac{\kappa_{7}}{\kappa_{8}}=1$
holds.

    4. (d)

The set of positive equilibria is
$\\{(at^{\alpha},bt^{\beta},ct^{\gamma},dt^{\delta})\colon t>0\\}$ with some
$a$, $b$, $c$, $d$ positive numbers and $\alpha$, $\beta$, $\gamma$, $\delta$
real numbers.

    5. (e)

The formula $y^{*}(t)$ in (i) above is linear in $t$.

    6. (f)

The formula $w^{*}(t)$ in (i) above is linear in $t$.

  4. (iv)

There exist rate constants and a stoichiometric class with an asymptotically
stable equilibrium and two limit cycles (an unstable and a stable one).

  5. (v)

The system is permanent.

###### Proof.

One proves (i) by a direct calculation.

Since for any $a$, $b$, $c$, $d$ positive numbers the function
$t\mapsto\frac{at^{2}+bt}{ct+d}$ is strictly increasing on $(0,\infty)$,
statement (ii) follows immediately from the parametrisation in (i).

We now prove (iii). The equivalence of (a), (b), (c) can be shown in the same
way as in Proposition 11. The equivalence of (c), (d), (e), (f) follows by a
direct computation, using the parametrisation in (i).

Statement (iv) is an immediate consequence of part (viii) in Theorem 10 and
[3, Theorem 1].

To prove (v), we use Theorem 4. For any fixed $c>0$ consider the mass-action
system (LABEL:eq:F_homog_ode) in the tetrahedron

$\displaystyle\Delta=\\{(x,y,z,w)\in\mathbb{R}^{4}_{\geq 0}\colon
x+y+z+w=c\\}.$

The only equilibrium on the boundary of $\Delta$ is
$E_{\mathsf{W}}=(0,0,0,c)$. If we eliminate $w$ by $w=c-x-y-z$, the $3\times
3$ Jacobian matrix at the equilibrium $E_{\mathsf{W}}$ is given by

$\displaystyle
J=\begin{bmatrix}-\kappa_{4}+c\kappa_{5}&0&c\kappa_{2}+\kappa_{3}\\\
0&-\kappa_{8}&c\kappa_{2}+\kappa_{7}\\\
\kappa_{4}&\kappa_{8}&-c\kappa_{2}-\kappa_{3}-\kappa_{7}\end{bmatrix}.$

The determinant of this matrix is
$c\kappa_{8}(\kappa_{2}\kappa_{4}+\kappa_{3}\kappa_{5})>0$. Also, $J$ is
quasipositive, i.e., its off-diagonal terms are nonnegative. Hence, the
stability modulus of $J$ (i.e., its largest (or rightmost) eigenvalue) is
positive. This shows that the equilibrium $E_{\mathsf{W}}$ is not saturated
(see [19]), and therefore no solution from
$\Delta^{\circ}=\Delta\cap\mathbb{R}^{4}_{+}$ converges to $E_{\mathsf{W}}$,
and $\\{E_{\mathsf{W}}\\}$ is an isolated invariant set in $\Delta$. Since
$\\{E_{\mathsf{W}}\\}$ is the maximal invariant set within the boundary of
$\Delta$, Theorem 4 with $M=\\{E_{\mathsf{W}}\\}$ shows permanence. ∎

##  4\. Wilhelm–Heinrich oscillator

In 1995, Thomas Wilhelm and Reinhart Heinrich presented the smallest
bimolecular chemical reaction system with an Andronov–Hopf bifurcation [27].
Namely, they introduced the mass-action system It is minimal in the sense that
it has the smallest number of terms on the r.h.s. of the differential equation
(two monomials in each line; it is a $3$-dimensional S-system) and there is
only one quadratic term. In 1996, the authors studied the Andronov–Hopf
bifurcation [28] and in 2012, Hal Smith analysed the system further [25]. We
summarise their findings in the following theorem.

######  Theorem 13.

For the mass-action system (LABEL:eq:WH_original_ode) the following statements
hold.

  1. (i)

For $\kappa_{1}-\kappa_{4}>0$ there exists a unique positive equilibrium, it
is given by

$\displaystyle(x^{*},y^{*},z^{*})=\frac{\kappa_{1}-\kappa_{4}}{\kappa_{2}}\left(\frac{\kappa_{3}}{\kappa_{4}},1,\frac{\kappa_{3}}{\kappa_{5}}\right).$

  2. (ii)

The system is competitive.

  3. (iii)

All forward trajectories are bounded, and for $\kappa_{1}-\kappa_{4}>0$ the
system is permanent.

  4. (iv)

The unique positive equilibrium undergoes a supercritical Andronov–Hopf
bifurcation at $\kappa_{1}-\kappa_{4}=\kappa_{3}+\kappa_{5}$.

  5. (v)

The unique positive equilibrium is globally asymptotically stable for
$0<\kappa_{1}-\kappa_{4}\leq\kappa_{3}+\kappa_{5}$.

  6. (vi)

For $\kappa_{1}-\kappa_{4}>\kappa_{3}+\kappa_{5}$, the unique positive
equilibrium is unstable, it is surrounded by a stable limit cycle, and there
are only finitely many periodic orbits.

As remarked by Hal Smith, the main open problem about the dynamics of the
mass-action (LABEL:eq:WH_original_ode) is the number of limit cycles in case
$\kappa_{1}-\kappa_{4}>\kappa_{3}+\kappa_{5}$. He says, numerical simulations
suggest there is only one.

We now homogenise the network in (LABEL:eq:WH_original_ode) and obtain the
mass-conserving mass-action system By [3, Theorem 1], the mass-action system
(LABEL:eq:WH_homog_ode) admits a stable limit cycle, i.e., for some rate
constants and in some stoichiometric class there exists a stable limit cycle.
In fact, as per the following theorem, for some rate constants in some
stoichiometric class the unique positive equilibrium is asymptotically stable
and is surrounded by two limit cycles (an unstable and a stable one). Note
that this behaviour is precluded for the Wilhelm–Heinrich oscillator
(LABEL:eq:WH_original_ode).

We now analyse the homogenised Wilhelm–Heinrich oscillator
(LABEL:eq:WH_homog_ode). In particular, we show that the Andronov–Hopf
bifurcation could be supercritical, subcritical, or degenerate. W.l.o.g. we
assume $\kappa_{1}=1$ and for better readability we introduce $p=\kappa_{2}$,
$q=\kappa_{3}$, $r=\kappa_{4}$, $s=\kappa_{5}$. Thus, we study the
differential equation

$\displaystyle\begin{split}\dot{x}&=xw-rx-pxy,\\\ \dot{y}&=sz-qy,\\\
\dot{z}&=rx-sz,\\\ \dot{w}&=pxy+qy-xw.\end{split}$ (6)

######  Theorem 14.

For the mass-action system (6) the following statements hold.

  1. (i)

In the stoichiometric class $x+y+z+w=c$ there is no positive equilibrium if
$c\leq r$ and there is exactly one positive equilibrium if $c>r$. The set of
positive equilibria is given by

$\displaystyle(x^{*}(t),y^{*}(t),z^{*}(t),w^{*}(t))=t\left(1,\frac{r}{q},\frac{r}{s},\frac{pr}{q}\right)+(0,0,0,r)\text{
for }t>0.$ (7)

  2. (ii)

The positive equilibrium (7) is asymptotically stable on $\mathcal{H}_{s}$,
undergoes an Andronov–Hopf bifurcation at $\mathcal{H}_{c}$, and is unstable
on $\mathcal{H}_{u}$, where

$\displaystyle\mathcal{H}_{s}$
$\displaystyle=\\{(p,q,r,s,t)\in\mathbb{R}^{5}_{+}\colon h(p,q,r,s,t)>0\\},$
$\displaystyle\mathcal{H}_{c}$
$\displaystyle=\\{(p,q,r,s,t)\in\mathbb{R}^{5}_{+}\colon h(p,q,r,s,t)=0\\},$
$\displaystyle\mathcal{H}_{u}$
$\displaystyle=\\{(p,q,r,s,t)\in\mathbb{R}^{5}_{+}\colon h(p,q,r,s,t)<0\\}$

with $h(p,q,r,s,t)=(q+r+s)t^{2}+[(q+s)^{2}-prs]t+qs(q+s)$.

  3. (iii)

On $\mathcal{H}_{c}$, up to a positive factor, the first focal value, $L_{1}$,
equals

$\displaystyle(q+r+s)t^{2}-(q+s)(q+r+s)t-2qs(q+s).$ (8)

In particular, the Andronov–Hopf bifurcation could be supercritical,
subcritical, or degenerate.

  4. (iv)

For fixed $p,q,r,s>0$ there exists exactly $\begin{cases}0\\\ 1\\\
2\end{cases}t>0$ with

$\displaystyle(p,q,r,s,t)\in\mathcal{H}_{c}\text{ if }prs\begin{cases}<\\\
=\\\ >\end{cases}(q+s)^{2}+2\sqrt{sq(q+s)(q+r+s)}.$

  5. (v)

For $p,q,r,s>0$ with $prs>(q+s)^{2}+2\sqrt{sq(q+s)(q+r+s)}$, the first focal
value is negative at the smaller $t$ with $(p,q,r,s,t)\in\mathcal{H}_{c}$.

  6. (vi)

There exist rate constants and there exists a stoichiometric class such that
the positive equilibrium is asymptotically stable and is surrounded by two
limit cycles (an unstable and a stable one).

  7. (vii)

For $c>r$ the system is permanent in the stoichiometric class $x+y+z+w=c$.

###### Proof.

Statement (i) follows by a direct calculation.

To prove statement (ii), we eliminate $w$ using the conservation law
$x+y+z+w=c$, where $c=x^{*}(t)+y^{*}(t)+z^{*}(t)+w^{*}(t)$. Thus, we are left
with a differential equation in $3$ variables. The $3\times 3$ Jacobian
matrix, $J$, at $(x^{*}(t),y^{*}(t),z^{*}(t))$ is given by

$\displaystyle J=\begin{bmatrix}-t&-(p+1)t&-t\\\ 0&-q&s\\\
r&0&-s\end{bmatrix}.$

The characteristic polynomial of $J$ is

$\displaystyle\lambda^{3}+(q+s+t)\lambda^{2}+(qs+qt+rt+st)\lambda+(qr+qs+rs+prs)t.$

It has a real and negative root. By the Routh-Hurwitz criterion one finds that
the other two eigenvalues

  * •

have negative real part if $(p,q,r,s,t)\in\mathcal{H}_{s}$,

  * •

are purely imaginary if $(p,q,r,s,t)\in\mathcal{H}_{c}$,

  * •

have positive real part if $(p,q,r,s,t)\in\mathcal{H}_{u}$.

This concludes the proof of (ii).

To compute the first focal value on $\mathcal{H}_{c}$, we first bring the
system to canonical form. With

$\displaystyle\varrho$ $\displaystyle=-q-s-t,$ $\displaystyle\omega$
$\displaystyle=\sqrt{sq+qt+rt+st},$ $\displaystyle T$
$\displaystyle=\begin{bmatrix}0&s+t&s\\\
-\frac{rs}{\omega}&\frac{q(s+t)}{\omega}&-\frac{st}{\omega}\\\
rs&q(q+s)+t(q+r+s)&-s(q+s)\end{bmatrix}$

one confirms that

$\displaystyle TJT^{-1}=\begin{bmatrix}0&-\omega&0\\\ \omega&0&0\\\
0&0&\varrho\end{bmatrix}\text{ for }(p,q,r,s,t)\in\mathcal{H}_{c}.$

Application of formula (12) to the transformed system then provides the first
focal value. Omitting all the factors that are positive for all
$(p,q,r,s,t)\in\mathbb{R}^{5}_{+}$, we are left with the expression (8). With
$p=8$, $q=1$, $r=2$, $s=1$, $t=\frac{3\pm\sqrt{7}}{2}$ we have
$(p,q,r,s,t)\in\mathcal{H}_{c}$, and by (8) we have $L_{1}<0$ for
$t=\frac{3-\sqrt{7}}{2}$ and $L_{1}>0$ for $t=\frac{3+\sqrt{7}}{2}$.
Additionally, with $p=3(1+\sqrt{2})$, $q=1$, $r=2$, $s=1$, $t=1+\sqrt{2}$ we
have $(p,q,r,s,t)\in\mathcal{H}_{c}$ and $L_{1}=0$. Thus, the Andronov–Hopf
bifurcation could indeed be supercritical, subcritical, or degenerate. This
concludes the proof of (iii).

For fixed $p,q,r,s>0$ the discriminant of the quadratic function $t\mapsto
h(p,q,r,s,t)$ has the same sign as $prs-(q+s)^{2}-2\sqrt{sq(q+s)(q+r+s)}$.
This proves (iv).

To prove (v), we first introduce the parameters

$\displaystyle a$ $\displaystyle=q+r+s,$ $\displaystyle b$
$\displaystyle=prs-(q+s)^{2},$ $\displaystyle c$ $\displaystyle=qs(q+s),$
$\displaystyle d$ $\displaystyle=qs.$

With this, for fixed $p$, $q$, $r$, $s>0$, the quadratic function $t\mapsto
h(p,q,r,s,t)$ has exactly two positive roots provided that $b>0$ and
$b^{2}-4ac>0$. The smaller one of these two roots equals
$\frac{b-\sqrt{b^{2}-4ac}}{2a}$. In terms of the new parameters, formula (8)
becomes $at^{2}-\frac{ac}{d}t-2c$, which is negative for
$t<\frac{\frac{ac}{d}+\sqrt{\frac{a^{2}c^{2}}{d^{2}}+8ac}}{2a}$. Thus,
statement (v) follows once we show that

$\displaystyle
b-\sqrt{b^{2}-4ac}<\frac{ac}{d}+\sqrt{\frac{a^{2}c^{2}}{d^{2}}+8ac}$ (9)

holds for all $a>0$, $b>0$, $c>0$, $d>0$ with $b^{2}-4ac>0$. To show the
inequality (9), estimate the r.h.s. from below by $\sqrt{4ac}$ and add
$\sqrt{b^{2}-4ac}$ to both sides. Then we arrive at
$b<\sqrt{4ac}+\sqrt{b^{2}-4ac}$, which holds under the stated assumptions on
$a$, $b$, $c$, because $\sqrt{A^{2}+B^{2}}<A+B$ holds for all $A,B>0$.

To prove (vi), set $q=1$, $r=2$, $s=1$ and keep $p$ and $t$ as parameters.
Since $h(p,1,2,1,t)=-2(2t^{2}+(2-p)t+1)$, an Andronov–Hopf bifurcation happens
at the curve $2t^{2}+(2-p)t+1=0$ in the $(t,p)$-plane. Along this curve, the
first focal value vanishes at $(t^{*},p^{*})=(1+\sqrt{2},3(1+\sqrt{2}))$,
while it is negative (respectively, positive) for $t<t^{*}$ (respectively,
$t>t^{*}$), see the left panel in Figure 2. Using e.g. MATCONT [7], one
computes the second focal value at $(t^{*},p^{*})$ and finds it is negative.
Thus, the positive equilibrium is asymptotically stable. Perturb slightly
$(t,p)$ along the Andronov–Hopf curve in the direction $t>t^{*}$, with this
the first focal value becomes positive, the equilibrium is repelling on the
center manifold, and is surrounded by a stable limit cycle that is born via a
Bautin bifurcation. Finally, perturb away from the Andronov–Hopf curve in such
a way that the real parts of the nonreal eigenvalues at the equilibrium become
negative. With this, an unstable limit cycle is created via a subcritical
Andronov–Hopf bifurcation.

To prove (vii), we use Theorem 4. For any fixed $c>r$ consider the mass-action
system (6) in the tetrahedron

$\displaystyle\Delta=\\{(x,y,z,w)\in\mathbb{R}^{4}_{\geq 0}\colon
x+y+z+w=c\\}.$

The only equilibrium on the boundary of $\Delta$ is
$E_{\mathsf{W}}=(0,0,0,c)$. One finds that the eigenvalues of the Jacobian
matrix at $E_{\mathsf{W}}$ within $\Delta$ are $-q$, $-s$, $c-r$, two of them
are negative and one is positive for $c>r$ ($E_{\mathsf{W}}$ undergoes a
transcritical bifurcation at $c=r$, giving birth to the unique positive
equilibrium (7)). Notice that the face $\\{(x,y,z,w)\in\Delta\colon x=0\\}$ is
forward invariant and every solution starting on this face converges to
$E_{\mathsf{W}}$. Thus, this face is the stable manifold of $E_{\mathsf{W}}$,
and we have $W^{s}(E_{\mathsf{W}})\subseteq\partial\Delta$. Since the third
eigenvalue at $E_{\mathsf{W}}$ is positive, no solution starting in the
interior of $\Delta$ can converge to $E_{\mathsf{W}}$. By Theorem 4 with
$M=\\{E_{\mathsf{W}}\\}$, the system is then permanent. ∎

|
---|---
Figure 2. Left panel: the bifurcation diagram of the mass-action system (6)
for $q=1$, $r=2$, $s=1$, while $p$ and $t$ are kept as parameters. Right
panel: illustrated are the two limit cycles of the mass-action system (6) with
$p=8$ and $t=2.84$ (notice that these parameters fall slightly to the right of
the orange curve of the left panel).

##  5\. Wilhelm oscillator

In 2009, Thomas Wilhelm mentioned a bimolecular reaction network with $3$
species and $4$ reactions that admit a supercritical Andronov–Hopf
bifurcation, and thus, a stable limit cycle [26, Discussion]. Namely, he
introduced the mass-action system

Since two-species bimolecular mass-action systems do not admit a limit cycle
[23], one needs at least three species for the existence of a limit cycle. For
a three-species bimolecular mass-action system to admit a limit cycle, the
underlying network must have a three-dimensional stoichiometric subspace [24].
A three-species mass-action system with only three reactions and a three-
dimensional stoichiometric subspace is of deficiency-zero, and hence, by
Theorem 1, do not admit any periodic solution. Consequently, in this sense,
the bimolecular mass-action system (LABEL:eq:W_original_ode) is a minimal one
with a limit cycle.

We now analyse the Wilhelm oscillator (LABEL:eq:W_original_ode). In
particular, we show that the Andronov–Hopf bifurcation is always
supercritical.

######  Theorem 15.

For the mass-action system (LABEL:eq:W_original_ode) the following statements
hold.

  1. (i)

The system is competitive.

  2. (ii)

There exists a unique positive equilibrium, it is given by

$\displaystyle(x^{*},y^{*},z^{*})=\left(\sqrt{\frac{2\kappa_{4}}{\kappa_{2}}}\frac{\kappa_{1}}{\kappa_{3}},\frac{4\kappa_{1}\kappa_{4}}{\kappa_{3}^{2}},\frac{\kappa_{1}}{\kappa_{3}}\right).$

  3. (iii)

The unique positive equilibrium undergoes a supercritical Andronov–Hopf
bifurcation at $\kappa_{3}=4(\kappa_{4}+\sqrt{2\kappa_{2}\kappa_{4}})$, the
first focal value equals

$\displaystyle-\frac{r^{2}(3q^{4}+11q^{3}+17q^{2}+12q+4)}{q^{2}(q^{2}+5q+4)^{2}(q^{2}+8q+4)},$
(10)

where $q=\sqrt{\frac{2\kappa_{4}}{\kappa_{2}}}$ and
$r=\frac{\kappa_{3}}{\kappa_{1}}$.

  4. (iv)

The unique positive equilibrium is

$\begin{cases}\text{asymptotically stable}\\\ \text{unstable}\end{cases}$ for
$\kappa_{3}\begin{cases}\leq\\\
>\end{cases}4(\kappa_{4}+\sqrt{2\kappa_{2}\kappa_{4}})$.

  5. (v)

The unique positive equilibrium is surrounded by a stable limit cycle for
$\kappa_{3}$ slightly larger than
$4(\kappa_{4}+\sqrt{2\kappa_{2}\kappa_{4}})$.

###### Proof.

With $\widetilde{x}=-x$, the differential equation (LABEL:eq:W_original_ode)
in the new coordinates $(\widetilde{x},y,z)$ takes the form

$\displaystyle\dot{\widetilde{x}}$
$\displaystyle=2\kappa_{2}\widetilde{x}^{2}-\kappa_{3}yz,$
$\displaystyle\dot{y}$ $\displaystyle=\kappa_{1}y-\kappa_{3}yz,$
$\displaystyle\dot{z}$
$\displaystyle=\kappa_{2}\widetilde{x}^{2}-2\kappa_{4}z^{2}$

with state space $\mathbb{R}_{-}\times\mathbb{R}_{+}\times\mathbb{R}_{+}$. The
Jacobian matrix of the r.h.s. equals

$\displaystyle\begin{bmatrix}4\kappa_{2}\widetilde{x}&-\kappa_{3}z&-\kappa_{3}y\\\
0&\kappa_{1}-\kappa_{3}z&-\kappa_{3}y\\\
2\kappa_{2}\widetilde{x}&0&-4\kappa_{4}z\end{bmatrix},$

a matrix whose off-diagonal entries are nonpositive on
$\mathbb{R}_{-}\times\mathbb{R}_{+}\times\mathbb{R}_{+}$, the system is thus
competitive, proving (i).

Statement (ii) follows by direct calculation.

To prove (iii), (iv), and (v), we multiply the r.h.s. of the differential
equation (LABEL:eq:W_original_ode) by the positive constant
$\frac{\kappa_{3}}{\kappa_{1}\kappa_{2}}$ and introduce the new parameters
$p=\frac{\kappa_{3}}{\kappa_{2}}$, $q=\sqrt{\frac{2\kappa_{4}}{\kappa_{2}}}$,
and $r=\frac{\kappa_{3}}{\kappa_{1}}$. With this, the differential equation to
be studied is

$\displaystyle\dot{x}$ $\displaystyle=pryz-2rx^{2},$ $\displaystyle\dot{y}$
$\displaystyle=py-pryz,$ $\displaystyle\dot{z}$
$\displaystyle=rx^{2}-q^{2}rz^{2}$

with the unique positive equilibrium being
$(x^{*},y^{*},z^{*})=\left(\frac{q}{r},\frac{2q^{2}}{pr},\frac{1}{r}\right)$.
The Jacobian matrix, $J$, at $(x^{*},y^{*},z^{*})$ is given by

$\displaystyle J=\begin{bmatrix}-4q&p&2q^{2}\\\ 0&0&-2q^{2}\\\
2q&0&-2q^{2}\end{bmatrix}.$

The characteristic polynomial of $J$ equals
$\lambda^{3}+2q(q+2)\lambda^{2}+4q^{3}\lambda+4pq^{3}$. By the Routh-Hurwitz
criterion, all the three eigenvalues of $J$ have negative real part if and
only if $p<2q(q+2)$. For $p=2q(q+2)$, the eigenvalues are $\\{-p,\omega
i,-\omega i\\}$, where $\omega=\sqrt{\frac{\det
J}{\operatorname{tr}J}}=2q^{\frac{3}{2}}$. When $p>2q(q+2)$, one eigenvalue is
real and negative, the other two eigenvalues have positive real part. Since
$p\lesseqgtr 2q(q+2)$ if and only if $\kappa_{3}\lesseqgtr
4(\kappa_{4}+\sqrt{2\kappa_{2}\kappa_{4}})$, the statements (iv) and (v) will
follow once we prove that the first focal value is negative.

Thus, it is left to show that the first focal value equals (10). First, we
bring the Jacobian matrix to canonical form. Notice that when $p=2q(q+2)$,

$\displaystyle TJT^{-1}=\begin{bmatrix}0&-2q^{\frac{3}{2}}&0\\\
2q^{\frac{3}{2}}&0&0\\\ 0&0&-p\end{bmatrix},$

where the transformation matrix $T$ and its inverse, $T^{-1}$, are given by

$\displaystyle T=\begin{bmatrix}0&p&-2q^{2}\\\
2q^{\frac{3}{2}}&0&4q^{\frac{3}{2}}\\\ -1&1&q\end{bmatrix}\text{ and
}T^{-1}=\frac{1}{2q^{2}+p(q+2)}\begin{bmatrix}2&\frac{p+2q}{2\sqrt{q}}&-2p\\\
q+2&\sqrt{q}&2q^{2}\\\ -1&\frac{p}{2q^{\frac{3}{2}}}&p\end{bmatrix}.$

Application of the formula (12) gives that the first focal value indeed equals
(10). ∎

One way to construct a mass-conserving bimolecular mass-action system with a
limit cycle is to homogenise the Wilhelm oscillator. This leads to the mass-
action system

Since a limit cycle that is born via an Andronov–Hopf bifurcation is
hyperbolic (at least when the parameters are close enough to the bifurcation
point in parameter space), the homogenised Wilhelm oscillator (5) does admit a
stable limit cycle by [3, Theorem 1].

At the beginning of this section we explained in what sense the Wilhelm
oscillator (LABEL:eq:W_original_ode) is minimal. We now argue that the
homogenised Wilhelm oscillator (5) is minimal among the mass-conserving
bimolecular mass-action systems with a limit cycle. Since three-species
bimolecular mass-action systems with a two-dimensional stoichiometric subspace
do not admit a limit cycle [24], one needs at least four species. Since no
bimolecular mass-action system with a two-dimensional stoichiometric subspace
admits a limit cycle, the stoichiometric subspace has to be three-dimensional
at least. A four-species mass-action system with only three reactions and a
three-dimensional stoichiometric subspace is of deficiency-zero, and hence,
does not admit any periodic solution, see Theorem 1. Consequently, in this
sense, the mass-conserving bimolecular mass-action system (5) is a minimal one
with a limit cycle.

We now analyse the homogenised Wilhelm oscillator (5). In particular, we show
that the Andronov–Hopf bifurcation is always supercritical.

######  Theorem 16.

For the mass-action system (5) the following statements hold.

  1. (i)

The phase portrait in the stoichiometric classes $x+y+z+w=c_{1}$ and
$x+y+z+w=c_{2}$ are identical (up to a scaling of all the variables by
$\frac{c_{2}}{c_{1}}$).

  2. (ii)

In every positive stoichiometric class there exists a unique positive
equilibrium, the set of positive equilibria is given by

$\displaystyle(x^{*}(t),y^{*}(t),z^{*}(t),w^{*}(t))=t(q,2q^{2}/p,1,r)\text{
for }t>0,$

where $p=\frac{\kappa_{3}}{\kappa_{2}}$,
$q=\sqrt{\frac{2\kappa_{4}}{\kappa_{2}}}$, and
$r=\frac{\kappa_{3}}{\kappa_{1}}$.

  3. (iii)

The unique positive equilibrium undergoes a supercritical Andronov–Hopf
bifurcation at the surface

$\displaystyle\mathcal{H}_{c}=\left\\{(p,q,r)\in\mathbb{R}^{3}_{+}\colon
r>1,q<r(r-1),p=2q(q+2)\frac{r^{2}+(q+2)r+q}{r(r-1)-q}\right\\}.$

  4. (iv)

The unique positive equilibrium is asymptotically stable (respectively,
unstable) for $(p,q,r)\in\mathcal{H}_{s}\cup\mathcal{H}_{c}$ (respectively,
$(p,q,r)\in\mathcal{H}_{u}$), where

$\displaystyle\mathcal{H}_{u}$
$\displaystyle=\left\\{(p,q,r)\in\mathbb{R}^{3}_{+}\colon
r>1,q<r(r-1),p>2q(q+2)\frac{r^{2}+(q+2)r+q}{r(r-1)-q}\right\\},$
$\displaystyle\mathcal{H}_{s}$
$\displaystyle=\mathbb{R}^{3}_{+}\setminus(\mathcal{H}_{c}\cup\mathcal{H}_{u}).$

  5. (v)

There exists a stable limit cycle for parameters in $\mathcal{H}_{u}$ that are
close enough to $\mathcal{H}_{c}$.

  6. (vi)

The system is permanent.

###### Proof.

Statement (i) is an immediate consequence of the fact that the r.h.s. is a
homogeneous polynomial (of degree two).

Statement (ii) follows by a direct calculation.

To prove (iii), (iv), and (v), we multiply the r.h.s. of the differential
equation (5) by the positive constant
$\frac{\kappa_{3}}{\kappa_{1}\kappa_{2}}$ and write the differential equation
using the parameters $p$, $q$, $r$:

$\displaystyle\dot{x}$ $\displaystyle=pryz-2rx^{2},$ $\displaystyle\dot{y}$
$\displaystyle=pyw-pryz,$ $\displaystyle\dot{z}$
$\displaystyle=rx^{2}-q^{2}rz^{2},$ $\displaystyle\dot{w}$
$\displaystyle=rx^{2}+q^{2}rz^{2}-pyw.$

We pick the particular stoichiometric class that contains the equilibrium
$(x^{*}(t),y^{*}(t),z^{*}(t),w^{*}(t))$ with $t=1$ (by (i), this step does not
restrict generality). Next, we eliminate $w$ using the conservation law
$x+y+z+w=c$, where $c=x^{*}(1)+y^{*}(1)+z^{*}(1)+w^{*}(1)$. Thus, we are left
with a differential equation in $3$ variables. The $3\times 3$ Jacobian
matrix, $J$, at $(x^{*}(1),y^{*}(1),z^{*}(1))$ is given by

$\displaystyle J=\begin{bmatrix}-4qr&pr&2q^{2}r\\\
-2q^{2}&-2q^{2}&-2q^{2}(r+1)\\\ 2qr&0&-2q^{2}r\end{bmatrix}.$

The characteristic polynomial of $J$ equals
$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$ with

$\displaystyle a_{2}$ $\displaystyle=2q(q+2r+qr),$ $\displaystyle a_{1}$
$\displaystyle=2q^{2}r(p+2q(2+q+r)),$ $\displaystyle a_{0}$
$\displaystyle=4q^{3}r^{2}(2q^{2}+p(1+q+r)).$

One confirms by the Routh-Hurwitz criterion that all the three eigenvalues of
$J$ have negative real part if and only if $(p,q,r)\in\mathcal{H}_{s}$. The
eigenvalues are $\\{\omega i,-\omega i,\varrho\\}$ with
$\varrho=\operatorname{tr}J$ and $\omega=\sqrt{\frac{\det
J}{\operatorname{tr}J}}$ if and only if $(p,q,r)\in\mathcal{H}_{c}$. Finally,
one eigenvalue is real and negative, while the other two eigenvalues have
positive real part if and only if $(p,q,r)\in\mathcal{H}_{u}$. Once we prove
that the first focal value is negative, all of the statements (iii), (iv), and
(v) follow.

We now bring the Jacobian matrix to canonical form. With

$\displaystyle\varrho$ $\displaystyle=-2q(qr+q+2r),$ $\displaystyle\omega$
$\displaystyle=qr\sqrt{2\frac{2q^{2}+p(q+r+1)}{qr+q+2r}},$ $\displaystyle T$
$\displaystyle=\begin{bmatrix}r&0&q+2r\\\
-\frac{2q^{2}r}{\omega}&-\frac{pr^{2}}{\omega}&\frac{2q^{2}r(q+r)}{\omega}\\\
q-r(r-1)&q(r+1)+r(r+2)&q(1+r+r^{2})\end{bmatrix}$

one confirms that

$\displaystyle TJT^{-1}=\begin{bmatrix}0&-\omega&0\\\ \omega&0&0\\\
0&0&\varrho\end{bmatrix}\text{ for }(p,q,r)\in\mathcal{H}_{c}.$

Application of formula (12) to the transformed system then provides the first
focal value. Since it is a long expression, we omit it. That it is negative
everywhere on $\mathcal{H}_{c}$ is confirmed by a computer algebra software
like Mathematica.

To prove (vi), we use Theorem 5. For any fixed $c>0$ consider the mass-action
system (5) in the tetrahedron

$\displaystyle\Delta=\\{(x,y,z,w)\in\mathbb{R}^{4}_{\geq 0}\colon
x+y+z+w=c\\}.$

The only equilibria on the boundary of $\Delta$ are $E_{\mathsf{Y}}=(0,c,0,0)$
and $E_{\mathsf{W}}=(0,0,0,c)$. Let $M_{1}=\\{E_{\mathsf{W}}\\}$ and
$M_{2}=\\{E_{\mathsf{Y}}\\}$. The maximal invariant subset of $\partial\Delta$
is $M=\\{(x,y,z,w)\in\Delta\colon x=z=0\\}$, the edge connecting
$E_{\mathsf{Y}}$ and $E_{\mathsf{W}}$. The flow there goes from
$E_{\mathsf{W}}$ to $E_{\mathsf{Y}}$. At $E_{\mathsf{Y}}$, the eigenvalues
within $\Delta$ are $-c\kappa_{1}$, $0$, $0$. There is a one-dimensional
stable manifold (the edge $M\setminus\\{E_{\mathsf{W}}\\}$) and a two-
dimensional center manifold. Since
$\dot{x}+2\dot{z}=z(\kappa_{3}y-4\kappa_{4}z)$, which is positive in a
neighborhood of $E_{\mathsf{Y}}$, the flow goes away from $E_{\mathsf{Y}}$ on
this center manifold. At $E_{\mathsf{W}}$, the eigenvalues within $\Delta$ are
$c\kappa_{1}$, $0$, $0$. There is a one-dimensional unstable manifold (the
edge $M\setminus\\{E_{\mathsf{Y}}\\}$) and a two-dimensional center manifold
(the face $\\{(x,y,z,w)\in\Delta\colon y=0\\}$) which is the stable set
$W^{s}(E_{\mathsf{W}})\subseteq\partial\Delta$. Thus, Theorem 5 applies and
shows permanence. ∎

We conclude this section by collecting a few open questions:

  1. (i)

Is the mass-action system (LABEL:eq:W_original_ode) permanent? The difficulty
is that the $y$-axis is invariant and the flow there goes to infinity.

  2. (ii)

For the mass-action systems (LABEL:eq:W_original_ode) and (5), does local
stability of the positive equilibrium imply its global stability?

  3. (iii)

Is it true for the mass-action systems (LABEL:eq:W_original_ode) and (5) that
whenever a limit cycle exists, it is unique?

##  Appendix A Stability of an equilibrium

In this section we collect a few known results related to the stability of an
equilibrium of a three-dimensional differential equation, cf. [22, Section
5.4]. In particular, we present formula (12) below, which is used in the main
text a couple of times.

We start by recalling the Routh-Hurwitz criterion for a cubic polynomial
$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$:

  1. (i)

all three roots have negative real part if and only if $a_{0}>0$, $a_{2}>0$,
$a_{2}a_{1}>a_{0}$,

  2. (ii)

the three roots are $\omega i$, $-\omega i$, $\varrho$ with $\varrho<0$ and
$\omega>0$ if and only if $a_{0}>0$, $a_{2}>0$, $a_{2}a_{1}=a_{0}$.

We now sketch how to bring a matrix $J\in\mathbb{R}^{3\times 3}$ with
eigenvalues $\omega i$, $-\omega i$, $\varrho$ (where $\varrho$ and $\omega$
are nonzero real numbers) to canonical form. Since $\omega i$ is an eigenvalue
of $J$ (and of $J^{\top}$), the real number $-\omega^{2}$ is an eigenvalue of
$J^{2}$ (and of $(J^{\top})^{2}$). Let
$v_{1},v_{2},v_{3}\in\mathbb{R}^{3}\setminus\\{0\\}$ be such that

$\displaystyle(J^{\top})^{2}v_{1}$ $\displaystyle=-\omega^{2}v_{1},$
$\displaystyle J^{\top}v_{1}$ $\displaystyle=-\omega v_{2},$ $\displaystyle
J^{\top}v_{3}$ $\displaystyle=\varrho v_{3}.$

Let the $i$th row of the matrix $T\in\mathbb{R}^{3\times 3}$ be
$v_{i}^{\top}$. Then it is straightforward to check that

$\displaystyle TJT^{-1}=\begin{bmatrix}0&-\omega&0\\\ \omega&0&0\\\
0&0&\varrho\end{bmatrix}.$

With $\omega,\varrho\in\mathbb{R}\setminus\\{0\\}$, for the differential
equation

$\displaystyle\begin{split}\dot{x}&=-\omega
y+\frac{1}{2}\left[f_{200}x^{2}+2f_{110}xy+f_{020}y^{2}+2f_{101}xz+2f_{011}yz+\cdots\right]+\\\
&\quad+\frac{1}{6}\left[f_{300}x^{3}+3f_{120}xy^{2}+\cdots\right]+\text{h.o.t.},\\\
\dot{y}&=\omega
x+\frac{1}{2}\left[g_{200}x^{2}+2g_{110}xy+g_{020}y^{2}+2g_{101}xz+2g_{011}yz+\cdots\right]+\\\
&\quad+\frac{1}{6}\left[3g_{210}x^{2}y+g_{030}y^{3}+\cdots\right]+\text{h.o.t.},\\\
\dot{z}&=\varrho
z+\frac{1}{2}[h_{200}x^{2}+2h_{110}xy+h_{020}y^{2}+\cdots]+\text{h.o.t.},\end{split}$
(11)

the origin is an equilibrium with eigenvalues $\omega i$, $-\omega i$,
$\varrho$. The first focal value, $L_{1}$, there is given by

$\displaystyle\begin{split}&L_{1}=f_{300}+f_{120}+g_{030}+g_{210}+\\\
&\quad+\frac{1}{\omega}\left[f_{110}(f_{200}+f_{020})-g_{110}(g_{200}+g_{020})+f_{020}g_{020}-f_{200}g_{200}\right]-\\\
&\quad-\frac{h_{200}}{\varrho\left(\varrho^{2}+4\omega^{2}\right)}\left[(3\varrho^{2}+8\omega^{2})f_{101}-2\varrho\omega
f_{011}-2\varrho\omega g_{101}+(\varrho^{2}+8\omega^{2})g_{011}\right]-\\\
&\quad-\frac{2h_{110}}{\varrho^{2}+4\omega^{2}}\left[2\omega f_{101}+\varrho
f_{011}+\varrho g_{101}-2\omega g_{011}\right]-\\\
&\quad-\frac{h_{020}}{\varrho\left(\varrho^{2}+4\omega^{2}\right)}\left[(\varrho^{2}+8\omega^{2})f_{101}+2\varrho\omega
f_{011}+2\varrho\omega
g_{101}+(3\varrho^{2}+8\omega^{2})g_{011}\right].\end{split}$ (12)

We now briefly explain how to derive this formula. Let us consider the
quadratic approximation of the center manifold:

$\displaystyle z=\frac{1}{2}[c_{20}x^{2}+2c_{11}xy+c_{02}y^{2}]+\text{h.o.t.}$
(13)

To find $c_{20}$, $c_{11}$, $c_{02}$, differentiate both sides of (13) w.r.t.
time and thus find that

$\displaystyle\begin{split}\dot{z}&=c_{20}x\dot{x}+c_{11}(\dot{x}y+x\dot{y})+c_{02}y\dot{y}+\text{h.o.t.}=\\\
&=-\omega c_{20}xy+\omega c_{11}(-y^{2}+x^{2})+\omega
c_{02}xy+\text{h.o.t.}\end{split}$ (14)

On the other hand, from (11), we find that on the center manifold

$\displaystyle\dot{z}=\frac{1}{2}[(\varrho c_{20}+h_{200})x^{2}+2(\varrho
c_{11}+h_{110})xy+(\varrho c_{02}+h_{020})y^{2}]+\text{h.o.t.}$ (15)

Equating the coefficients of $x^{2}$, $xy$, $y^{2}$ in (14) and (15) gives a
system of three linear equations for the three unknowns $c_{20}$, $c_{11}$,
$c_{02}$. After solving that, plug in
$z=\frac{1}{2}[c_{20}x^{2}+2c_{11}xy+c_{02}y^{2}]$ to $\dot{x}$ and $\dot{y}$
in (11). This gives the dynamics on the two-dimensional center manifold and
one can then apply the planar first focal value formula to obtain (12).

## References

  * [1] D. F. Anderson, D. Cappelletti, J. Kim, and T. D. Nguyen.  Tier structure of strongly endotactic reaction networks.  Stochastic Processes and their Applications, 130(12):7218–7259, 2020.
  * [2] M. Banaji.  Splitting reactions preserves nondegenerate behaviours in chemical reaction networks, 2022.  https://arxiv.org/abs/2201.13105.pdf.
  * [3] M. Banaji, B. Boros, and J. Hofbauer.  Adding species to chemical reaction networks: preserving rank preserves nondegenerate behaviours, 2021.  https://arxiv.org/pdf/2112.06801.pdf.
  * [4] B. Boros.  Existence of positive steady states for weakly reversible mass-action systems.  SIAM Journal on Mathematical Analysis, 51(1):435–449, 2019.
  * [5] B. Boros and J. Hofbauer.  Permanence of weakly reversible mass-action systems with a single linkage class.  SIAM Journal on Applied Dynamical Systems, 19(1):352–365, 2020\.
  * [6] B. Boros and J. Hofbauer.  Limit cycles in mass-conserving deficiency-one mass-action systems, 2022\.  https://arxiv.org/abs/2202.10406.pdf.
  * [7] A. Dhooge, W. Govaerts, and Y. A. Kuznetsov.  MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs.  ACM Transactions on Mathematical Software, 29(2):141–164, 2003\.
  * [8] B. B. Edelstein.  Biochemical model with multiple steady states and hysteresis.  Journal of Theoretical Biology, 29(1):57–62, 1970.
  * [9] M. Feinberg.  Lectures on chemical reaction networks, 1979.  Notes of lectures given at the Mathematics Research Center, University of Wisconsin.
  * [10] M. Feinberg.  Chemical oscillations, multiple equilibria, and reaction network structure.  In Dynamics and Modelling of Reactive Systems, W. Stewart, W. Rey, and C. Conley, eds., pages 59–130, New York, NY, 1980. Academic Press.
  * [11] M. Feinberg.  Reaction network structure, multiple steady states, and sustained composition oscillations: A review of some results.  In Modelling of Chemical Reaction Systems, K. H. Ebert, P. Deuflhard P., W. Jäger, eds., Series in Chemical Physics, pages 56–68, Heidelberg, 1981. Springer.
  * [12] M. Feinberg.  Chemical reaction network structure and the stability of complex isothermal reactors - I. The Deficiency Zero and the Deficiency One Theorems.  Chemical Engineering Science, 42(10):2229–2268, 1987.
  * [13] M. Feinberg.  Chemical reaction network structure and the stability of complex isothermal reactors - II. Multiple steady states for networks of deficiency one.  Chemical Engineering Science, 43(1):1–25, 1988.
  * [14] M. Feinberg.  The existence and uniqueness of steady states for a class of chemical reaction networks.  Archive for Rational Mechanics and Analysis, 132(4):311–370, 1995\.
  * [15] M. Feinberg.  Multiple steady states for chemical reaction networks of deficiency one.  Archive for Rational Mechanics and Analysis, 132(4):371–406, 1995\.
  * [16] M. Feinberg.  Foundations of chemical reaction network theory, volume 202 of Applied Mathematical Sciences.  Springer, Cham, 2019.
  * [17] M. Gopalkrishnan, E. Miller, and A. Shiu.  A geometric approach to the global attractor conjecture.  SIAM Journal on Applied Dynamical Systems, 13(2):758–797, 2014\.
  * [18] M. W. Hirsch and H. Smith.  Monotone dynamical systems.  In Handbook of differential equations: ordinary differential equations, volume 2, pages 239–357. Elsevier, 2005.
  * [19] J. Hofbauer.  An index theorem for dissipative semiflows.  Rocky Mountain Journal of Mathematics, 20(4):1017–1031, 1990.
  * [20] J. Hofbauer and J. W.-H. So.  Uniform persistence and repellors for maps.  Proceedings of the American Mathematical Society, 107(4):1137–1142, 1989.
  * [21] F. Horn.  Necessary and sufficient conditions for complex balancing in chemical kinetics.  Archive for Rational Mechanics and Analysis, 49(3):172–186, 1972/73.
  * [22] Y. A. Kuznetsov.  Elements of Applied Bifurcation Theory, volume 112 of Applied Mathematical Sciences.  Springer-Verlag, New York, third edition, 2004.
  * [23] G. Póta.  Two-component bimolecular systems cannot have limit cycles: A complete proof.  Journal of Chemical Physics, 78(3):1621–1622, 1983.
  * [24] G. Póta.  Irregular behaviour of kinetic equations in closed chemical systems. oscillatory effects.  Journal of the Chemical Society, Faraday Transactions 2, 81(1):115–121, 1985.
  * [25] H. L. Smith.  Global dynamics of the smallest chemical reaction system with Hopf bifurcation.  Journal of Mathematical Chemistry, 50(4):989–995, 2012.
  * [26] T. Wilhelm.  The smallest chemical reaction system with bistability.  BMC Systems Biology, 3:90, 2009.
  * [27] T. Wilhelm and R. Heinrich.  Smallest chemical reaction system with Hopf bifurcation.  Journal of Mathematical Chemistry, 17:1–14, 1995.
  * [28] T. Wilhelm and R. Heinrich.  Mathematical analysis of the smallest chemical reaction system with Hopf bifurcation.  Journal of Mathematical Chemistry, 19(2):111–130, 1996.
  * [29] P. Y. Yu and G. Craciun.  Mathematical analysis of chemical reaction systems.  Israel Journal of Chemistry, 58(6–7):733–741, 2018.
  * [30] H.-R. Zhu and H. L. Smith.  Stable periodic orbits for a class of three dimensional competitive systems.  Journal of Differential Equations, 110(1):143–156, 1994.