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

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


Introduction and outline of the main result
An important theme in the mathematical study of chemical reaction networks (CRNs) relates to how network structure influences network dynamics. The results in this direction sometimes allow us to infer detailed information on dynamical behaviours of reaction networks using only graph theory and linear algebra. We may, for example, be able to conclude from basic computations that a CRN has very simple behaviour, such as the convergence of all initial conditions to an equilibrium, regardless of parameter values. In the opposite direction, we may be able to state, without numerical simulation, that a CRN admits some interesting behaviour such as stable oscillation, and even know a priori how to choose parameter values to obtain this behaviour.
Amongst more complicated behaviours which occur in ordinary differential equation models of CRNs, the most well-studied are multistationarity and oscillation. The history of study of these behaviours in the context of biological modelling is reviewed in [40]. Crucially, multistationarity and oscillation are not just of abstract interest, but may be of functional importance in biological switching and signalling processes [24,28,31,33,29,12,30]. For this reason, results which tell us which structures in a CRN guarantee such behaviours are of considerable interest.
A family of theorems termed "inheritance results" tell us when a CRN is guaranteed to exhibit some dynamical behaviour simply because of the presence of a certain subnetwork. Examples of such results can be found in [23,18,8,3,4]. Such conditions which guarantee nontrivial behaviours based on CRN structure are dual to claims which rule out nontrivial behaviours in a CRN. There is a large classical and modern literature on conditions which preclude certain dynamical behaviours in CRNs. Examples include [22,16,2,13,37,6,17,1]. Each new result in either direction helps to narrow the gap between conditions guaranteeing, and conditions ruling out, nontrivial behaviours in CRNs.
The main result of this paper is a new inheritance result, which is both natural and relatively straightforward to prove. The theorem is most simply stated in terms of the rank of a CRN. Each reaction of a CRN on n chemical species defines a reaction vector, a real (often integer) n-vector whose kth entry tells us the net production of species k in the reaction. The span of these vectors is the stoichiometric subspace of the CRN, whose dimension is defined to be the rank of the CRN. In any model of the CRN, numbers or concentrations of the chemical species are confined to affine subspaces parallel to the stoichiometric subspace. The nonnegative portions of these affine subspaces are termed the stoichiometric classes of the system. The rank of a CRN figures in various aspects of the theory. For example, it plays a crucial role in the original results of deficiency theory [16], and indeed in the definition of the deficiency of a reaction network. More recently, various results have been proved for CRNs of sufficiently low rank, regardless of how many species or reactions are involved. Characterisations of rank-1 CRNs admitting multistationarity are given in [27]. In rank-2 CRNs, the Poincaré-Bendixson theorem [41,Chapter 9] can be used to rule out chaos or guarantee oscillation; and the same can sometimes be extended to rank-3 CRNs, a fact exploited in the analysis in [39] and [11]. In [32], the famous "global attractor conjecture" is proved for rank-3 CRNs.
Stated informally, we have the following complementary inheritance results involving the rank of a CRN. Both are simple applications of regular perturbation theory: 1. Adding new reactions into a CRN without changing its rank preserves its capacity for nontrivial behaviours including nondegenerate multistationarity and oscillation. This was proved in previous work (Theorem 1 in [8] and Theorem 1 in [3]). 2. Adding new species into a CRN without changing its rank preserves its capacity for nontrivial behaviours including nondegenerate multistationarity and oscillation. This is the content of Theorem 1 here.
These claims demonstrate that in CRN theory we often have complementary results where we can interchange "species" and "reactions". Note, however, that adding species into a CRN while preserving its rank can result in some fairly fundamental changes to the CRN. For example, stoichiometric classes of the enlarged CRN may be bounded even if those of the original were unbounded. This is in contrast to adding linearly dependent reactions, a process which leaves stoichiometric classes unchanged.
Taken together, the above claims imply that building a CRN without altering its rank preserves its capacity for nondegenerate dynamical behaviours. This claim is stated more formally as Corollary 1 later. Note that the enlarged CRNs may, of course, admit more complicated and interesting behaviours than the original, as we shall see by example.
After some preliminary definitions, we will present the statement and proof of the main theorem and several remarks on its implications and generalisations. This is followed by several examples which demonstrate both the main result and its limitations.

Preliminaries
We present some key notions briefly. A much more expansive treatment of the main background can be found in previous work [7,8,3]. We consider a CRN to be an ordered set of chemical species and an ordered set of reactions, where the orderings are arbitrary but fixed. Each reaction is an ordered pair of complexes, namely formal linear combinations of chemical species. The coefficient of each species in a complex is taken to be nonnegative, and often to be an integer, although the latter is not required here. The pair of complexes which define a reaction are termed the "reactant complex" and "product complex" of the reaction.
Positive sets in Euclidean space. The positive orthant in R n is denoted by R n + and defined as {x ∈ R n : x i > 0 for i = 1, . . . , n}. A set in R n is termed positive if it lies in R n + . The closure of R n + is denoted by R n ≥0 and referred to as the nonnegative orthant in R n , namely R n ≥0 = {x ∈ R n : x i ≥ 0 for i = 1, . . . , n}.
Stoichiometric matrix and stoichiometric classes. Each reaction is associated with a reaction vector whose kth entry is the net production of the kth species in the reaction: this is just the stoichiometric coefficient of the kth species in the product complex minus its stoichiometric coefficient in the reactant complex. The stoichiometric matrix of a CRN is the matrix whose jth column is the jth reaction vector of the CRN. Given a CRN R on n species with stoichiometric matrix Γ, the span of the columns of Γ (a linear subspace of R n ) is denoted by im Γ, and is termed the stoichiometric subspace of R. The rank of Γ is termed the rank of R.
The nonnegative parts of cosets of im Γ, namely sets of the form (x + im Γ) ∩ R n ≥0 (x ∈ R n ≥0 ), are the stoichiometric classes of R. The positive parts of stoichiometric classes, namely sets of the form (x + im Γ) ∩ R n + (x ∈ R n + ), are the positive stoichiometric classes of R.
Ordinary differential equations (ODEs) and rate functions. Any system of ODEs describing the evolution of a CRN with stoichiometric matrix Γ takes the formẋ = Γv(x). The function v is termed the rate function of the CRN and its jth component tells us how the rate of the jth reaction depends on the concentrations of the chemical species. In this paper, v is assumed, as a minimum, to be defined and continously differentiable on R n + .
Nondegenerate and linearly stable limit sets. Consider a CRN of rank r with stoichiometric matrix Γ. Let S be some coset of im Γ containing a positive equilibrium (resp., periodic orbit) O. Since the positive part of S is locally invariant, O has exactly r eigenvalues (resp., Floquet multipliers) relative to S. If none of these are equal to zero (resp., exactly one of these is equal to one), then we say that O is nondegenerate. (Note that this terminology differs from that in [3], where an invariant set was referred to as "nondegenerate" only if it was hyperbolic relative to its stoichiometric class.) In an abuse of terminology, we refer to O as hyperbolic if it is hyperbolic relative to its stoichiometric class, and linearly stable if it is linearly stable relative to its stoichiometric class.
Continuation of nondegenerate limit sets. Suppose U is some open region in R n , a > 0, and f : U × (−a, a) → R n is C 1 . Ifẋ = f (x, 0) has a nondegenerate equilibrium (resp., periodic orbit) on U , then the same is true forẋ = f (x, ε) for all ε sufficiently small. The same conclusion holds if we replace "nondegenerate" by "hyperbolic" or "linearly stable". The results follow from the implicit function theorem and the fact that eigenvalues of a matrix depend continuously on its entries. The claim for equilibria is an immediate consequence of the implicit function theorem (see, for example, Corollary 6.9 in [8]), while the details for periodic orbits are laid out in Section IV in [20]. An immediate consequence is that nondegenerate equilibria and periodic orbits ofẋ = Γv(x) on some positive stoichiometric class survive sufficiently small C 1 perturbations to v(x) (here v is assumed to be at least C 1 ). In fact, if we restrict attention to hyperbolic equilibria and periodic orbits, then the results for equilibria and periodic orbits are special cases of more general results on the persistence of normally hyperbolic invariant manifolds in [19] and [21]. Thus, the main result here generalises naturally to the case of invariant manifolds which are normally hyperbolic relative to their stoichiometric class.
Entrywise products and generalised monomials. The notation a • b denotes the entrywise product of matrices or vectors a and b assumed to have the same dimensions. Given x = (x 1 , . . . , x n ) t and a = (a 1 , . . . , a n ), x a is an abbreviation for the generalised monomial x a1 1 x a2 2 · · · x an n . Let A 1 , . . . , A m be the rows of an m × n matrix A. Each x Ai is then a generalised monomial, and x A denotes the vector of these monomials, namely, (x A1 , x A2 , . . . , x Am ) t .
Kinetics and admitted behaviours. When we restrict the rate function of a CRN R to some class of functions K, the pair (R, K) is referred to as a CRN with kinetics. We can think of (R, K) as a set of allowed ODE models of the CRN. We say that (R, K) admits some particular dynamical behaviour if this behaviour occurs in some allowed model, i.e., for some choice of rate function from K and on some stoichiometric class. Otherwise the CRN with kinetics forbids this behaviour. Different classes of kinetics for CRNs are discussed in detail in [3].
Power-law kinetics. Let X 1 , . . . , X n denote the chemical species of a CRN and x 1 , . . . , x n denote the concentrations of these species. If the ith reaction has power-law kinetics, this means that the ith rate function takes the form v i (x) = κ i x a , where κ i is a positive constant (termed the rate constant of the ith reaction), and a is a real (row) vector, termed the vector of exponents for the reaction. If all reactions of a CRN have power-law kinetics, we can stack these row vectors into a matrix A, termed the matrix of exponents of the reaction network, whose ijth entry tells us the exponent of species j in the rate function for reaction i. In this case, the rate function can be written briefly κ • x A . If A is fixed in advance, we say that the CRN has fixed power-law kinetics.
Mass action kinetics. Mass action kinetics is a special case of fixed power-law kinetics where a ij , the ijth entry in the matrix of exponents, is precisely the stoichiometric coefficient of species X j in the reactant complex of reaction i.
Enlarging CRNs and inheritance. We are often interested in claiming that whenever a CRN with kinetics admits some behaviour, then so does an enlarged CRN with kinetics in some related class. Inheritance results relate dynamical behaviours admitted in the enlarged CRNs to those admitted in the original CRNs.
Adding linearly dependent species. In this paper we are interested in an enlargement where R is obtained from R by adding a linearly dependent species into the reactions of R. In such a modification, the stoichiometric matrix of the CRN is unchanged except for the addition of a new row, and this new row is a linear combination of the existing rows of the stoichiometric matrix. We refer to the process of enlarging a CRN R by adding in some new, linearly dependent, species as lifting. The terminology is motivated by the fact that the addition of a new species increases the dimension of the state space by 1, with the original, lower dimensional, state space naturally embedded in the new state space.
Derived power-law kinetics. Suppose that (R, K) and (R , K ) are two CRNs, with R obtained by adding new species and/or reactions to R, and K, K being fixed power-law kinetics with matrices of exponents A and A respectively. We then say that K is derived from K, if the submatrix of A corresponding to the original species and reactions of R is precisely A.
Permanence. Consider a system of ODEs on some subset of R n and suppose that X ⊆ R n ≥0 is forward invariant for the system. The system is permanent on X if there exists a forward invariant, compact, positive set Z ⊆ X such that the forward trajectory of every positive initial condition in X eventually enters Z. In the context of CRNs we may think of X as a stoichiometric class: we may be interested in permanence on some or all stoichiometric classes.

The main result
Given a parameterised family of compact sets X ε in Euclidean space, with ε ∈ (0, a) for some a > 0, "X ε is close to X " will mean that given any δ > 0, there exists ε 1 ∈ (0, a] such that for all ε ∈ (0, ε 1 ) the Hausdorff distance between X ε and X is less than δ. Theorem 1. Let (R, K) be a CRN with fixed power-law kinetics. Let R be derived by adding to R a new linearly dependent species, and let K be any fixed power-law kinetics for R derived from K. Suppose that for some choice of kinetics from K, R has, on some stoichiometric class, at least 0 ≤ r 1 < ∞ positive, nondegenerate (resp., hyperbolic, resp., linearly stable) equilibria and at least 0 ≤ r 2 < ∞ positive, nondegenerate (resp., hyperbolic, resp., linearly stable) periodic orbits. Then, for some choice of kinetics from K , R has, on some stoichiometric class, at least r 1 positive, nondegenerate (resp., hyperbolic, resp., linearly stable) equilibria and at least r 2 positive, nondegenerate (resp., hyperbolic, resp., linearly stable) periodic orbits.
Proof. Let R include n chemical species and m irreversible reactions with stoichiometric matrix Γ ∈ R n×m . By hypothesis, there exists v ∈ K, such that the associated ODE systeṁ has r 1 nondegenerate equilibria and r 2 nondegenerate periodic orbits on some positive stoichiometric class S + . The rate function v is analytic on R n + , and hence certainly C 1 on R n + which is, in fact, the only assumption we will need about v. Let S be the affine hull of S + , i.e., the coset of im Γ containing S + .
Let O refer to a nondegenerate positive equilibrium (resp., periodic orbit) on S. Choose Z ⊆ S to be compact and positive, with O ⊆ Z o , the relative interior of Z in S. We may assume that similar sets are constructed around each of the r 1 + r 2 nondegenerate equilibria and periodic orbits of (1), and that these sets are pairwise disjoint.
The hypothesis that the added species is linearly dependent implies that there exists c ∈ R n such that the stoichiometric matrix of R takes the form Define ε to be a positive parameter to be controlled, and let α j ∈ R (j = 1, . . . , m) be any real numbers. Denoting the concentration of the new species by y, set the rate of the jth reaction This choice of reaction rate corresponds to giving the new species exponent α j in reaction j, and multiplying the original rate constant of the jth reaction by ε αj . In brief notation, R now has rate function v (x, y, ε) Note that R has a new conservation law of the form −c t x + y = constant. For any fixed ε > 0, we will focus our attention on the invariant set H ε ⊆ R n+1 defined by setting this constant to be equal to 1 ε , namely, The map x is an affine bijection between R n and H ε . If c = 0, set ε 1 = 1, and otherwise set Then, for ε ∈ (0, ε 1 ), h ε (Z) is a compact, positive subset of H ε . The map h ε defines local coordinates on H ε which evolve according tȯ Here x refers to the local coordinate on H ε , rather than the original coordinate on R n : this should cause no confusion as we are identifying H ε with R n via h ε . Note that the right hand side of (3) is well-defined (and C 1 ) provided x is positive, and ε < 1 |c t x| , which certainly holds We wish to restrict our attention to Z o . We can, if desired, pass to local coordinates on S in a standard way (see the proofs of several results in [3]), but here this is unnecessary: we simply bear in mind that we are considering the restriction of (3) to Z o , with ε ∈ (−ε 1 , ε 1 ).
Since the vector field in (3) is a C 1 perturbation of that of (1), by regular perturbation theory, there exists ε 2 ∈ (0, ε 1 ] such that for each ε ∈ (0, ε 2 ), (3) has an equilibrium (resp., periodic orbit) O ε in Z o , which is nondegenerate, and close to O. (The details in the harder case where O is a periodic orbit are in Section IV in [20], for example.) If O is hyperbolic relative to S, then we can choose ε 2 to ensure that the linear stability type of O ε relative to S is the same as that of O. More precisely, (i) if O is an equilibrium with, relative to S, k 1 eigenvalues with positive real part, k 2 eigenvalues with negative real part, and no eigenvalues on the imaginary axis, then the same holds for O ε ; (ii) if O is a periodic orbit with, relative to S, k 1 Floquet multipliers inside the unit circle, k 2 outside the unit circle, and precisely one multiplier on the unit circle, then the same holds for O ε . As a special case, if O was linearly stable relative to S, then the same holds for O ε .
where x 0 is any element of S. Note that S ε has the same dimension as S. Clearly, is an equilibrium (resp., periodic orbit) of (2), and we have ensured (via the choice of ε 1 ) that As h ε is an affine bijection between S and S ε , the choice of ε 2 ensures that the linear stability type of O ε relative to S ε is the same as that of O ε relative to S.
We can repeat the same argument in a neighbourhood of each of a finite number of nondegenerate equilibria or periodic orbits of (1) on S + . By choosing ε * 2 to be the minimum of the values of ε 2 associated with each limit set, we can ensure that provided ε ∈ (0, ε * 2 ), R has at least r 1 positive, nondegenerate equilibria and at least r 2 positive, nondegenerate periodic orbits on S ε . Moreover, whenever one of the original limit sets was hyperbolic relative to S, we can ensure that the lifted limit set is of the same linear stability type relative to S ε . This completes the proof.
Several remarks are in order.
Remark 1 (Mass action kinetics). The result clearly holds if we insist that both R and R have mass action kinetics which is simply a special case of fixed power-law kinetics. In this case, in the proof of Theorem 1, α j is the stoichiometric coefficient of the new species in the reactant complex of the jth reaction of R .
Remark 2 (The proof is constructive). The proof of Theorem 1, as with other inheritance results based on perturbation theory, is constructive. It tells us how to set rate constants and how to choose a stoichiometric class in order to find the desired behaviour in the enlarged CRN R .
Remark 3 (The projected dynamics are close to the original). Consider some "lifted" bounded orbit of R such as O ε in the proof of Theorem 1. Its projection O ε onto x coordinates, can be made as close as we desire to the original orbit of R (namely, O) by choosing ε to be small. But this comes at the cost of large values of the new species concentration y on the lifted orbit, and small rate constants. The next remark indicates the limitations of the lifting process.
Remark 4 (We cannot always control the lifted dynamics over an entire stoichiometric class). The proof of Theorem 1 tells us the following: given any positive stoichiometric class, say S + , of the original CRN R, fixing the perturbation parameter ε at any positive value selects a positive stoichiometric class, say (S ε ) + , of the lifted CRN R . Assume that rate constants are fixed and let V and V ε refer to the original and lifted vector fields on S + and (S ε ) + respectively. Choosing ε to be small ensures that the projection of V ε onto S + is close to the original vector field V on S + on that portion of S ε where the concentration of the added species (denoted by y in the proof) is large. But, regardless of how small we choose ε to be, if there are regions of (S ε ) + where y is small, then in these regions V ε need not be close to V . The consequences are illustrated in the example of the Brusselator in Section 4.4, where the lifting process leads to a loss of permanence on every stoichiometric class. Note, however, that if S + , the original positive stoichiometric class of R, is itself bounded, then we can control the lifted vector field over the entirety of (S ε ) + .
Remark 5 (Normally hyperbolic invariant manifolds persist). Although Theorem 1 is phrased in terms of equilibria and periodic orbits, the result admits generalisation. Indeed, with the assumptions of the theorem, if O is any positive, compact, invariant manifold admitted by R and normally hyperbolic relative to its stoichiometric class S, then it survives C 1 perturbations [19,21], and hence is admitted by R . If, for example, R admits a k-dimensional torus on some positive stoichiometric class, and the torus is normally hyperbolic relative to this class, then the same holds for R .
Remark 6 (Bifurcations persist). Suppose R admits, on some positive stoichiometric class, a nondegenerate local bifurcation of an equilibrium or periodic orbit, unfolded nondegenerately by the rate constants. Then, for sufficiently small, fixed, ε > 0, R admits the same nondegenerate local bifurcation on some positive stoichiometric class as we vary the same combination of rate constants (note that rate constants of R and R are in natural one-to-one correspondence). Essentially, the nondegeneracy and transversality conditions associated with the bifurcation, allow us to continue the bifurcation as we vary ε. Moreover, these conditions continue to hold for sufficiently small ε. For a concrete example demonstrating the pesistence of a bifurcation, see Section 4.2.
Remark 7 (A generalisation of previous claims). For CRNs with power-law kinetics, Theorem 1 generalises the claims in Theorems 3 in [8] and [3], which treat the very special case where the added species figures only trivially in reactions, i.e., adds only a row of zeros to the stoichiometric matrix of the network. Note that, in that case, the new stoichiometric classes were bounded if and only if the original stoichiometric classes were bounded.
Remark 8 (Generalisations to other classes of kinetics). Phrasing the result in terms of powerlaw kinetics simplifies the proof, but is not key to it. The broad template of the proof can be applied to CRNs with other classes of kinetics. Remark 9 (CRNs with bounded stoichiometric classes do not have greatly restricted dynamics). One immediate consequence of Theorem 1 is that insisting a CRN has bounded stoichiometric classes, does not greatly restrict its behaviour. If a given CRN with unbounded stoichiometric classes admits some finite set of bounded nondegenerate limit sets on one of its stoichiometric classes, then we can always construct, by adding in a dependent species, a CRN with bounded stoichiometric classes which admits the same bounded nondegenerate limit sets on one of its stoichiometric classes. We see several instances of this in the examples presented in Section 4.
The next corollary tells us that it may be helpful to examine full-rank subnetworks of a CRN: finding nontrivial behaviours in these subnetworks is sufficient to ensure that they occur in the original CRN. We recall the definition of an induced subnetwork of a CRN from [3]: this is a CRN obtained by removing some reactions from a CRN, and/or some species from all the reactions in which they figure. In terms of the Petri-net graph of the CRN this corresponds to removing some vertices from the graph along with all their incident arcs.
Corollary 1. Let R be a CRN of rank r, and let R 0 be any rank-r induced subnetwork of R. If, for some fixed power-law kinetics, R 0 admits k 1 positive nondegenerate (resp., hyperbolic, resp., linearly stable) equilibria and k 2 positive nondegenerate (resp., hyperbolic, resp., linearly stable) periodic orbits, then the same holds for R with any derived power-law kinetics.
Proof. Clearly R can be built from R 0 by adding linearly dependent reactions and linearly dependent species to R 0 . The result is thus an immediate consequence of Theorem 1 above which deals with the case of adding linearly dependent species, and Theorems 1 in [8] and [3] which deal with the case of adding linearly dependent reactions. (Although Theorems 1 in [8] and [3] are stated in more restricted terms, the generalisations required are immediate.)

Examples
First, we introduce some terminology relevant to the examples below.
Homogeneous CRNs. The molecularity of a complex in a CRN is the sum of its stoichiometric coefficients. We call a CRN homogeneous if, for every reaction, the molecularities of the reactant complex and the product complex are equal. Clearly, this condition is equivalent to (1, 1, . . . , 1) t being an element of the kernel of Γ t . In particular, the stoichiometric classes of a homogeneous CRN are bounded: for the ODE associated with a homogeneous CRN, regardless of the precise nature of the kinetics, the quantity x 1 + x 2 + · · · + x n is conserved. A partial converse is also true: it is easily shown that if a CRN is endowed with any fixed power-law kinetics, and x 1 + x 2 + · · · + x n is constant along trajectories for some open set of rate constants, then the CRN is homogeneous.
Homogenisation of CRNs. Starting with an arbitrary CRN R, one can make it homogeneous by adding a new species with appropriate stoichiometric coefficients to the reactant or product complex of each reaction [14, Exercise 4 on page 29]. This operation preserves the rank of R, and the homogenisation can be carried out in multiple ways. By Theorem 1, if R has mass action kinetics, and nondegenerate multistationarity (resp., oscillation) occurs in R, then it also occurs in the homogenised CRN.
In all of the examples below, we start with a network that is not homogeneous, and then homogenise it. The examples in Sections 4.1 and 4.2 illustrate our main result, while those in Sections 4.3 and 4.4 demonstrate its limitations. The rank of each of these networks is one or two. For the homogenisation of some rank-three mass-action systems, consult [11]. For applications of Theorem 1 where the enlarged network is not homogeneous, but nevertheless has bounded stoichiometric classes, see [10].

Schlögl model: a single-species CRN with multiple equilibria
We use the reversible version of the Schlögl model [34,35] to demonstrate the use of Theorem 1 for guaranteeing the existence of multiple nondegenerate equilibria in an enlarged reaction network. It also provides some insight into the proof of Theorem 1.
Consider the following single-species mass action system and its associated differential equation: It has 3 positive nondegenerate equilibria: at x = 1, x = 2, and x = 3. The first and third are linearly stable, while the one at x = 2 is linearly unstable.
We now homogenise the network. Note that the simplest choice of homogenisation gives the enlarged network Y X, 2X + Y 3X. However, we choose the slightly more complicated homogenised network 2Y X+Y, 2X+Y 3X in order to demonstrate two points: that there are many ways to homogenise a CRN; and that the modified rate constants given in the proof of Theorem 1 can have nonlinear dependence on a perturbation parameter ε as seen below.
Theorem 1 now tells us that this network must admit three nondegenerate equilibria on some stoichiometric class, two linearly stable and one linearly unstable. To see why, we follow the proof of Theorem 1 and obtain the following mass action system and its associated differential equation, dependent on a new parameter ε: = −x 3 + 6εx 2 y − 11εxy + 6ε 2 y 2 , y = +x 3 − 6εx 2 y + 11εxy − 6ε 2 y 2 .
In the homogenised system, the stoichiometric subspace remains 1-dimensional, and the quantity x+y is conserved. We now restrict attention to the stoichiometric class defined by x+y = 1 ε and replace y by 1 ε − x. The dynamics of x for 0 < x < 1 ε is then given bẏ On any compact subinterval of (0, ∞), the vector field −x 3 +6x 2 (1−εx)−11x(1−εx)+6(1−εx) 2 converges uniformly to −x 3 + 6x 2 − 11x + 6 as ε → 0. It is not hard to see that for all sufficiently small ε the lifted system with the scaled rate constants has 3 positive equilibria in the stoichiometric class x+y = 1 ε , two of which are stable, while one is unstable. Note, however, that for any fixed rate constants, each stoichiometric class defined by x + y = c with c > 0 being large enough, has a unique positive equilibrium. This illustrates that in order to obtain the desired behaviour we must simultaneously choose the rate constants and a stoichiometric class of the lifted system.
Let us now homogenise the LVA. The resulting network and its associated mass action differential equation take the form Theorem 1 tells us that the homogenised LVA must admit a linearly stable periodic orbit on some stoichiometric class. Let us use these systems to demonstrate the arguments in the proof of Theorem 1 and some subsequent remarks. Consider again the ODE systems associated with the LVA,ẋ = κ 1 x 2 − κ 2 x 3 − κ 3 xy, and with the homogenised LVA, where we have now added a prime to the rate constants to distinguish them from those of the original LVA, If we set (κ 1 , κ 2 , κ 3 , κ 4 ) = (εκ 1 , κ 2 , κ 3 , κ 4 ), and restrict attention to the stoichiometric class on which x + y + z = 1 ε , then, using this equation to eliminate z, we find that x and y in the homogenised LVA evolve according tȯ We see immediately that the vector field in (6) converges on compact sets to that in (4) as ε → 0.
In particular, suppose we fix some values of κ 1 , κ 2 , κ 3 , κ 4 such that κ 1 κ 3 − 2κ 2 κ 4 is positive and sufficiently small to ensure that (4) has a linearly stable, positive periodic orbit O. Then, for sufficiently small ε > 0, (6) has a linearly stable, positive, periodic orbit O ε . Moreover, provided ε is sufficiently small, z = 1 ε − x − y remains positive as x and y vary along this periodic orbit, and so the lifted system (5) has a positive periodic orbit O ε on the stoichiometric class defined by x + y + z = 1 ε . Clearly, since O ε is linearly stable for (6), O ε is linearly stable relative to this class.

Lotka reactions
We now provide an example where the assumptions of Theorem 1 are violated, and periodic solutions are not preserved by lifting.
Consider the Lotka reactions and the associated mass action differential equation: In this case, the unique positive equilibrium κ3 κ2 , κ1 κ2 is surrounded by a continuum of periodic orbits; these are level sets of the nonlinear first integral x −κ3 y −κ1 e κ2(x+y) . Since these periodic orbits are degenerate, Theorem 1 does not apply.
Indeed, adding a new species, Z, to some of the reactions while preserving the rank of the network can lead to the destruction of all of the periodic orbits. Consider the following network and its associated mass action differential equation: The set of positive equilibria is κ3 κ2 , κ1 κ2 t, t : t > 0 . Thus, the stoichiometric classes x+y+z = c with c ≤ κ3 κ2 have no positive equilibria, while those with c > κ3 κ2 have a unique positive equilibrium. Note that the divergence of the vector field after division by xyz equals − κ3 xz 2 . Since this quantity is negative on R 3 + , there is no periodic orbit that lies entirely in the positive orthant [36, Satz 1] (see also Remark (v) following the proof of Theorem 2.3 in [26]). In fact, every positive equilibrium is globally asymptotically stable within its positive stoichiometric class: this follows via the Poincaré-Bendixson Theorem, since it can be shown that on any stoichiometric class which includes a positive equilibrium no positive initial condition has omega limit set intersecting the boundary of the nonnegative orthant.

Brusselator
Our final example demonstrates that while the lifted CRN must admit the nondegenerate behaviours of the original CRN, it may also allow other behaviours not seen in the original CRN, such as multiple periodic orbits and homoclinic orbits. Looked at from another angle, omitting a single species from a CRN, even without changing its rank, can result in the loss of many different nontrivial behaviours.
Consider the Brusselator and its associated mass action differential equation: At any fixed values of the rate constants the system has a unique positive equilibrium, (x * , y * ) = κ1 κ2 , κ2κ3 κ1κ4 . It can be shown that the system is permanent. Moreover, it is known that (x * , y * ) is globally asymptotically stable for κ 3 ≤ κ 2 + , while it is repelling for κ 3 > κ 2 + κ 2 1 κ4 κ 2 2 and is surrounded by a periodic orbit that is born via a supercritical Andronov-Hopf bifurcation. Moreover, this periodic orbit is unique and attracts every positive initial condition, except (x * , y * ) (see [42, Example 5 on page 135]).
Let us now homogenise the Brusselator. The resulting CRN and its mass action differential equation take the form By Theorem 1, there exist rate constants and a stoichiometric class where the new system has a stable periodic orbit. In fact, given any rate constants satisfying κ 3 > κ 2 + κ 2 1 κ4 κ 2 2 and such that the Brusselator has a linearly stable periodic orbit, the proof of Theorem 1 gives us a procedure for finding a linearly stable periodic orbit in the lifted system. Moreover, by Remark 6, a supercritical Andronov-Hopf bifurcation must occur on some stoichiometric class as rate constants are varied in the lifted system. Note, however, that the global behaviour of the homogenised Brusselator is quite different from that of the Brusselator. The first interesting difference is that whereas the original system was permanent, the lifted system is not permanent on any stoichiometric class intersecting the positive orthant; we show in Appendix A that for all rate constants and all positive values of the parameter c, the stoichiometric class satisfying x + y + z = c includes a boundary equilibrium (0, c, 0) which is asymptotically stable relative to this class. Thus we cannot control the vector field at all points on any of these stoichiometric classes, even though these classes are compact (see Remark 4 following the proof of Theorem 1).
Aside from the loss of permanence, the homogeneous Brusselator admits several nontrivial behaviours forbidden in the original. In Appendix A we show that for the homogenised network both supercritical and subcritical Andronov-Hopf bifurcations can occur on some stoichiometric classes as we vary rate constants. Furthermore, various interesting co-dimension two bifurcations take place: a generic Bautin bifurcation [25,Section 8.3] and a generic Bogdanov-Takens bifurcation [25,Section 8.4] can both occur. This implies that the lifted system admits, for various choices of the rate constants, the following behaviours not seen in the original Brusselator: • a fold bifurcation of equilibria; • an unstable positive equilibrium surrounded by a stable and an unstable periodic orbit; • a fold bifurcation of periodic orbits; • a stable positive equilibrium surrounded by an unstable periodic orbit or a homoclinic orbit; • a homoclinic bifurcation.
It is not unexpected that the lifted system admits richer dynamical behaviour; but it is worth noting that if we restrict attention to any stoichiometric class and consider the evolution of concentrations of X and Y, then the (2-dimensional) lifted vector field differs from the original only by the addition of linear terms (see equation (A.1) in Appendix A). The addition of linear terms to a 2D differential equation can thus quite dramatically increase the complexity of its behaviour.

Discussion and conclusions
Inheritance results tell us how we might enlarge a CRN while preserving its capacity for various dynamical behaviours. The main result in this paper adds a simple but important inheritance result relevant to the study of both multistationarity and oscillation in CRNs. Some inheritance results, including the main result of this paper, are gathered in [5]. Taken together, these results provide a powerful tool for predicting nontrivial behaviours in a CRN based on its subnetworks.
It is useful to think of inheritance results in terms of partial orders. For example, suppose a CRN R with mass action kinetics admits linearly stable oscillation. The same then holds for all CRNs greater than R in the partial order defined by available inheritance results. We may then look for CRNs admitting linearly stable oscillation which are minimal with respect to this partial order in order to gain insight into the capacity for stable oscillation in larger and more complex CRNs [3]. Such a program provides a rigorous basis for claims about "motifs", namely small subnetworks which are at the root of certain behaviours in biological systems. Systematically identifying minimal CRNs with prescribed behaviours, followed by the development of algorithms to test for their presence in larger CRNs, is a natural avenue for future work.
While inheritance results are often phrased in terms of enlarging CRNs, they can also be seen in terms of modelling choices. The main result here is relevant to choices which might affect conservation laws. In physically realistic systems of chemical reactions occurring in a closed environment we expect numbers of atoms of each element to be conserved: this is the so-called law of atomic balance [14]. However, it is common in modelling CRNs to omit some species from reactions, particularly when they are considered to be present in abundance, or their concentration is subject to external control. As an example, when reactions involving ATP and ADP occur in biochemical models, inorganic phosphate and water are often omitted from the equations. Such omissions tend to destroy physical conservation laws.
It is natural to worry that omitting species from reactions might introduce fundamentally new, and unrealistic, behaviours into the system. Theorems 4 in [8] and [3] provide some reassurance that this will not occur if we omit species whose concentration is sufficiently strongly controlled by external processes. In a similar way, Theorem 1 here provides some reassurance in the case where the omitted species are linearly dependent on the others. In this case, the omissions cannot introduce the capacity for behaviours such as stable oscillation or multistability. In fact, any behaviour of the simplified CRN occurring on a compact set, and which is robust in the sense that it survives C 1 perturbation, can also be obtained on some stoichiometric class of the larger CRN for appropriate choices of parameters. We can, of course, still lose interesting behaviours by omitting dependent species, as illustrated by the example of the Brusselator in Consider the homogenised Brusselator and its associated mass action differential equation: The stoichiometric classes which intersect the positive orthant are P c = {(x, y, z) ∈ R 3 ≥0 : x + y + z = c} for c > 0. Note that for each positive c, the corner (0, c, 0) of the triangle P c is an equilibrium of this system (and this is the only boundary equilibrium in P c ).

Appendix A.1. The homogenised Brusselator is not permanent
Let us now restrict the dynamics to P c . After elimination of z by the conservation law x+y+z = c, we obtain the following ODE system on {(x, y) ∈ R 2 ≥0 : x + y ≤ c}: From here on, we will focus on this system, which can be seen as describing, in local coordinates, the dynamics of the original ODEs restricted to a particular stoichiometric class, parameterised by c. Note that parameters of the system are now the four original rate constants κ 1 , κ 2 , κ 3 , κ 4 , along with c. For the purposes of bifurcation analysis we assume, however, that c is fixed.
Observe that the Jacobian matrix of (A.1) at the corner equilibrium (0, c), denoted by J 0 , equals Since det J 0 > 0, tr J 0 < 0, and both eigenvalues are real and negative. Therefore, the boundary equilibrium (0, c) is asymptotically stable. In particular, for all (κ 1 , κ 2 , κ 3 , κ 4 , c) ∈ R 5 + , (A.1) is not permanent, and since c > 0 was arbitrary, the same clearly holds for the homogenised Brusselator on each positive stoichiometric class. By contrast, the Brusselator was permanent for all positive values of the rate constants.

Appendix A.2. Local bifurcation analysis
We proceed as follows. We first parameterise the set of equilibria which simplifies many calculations. The parameter used, denoted by t, can replace the parameter c in many calculations: along the branch of equilibria we consider, t and c are in one-to-one correpondence. We write down necessary conditions for fold and Andronov-Hopf bifurcations to occur, and confirm that both supercritical and subcritical Andronov-Hopf bifurcations can occur. We also confirm that two codimension-2 bifurcations can occur on stoichiometric classes as we vary the rate constants: a generalised Andronov-Hopf bifurcation, also known as a Bautin bifurcation; and a Bogdanov-Takens bifurcation. We can check that apart from along exceptional sets, nondegeneracy and transversality conditions hold for all the bifurcations. See [9] for full details of all the calculations.
In the positive stoichiometric classes P c with 2 positive equilibria, the equilibrium corresponding to the value of t < t * , which is closer to the corner equilibrium (0, c), is a saddle for all values of the rate constants. For the purposes of local bifurcation analysis, we thus focus attention on equilibria satisfying t ≥ t * . As values of t and c are in one-to-one correspondence for t ≥ t * and c ≥ c * , we can thus pass back and forth between parameters c and t, and do so, sometimes without comment, in calculations to follow.
Fold bifurcations of equilibria. These potentially occur when det J(t) = 0, namely on the set T in parameter space where t = t * : Equivalently, fold bifurcations potentially occur when c = c * , namely along: For any fixed value of c, T tells us combinations of rate constants for which we expect a fold bifurcation to occur on the stoichiometric class parameterised by c. It is straightforward to confirm that, provided κ 3 = (κ 1 +κ 2 ) 2 /κ 2 , the fold bifurcation is nondegenerate, and is unfolded nondegenerately by the rate constants κ i (any of κ 1 , κ 2 , κ 3 or κ 4 serves to unfold the bifurcation nondegenerately for all rate constants). The degenerate case where κ 3 = (κ 1 + κ 2 ) 2 /κ 2 will be of importance later.
Andronov-Hopf bifurcations. These potentially occur when tr J(t) = 0 and det J(t) > 0, namely, when Combining these conditions, Andronov-Hopf bifurcations potentially occur along  . The black curve is where ab = (1 + a) 2 , namely the boundary of H where the determinant of the Jacobian matrix vanishes and we expect Bogdanov-Takens bifurcations to occur; notice that it touches only the L 1 > 0 region, and we can indeed confirm that all Bogdanov-Takens bifurcations are subcritical. Along the red curve, where L 1 = 0 we can confirm that L 2 > 0. This is where we expect Bautin bifurcations to occur and, in fact, we find these bifurcations to be nondegenerate away from an exceptional set.
Note that the necessary condition for Andronov-Hopf bifurcation κ 2 κ 3 > (κ 1 +κ 2 ) 2 is equivalent to ab > (1 + a) 2 . Since Q(a, b) is positive whenever this condition is satisfied, it suffices to investigate the sign of P (a, b). As a and b vary along H, P (a, b), and hence L 1 , can be positive, negative or zero along the bifurcation set, corresponding to subcritical, supercritical, and degenerate Andronov-Hopf bifurcations respectively, as shown in Figure  We can also confirm that wherever L 1 = 0, the parameters unfold the bifurcation nondegenerately.
Indeed, the nondegeneracy conditions for the Bautin bifurcation can be checked and found to hold on GH apart from along an exceptional set: 1. Along GH, we can compute the second focal value, L 2 , and find it is always positive. The calculations are lengthy and are omitted, but can be found in the supporting documentation [9]. 2. The parameters (κ 3 , κ 4 ) unfold the bifurcation nondegenerately. Choosing these as our bifurcation parameters, then for fixed choices of the remaining parameters, we can confirm Thus, away from the exceptional set, the Bautin bifurcation is nondegenerate. For parameter values near to GH, on certain stoichiometric classes an unstable positive equilibrium is surrounded by a pair of periodic orbits, one stable and one unstable. As we vary rate constants, a fold bifurcation of periodic orbits must occur on some stoichiometric class.
We can confirm that the B-T bifurcation is always nondegerate, subcritical, and unfolded nondegenerately by parameters (κ 3 , κ 4 ). This involves checking conditions (BT.0), (BT.1), (BT.2), and (BT.3) in [25,Theorem 8.4]. The computations can be carried out explicitly, but are lengthy, and so are omitted here. They can be found in the supporting documentation [9]. Thus for parameter values close to BT, on some stoichiometric classes, the system can have a stable equilibrium surrounded by an unstable periodic orbit. Furthermore, homoclinic bifurcations can occur.
In Figure A.2, we fix κ 1 = 2, κ 2 = 4, and c = 6, and depict the curves corresponding to T and H, and the points corresponding to BT and GH, in the (κ 3 , κ 4 ) plane. This completes our analysis of the local bifurcations of equilibria in the homogeneous Brusselator.