Conservation Laws in Biochemical Reaction Networks

We study the existence of linear and non-linear conservation laws in biochemical reaction networks with mass-action kinetics. It is straightforward to identify the linear conservation laws as they are related to the left null-space of the stoichiometry matrix. The non-linear conservation laws are much more diﬃcult to study and so far have rarely been considered in the context of mass-action reaction networks. Our aim is to give structural conditions – that is, parameter independent conditions, on a reaction network to guarantee the existence of non-linear conservation laws of some type. We do so by means of Darboux theory of integrability. We show that F ( x ) = x i is a Darboux polynomial if and only if the reaction network fulﬁl a certain structural condition. Furthermore, this allows us to conclude that a speciﬁc type of a non-linear ﬁrst integral (similar to that of the Lotka-Volterra system) only exists if the reaction network fulﬁls the same structural condition. The existence of such a ﬁrst integral generally implies that the system is persistent and has stable steady states. We illustrate our results by examples.


Introduction
Typical models of biochemical reaction networks are systems of polynomial ordinary differential equations (ODEs).Generally, these ODE systems are difficult to solve and analyse, and little can be said about their solutions, exactly or qualitatively, except in special cases [8,24].However, in many applications it is important to know whether the solutions oscillate, are attracted towards steady states, or are persistent.
To study the dynamics, a potentially promising strategy would be to look for first integrals (i.e.conservation laws).First integrals are quantities depending on the variables of the system, that are conserved over time.Typically, differential systems (equivalently vector fields) do not admit first integrals, however, when they do, and can be found, the benefit can be significant [24,15,3].In connection with reaction networks, first integrals have been used to demonstrate oscillatory behaviour [24], chaos [23] and Turing instability [22] in particular examples.
The ODE system of a reaction network is composed of two parts: a stoichiometric (structural) representation of the reactions and a specification of the rates by which the reactions occur.Only the second part is parameter dependent.Our aim is twofold.First of all, we are interested in finding sufficient and/or necessary structural conditions -that is, conditions that are independent of the parameters -for the existence of non-linear first integrals.Secondly, for a given reaction network, we ask, what can we learn about the qualitative behaviour of the solutions from the existence of non-linear first integrals?We make use of Darboux theory of integrability [4,5], which is a powerful tool in determining the non-linear first integrals [18,16].
In relation to the first aim, we show that F (x) = x i is a Darboux polynomial (see Definition 11) if and only if the reaction network fulfils a certain structural condition on the form of the reactions.This in turn allows us to conclude that a particular type of non-linear first integral only exists for a class of reaction networks defined by the structural condition.This class contains a standard reaction network-representation of the Lotka-Volterra model.We find this re-sult appealing because the condition is easy to check, simply by looking at the reactions.However, the condition is only necessary and not sufficient for the existence of a first integral of the desired form.We show this by example and demonstrate the first integral might not exist at all, or might exist for some or all parameter values.For the second aim, we show that the existence of a first integral of the desired form for a conservative system implies that it is persistent and has a unique non-attracting, but stable, steady state in each stoichiometric compatibility class (these classes are determined by linear first integrals).
The use of linear first integrals in reaction network theory is common.With some exceptions the non-linear first integrals have rarely been considered [1,17,20].Quadratic first integrals have been characterised in [20].Ideally, one would in general like to give conditions for the existence of first integrals of certain forms, irrespectively of the reaction rate constants.Our results might be seen as a first step in this direction.Also our results are in this sense analogous to other results in reaction network theory, for example, the deficiency zero and one theorems that give structural conditions for the existence of steady states [8], and the persistence criteria given in [2].
It is generally difficult to establish the existence or non-existence of Darboux first integrals.This typically depends on the reaction rate constants in intricate ways [12,13,17,11].We hope that this paper would inspire further work to characterise the existence and form of non-linear first integrals for reaction networks.As our results show, non-linear first integrals might be useful to study qualitative aspects of the dynamics.

Reaction networks with mass-action kinetics
Here we define a reaction network and show how to construct a system of polynomial ODEs describing the evolution of the species concentrations.

Reaction networks
Denote by N the set of non-negative integers, and by R >0 the set of positive (non-negative) real numbers.For In the following, Ω denote an open set of R n .
(b) C has p elements, called complexes, being a linear combination of species, (c) R has m elements, called reactions, and each reaction r ∈ R is an ordered pair of distinct complexes, written y → y , and therefore R ⊂ C × C.
In the following we assume that the sets S, C, and R are ordered.It is common to assume that C only contains complexes which are part of reactions, and S only contains species which are part of complexes.In that case, S and C are determined by R.

Stoichiometric subspace and compatibility class
Consider a reaction network (S, C, R).Identify each species A i with the i-th unit vector in N n , that is, the vector with one in the i-th entry and zero elsewhere.Thus, a complex y = α 1 A 1 +. ..+α n A n is identified with the element For the j-th reaction y j → y j , we let γ j = y j −y j ∈ R m be the net production of species in the reaction.The vector γ j is referred to as the j-th reaction vector.
Definition 4 (Stoichiometric compatibility class).The stoichiometric compatibility classes are the affine vector subspaces in R n defined by (x + S) ∩ R n ≥0 , where x ∈ R n ≥0 .

Reaction networks with mass-action kinetics
Let (S, C, R) be a reaction network.Denote by x i the concentration of the species A i and the vector of species concentrations by x(t) = (x 1 (t), . . ., x n (t)).
Before we show how to describe the time evolution of the vector of concentrations, we need some preliminary definitions and notation.
A rate function for a reaction y j → y j , j = 1, . . ., m, is a function v j (x) : R n ≥0 → R ≥0 that describes the instantaneous change in the species composition x = (x 1 , . . ., x n ) due to this reaction.
A kinetics V is an assignment to each reaction y j → y j , j = 1, . . ., m, of a rate function v j (x) : R n ≥0 → R ≥0 .We can think of V as a set of rate functions indexed by the elements of the reaction set, that is V = {v j (x) : j = 1, . . ., m}.
Definition 5 (ODE system of a reaction network).Let (S, C, R) be a reaction network.The evolution of the concentrations x(t) = (x 1 (t), . . ., x n (t)) under the kinetics V is given by the ODE system where v(x) = (v 1 (x), . . ., v m (x)) is the vector of the rate functions.
It follows from ( 2) that the evolution of the species concentrations is confined to the stoichiometric subspace Non-negativity of the solutions follows from the assumption that the rate functions are non-negative, in particular at the boundary of R n ≥0 [21].
We are particularly interested in mass-action kinetics.
Definition 6 (Mass-action kinetics).A kinetics is called mass-action if the rate functions take the form where y j = (α j1 , . . ., α jn ), and k j > 0 is a positive reaction rate constant.We denote the vector of reaction rate constants by κ = (k 1 , . . ., k m ).
Example 2 (Volpert's network, part 2).Consider the reaction network given in (1) with the mass-action kinetics.The corresponding system of ODEs is given by where

Linear conservation laws
Proposition 7. Let (2) be an ODE system for a reaction network.For any row vector ω = (ω 1 ω 2 . . .ω m ) such that ω Γ = 0, the quantity is a linear conservation law.Moreover, there are at least s = n − rank(Γ) independent linear conservation laws.
Proof.We have The left kernel of Γ has dimension s, hence there must be at least s independent linear conservation laws.
Remark 8.A linear conservation law is a linear first integral, see Definition 10 below.

Consider the reaction network
In this case, s = 1 as is not in the kernel of Γ.The two vectors ω 1 and ω 2 are linearly independent for all reaction rate constants.
)) ⊆ S for a reaction network with mass-action kinetics, as in (3).If the number of linkage classes equals the number of terminal linkage classes, then S κ = S and the number of linear conservation laws is s = n − rank(Γ).
In (5), there is one linkage class but two terminal linkage classes (A 1 , and A 2 ), hence the proposition does not apply.As a second example, consider the reaction network The function in the left kernel of Γ, that is, ω is not of the form in Proposition 7. Indeed Example 2 (see also Example 1) has three linkage classes and three terminal linkage classes, and hence there is only s = 1 conservation law for all κ, namely

Non-linear conservation laws via Darboux theory
The theory developed by Darboux [4,5] is one of the most useful tools for identifying the first integrals of polynomial differential equations.Here our aim is to review the main elements of the Darboux method and subsequently show how it can be applied to search for nontrivial (non-linear) conservation laws of the reaction networks with mass-action kinetics.
We start by considering an n-dimensional polynomial differential system where x = (x 1 , . . ., x n ) ∈ R n and P i (x) ∈ R[x] are polynomials in x.Whenever convenient we omit the argument x in P i = P i (x) and other functions of x.
For any system of differential equations (7) we define a derivation by Let In the examples, Ω can be taken to be R n >0 .
The degree of the system (7), or of the derivation (8), is defined as where deg(P ) denotes the total degree of a polynomial P .
Definition 10 (First integral).We say that a C 1 -function H : Ω → R, where Ω ⊂ R n is a first integral of system (7), or of the derivation D, if D(H) = 0, and H is not locally constant on any positive Lebesgue measure subset of Ω.
Definition 11 (Darboux polynomial).A Darboux polynomial of the system ( 7) for some cofactor Remark 12.A polynomial first integral is also a Darboux polynomial with cofactor K = 0.If the degree of the derivation (8) is d, then the degree of any cofactor in ( 9) is bounded by d − 1.
Definition 13 (Exponential factor).An exponential factor of the system (7) is where the degree of L is lower than the degree of D. Since the number of linearly independent polynomials in two variables of degree at most d − 1 is d+1 2 , there must be a linear combination of cofactors summing to zero, whenever r + s is strictly bigger than d+1 2 .Hence Theorem 15 applies.
Theorem 16.Assume the derivation D admits r Darboux polynomials F i with respective cofactors K i , i = 1, . . ., r, and s exponential factors E j with respective cofactors L j , j = 1, . . ., s.If there exist λ 1 , . . ., λ r , µ 1 , . . ., µ s in C, such that r i=1 then the Darboux function is a first integral of (7), provided Φ is not locally constant on any positive Lebesgue measure subset of its domain of definition.
Proof.We first note that if Φ(x) is a first integral, then the composition is another first integral, where F is any C 1 -function.Now we prove that log Φ is a first integral.We have The theorem follows from (10).

A particular class of reaction networks
We start by making a connection between the existence of certain Darboux polynomials and the structural form of a class of reaction networks.This form is shared by many reaction networks, in particular the standard representation of the Lotka-Volterra system.
Proof.The equivalence between (i) and (ii) follows from the definition of a Darboux polynomial.
Assume (iii).Consequently x i is a factor of (α ji − α ji )x yj , if the term is not zero.If α ji = 0 and α ji > 0, then x i is also a factor of (α ji − α ji )x yj = −α ji x yj by definition of mass-action kinetics.Hence P i = m j=1 (α ji − α ji )x yj factorises as P i = Kx i for some polynomial K ∈ C[x].If P i is zero, then K = 0 and (ii) is proved.
For the reverse statement, assume (ii), that is, P i = Kx i (which might be zero), for some K ∈ C[x], and α ji > 0 for some j = 1, . . ., m.We need to show that α ji > 0. Assume the opposite that α ji = 0. Then α ji k j x yj is a term in P i , not involving x i , by definition of mass-action kinetics.Hence it must cancel with terms from other reactions with the same monomial x yj , that is, 0 = {j : y j =yj } (α j i − α j i )k j x yj .Since x i is not a factor of x yj , we have α j i = 0, and hence 0 = {j : y j =yj } α j i k j x yj ≥ α ji k j x yj , implying α ji = 0, and we have reached a contradiction.Condition (iii) is independent of the parameters of the reaction network, and hence it is fulfilled for all κ or for none at all.Thus, we have given a structural characterisation of networks with Darboux polynomials F = x i for all i.System (6) fulfills all three statements of Proposition 17 with K = k 1 − k 2 .
Also system (5) fulfills all three statements of Proposition 17 for i = 1, that is,  4), it follows from Proposition 17, that the system admits the Darboux polynomials F i (x) = x i for i = 1, 2, 3 with cofactors Now we aim to check if the above cofactors are linearly dependent.Since this is so, and the network admits, according to Theorem 15, the Darboux first integral for any values of the reaction rate constants k 1 , k 2 , k 3 > 0. In addition, according to Proposition 7 and 9, there is one independent linear conservation law H 0 = Theorem 18.Let (S, C, R) be a reaction network with mass-action kinetics and reaction rate constants κ and let I ⊆ {1, . . ., n}.If is a first integral of the reaction network, then (ii) if α ji > 0 then α ji > 0 for all j = 1, . . ., m and i ∈ I.
Proof.If (i) is fulfilled, then x i is a Darboux polynomial of P i [6].Then (ii) follows from Proposition 17.
The proof shows that under the assumptions of the theorem, there is a linear combination of cofactors for the Darboux polynomials F i = x i and the exponential factors G i = e δixi , yielding the first integral H.The reverse statement is not true in general: system (6) fulfills condition (ii) in Theorem 18, but cannot have a conservation law of the form (i) as it is one-dimensional (unless in which case x(t) = x(0) for all t > 0).

Example 4 (Lotka-Volterra network). Consider the following version of the
(prey reproduce, predators eat prey and reproduce, predators die).Proposition 17(iii) is fulfilled, hence F 1 = x 1 and F 2 = x 2 are Darboux polynomials.The corresponding cofactors are In addition, E = exp(x 1 + x 2 ) is an exponential factor with cofactor L = is a first integral of the reaction network, which is well known [24].
Example 5 (Reversible network).Consider the following reaction network with mass-action kinetics The reaction network is said to be reversible because the second reaction is the first reversed.Proposition 17(iii) is fulfilled, and Darboux polynomials with cofactors but they are not linearly dependent and there is no first integral as in Theorem 18.In fact, since H = x 1 + x 2 is a linear conservation law, we might reduce the system to a one dimensional ODE, which cannot have a law of the form as in Theorem 18 (unless trajectories are constant).
It follows from the so-called deficiency zero theorem that the network in Example 5 has precisely one steady state in each stoichiometric compatibility class (given by H = x 1 + x 2 ) and that this steady state is asymptotically stable.
It also follows that there cannot be periodic orbits [8], which is also consequence of the fact that the system essentially is one dimensional.

First integrals and dynamics
In this section, we show by examples how to apply the Darboux method to determine non-linear conservation laws of mass-action reaction networks.Additionally, we demonstrate how first integrals might be used to study dynamical properties of the system.For that purpose we use different generalisations of We start with a definition of persistence and a lemma.
Definition 19 (Persistance).An ODE system for a reaction network is said to be persistent if any trajectory starting from a positive initial point is bounded away from the boundary of R n ≥0 , that is, Lemma 20.Let (S, C, R) be a reaction network with mass-action kinetics and reaction rate constants κ.Assume the system admits the following first integrals , and there are no other linearly independent linear first integrals.Then the system is persistent.
In addition, if δ i = 0 for all i, then there exists at least one stable positive steady state x * = (x * 1 , . . ., x * n ) in each stoichiometric compatibility class (determined by H k , k = 1, . . ., r), given by No trajectories apart from the constant trajectory x(t) = x * , t ≥ 0, are attracted towards x * , that is, there exists > 0, depending on the initial condition x(0), By assumption H k ≥ 0 for all k = 1, . . ., r.It follows that each concen- where and from H 0 that x i cannot get arbitrary close to zero as this would imply that at least one other concentration, say x j , would become arbitrarily large, contradicting that x j ≤ M for all j.Hence the system is persistent.
Choose an element i k ∈ I k for each k and define J k = I r \ {i k }.By inserting which attains its global maximum at x * = (x * 1 , . . ., x * n ), where Any trajectory starting out at x * must remain there, hence x * is a positive steady state in the given stoichiometric class.Since H 0 is continuous as a function of x, it follows that any trajectory that starts close to x * remains close to x * for all t > 0. Hence the steady state x * is stable.
For the last part, consider the right hand side of (11) as a continuous function of (x i ) i∈J , J = ∪ r k=1 J k (with fixed values of H 1 , . . ., H k ).Assume some trajectory x(t) converges towards x * as t → ∞.Hence by continuity, H 0 (x(t)) → H 0 (x * ) as t → ∞.Since H 0 is constant along trajectories, it follows that H 0 (x(t)) = H 0 (x * ), unless x(t) = x * .Remark 21.In [2] (see also [19]), a structural condition for persistence is introduced.A siphon Σ ⊆ S is a subset of the species set with the property that if A i ∈ Σ and α ij > 0 for some reaction R j , then there is a species A k ∈ Σ such that α kj > 0 for the same reaction.If all siphons of the network contain the species of a linear first integral, then the network is persistent.If a network fulfills the structural condition given in Theorem 18, then it does not satisfy this condition.Any species A i for which α ij > 0 constitutes a siphon, but Σ = {A i } will only contain the support of a linear first integral if the concentration x i of A i is constant through time.

Generalized Volpert network, example A
Consider the following chemical reaction network where i = 1, . . ., .This is a generalization of the system in Example 1 for which = 1.This general system has n = 2 + 1 species S = {A 1 , . . ., A 2 +1 }, There are m = 3 reactions.Proposition 22.For any values of the reaction rate constants k 1 , . . ., k 3 > 0 the polynomial differential system corresponding to the reaction network in (12) admits the following two first integrals of Darboux type x i , the latter being the only linear first integral of the system.The first integral is the only possible first integral of the form H 0 = n i=1 x λi i .
Proof.The differential equations associated with (12) are given by for i = 1, . . ., .In total we have 2 + 1 equations and 3 reaction rate constants.Since there are 2 + 1 linkage classes, 2 + 1 terminal linkage classes, and s = n − rank(Γ) = 1, then H 0 is the only linear conservation law according to Proposition 9, for all values of the reaction rate constants.
It follows from Proposition 17 that the system admits at least 2 +1 Darboux polynomials, namely F i = x i , where i = 1, . . ., 2 + 1, with the following 2 + 1 cofactors In order to find a Darboux first integral one needs to find λ 1 , . . ., λ 2 +1 such that This expression is zero only if each coefficient of x i is zero.It follows that any solution is proportional to Hence, according to Theorem 16, i is a first integral for some λ i ∈ C, then so is log(H 0 ).It follows that D(H 0 ) is a linear combination of the cofactors K i , hence λ i must take the form stated above.
The assumption of Lemma 20 is fulfilled for the reaction network above; 290 hence the conclusions apply.In particular the system is persistent.Volpert [24] showed that for = 1 (n = 3) any trajectory starting in a positive initial point is periodic.

Consider the following chemical reaction network
where i = 1, . . ., n.This is a generalization of the system in Example 1, for 295 which n = 3.The general system has n species, and p = 2n complexes, C = . . ., 2A n }.There are m = n reactions.
Proposition 23.For any values of the reaction rate constants κ, the polynomial differential system corresponding to the network (13) admits only one linear conservation law, namely Moreover, if n ≥ 3 is odd, then for any values of the reaction rate constants κ, there is an additional non-linear conservation law If n ≥ 4 is even and the following constraint on the rate constants is fulfilled, then there are two non-linear conservation laws of the form H 0,even = n i=1 x λi i with and for i = 1, . . ., n 2 .For general n ≥ 3, any first integral of the form H 0 = n i=1 x λi i is a product H 0 = (H 1 0,even ) γ H 2 0,even , where γ ∈ C, H 1 0,even has λ i chosen as in (15), and 300 H 2 0,even has λ i chosen as in (16).
Oppositely, if H 0 = n i=1 x λi i is a first integral for some λ i ∈ C, then so is log(H 0 ).It follows that D(H 0 ) is a linear combination of the cofactors K i , hence λ i must take the form stated above.
Lemma 20 is applicable for n odd.In order to apply Lemma 20 for n even, we need to choose two first integrals for which different λ i 's are zero and combine them into one first integral.Hence, the assumption of Lemma 20 is fulfilled and therefore also the conclusions hold.
and three species, S = {A 1 , A 2 , A 3 }.A reaction network (S, C, R) has a natural representation in terms of a (directed) graph, as in Example 1, where the vertices are the complexes and the (directed) edges between vertices are the reactions.The connected components of the undirected graph are called linkage classes and the strongly connected components of the directed graph are called terminal linkage classes.In Example 1 there are three linkage classes (e.g., A 1 + A 2 → 2A 2 ) and three terminal linkage classes (e.g., 2A 2 ).

Remark 14 .Theorem 15 . 1 n.
Notice that the polynomial F in the expression of E = exp(G/F ) is a Darboux polynomial.Moreover it can be easily deduced that D(G) = K F G + LF , where K F is the cofactor of F .The Darboux theory of integrability relates the number of Darboux polynomials and exponential factors with the existence of a Darboux first integral.A Darboux function is a product of complex powers of Darboux polynomials and exponential factors.Note that a rational function (a ratio of two polynomials) is a special case of a Darboux function.A Darboux first integral is a Darboux function that is a first integral according to Definition 10.Assume that a derivation D of degree d admits r Darboux polynomials F i , i = 1, . . ., r, and s exponential factors E j , j = 1, . . ., s.Let N = n+d−Then the following statements hold.(a) If r + s ≥ N + 1, then the derivation D admits a Darboux first integral.(b) If r + s ≥ N + n, then the derivation D admits a rational fist integral.Theorem 15(a) is due to Darboux [5, 4], whereas Theorem 15(b) is attributed to Jouanolou [14].Theorem 15 casts the burden of finding first integrals into finding Darboux polynomials.As an illustration of the bound in Theorem 15(a) consider the case n = 2, that is, a system in two variables.As noted earlier, for Darboux polynomials F i with cofactors K i and exponential factors E j with cofactors L j , we have max{deg K 1 , . . ., deg K r , deg L 1 , . . ., deg L s } ≤ d − 1.

Example 3 (
Volpert's network, part 3).Consider Volpert's reaction network in Example 1 and 2 with mass-action kinetics.According to Theorem 15 guarantee the existence of a conservation law.However, fewer might suffice if some linear combination of cofactors equals zero, see Theorem 15.Using (1) and (

Volpert's reaction
network (Example 1) that fall into the class of reaction networks characterised in Proposition 17.