# A Computational Approach to Complete Exact and Approximate Conservation Laws
of Chemical Reaction Networks

Aurélien Desoeuvres1, Alexandru Iosif7, Christoph Lüders6,
Ovidiu Radulescu1,∗, Hamid Rahkooy4, Matthias Seiß2, Thomas Sturm3,4,5

( 1 University of Montpellier and CNRS UMR5235 LPHI, Montpellier, France
2 University of Kassel, Kassel, Germany
3 CNRS, INRIA, and the University of Lorraine, Nancy, France
4 Max Planck Institute for Informatics, Saarbrücken, Germany
5 Saarland University, Saarbrücken, Germany
6 University of Bonn, Bonn, Germany
7 Saint Louis University Madrid Campus, Madrid, Spain
∗ corresponding author<EMAIL_ADDRESS>
)

###### Abstract

We consider chemical reaction networks (CRN) with monomial reaction rates. For
these models, we introduce and discuss the concepts of exact and approximate
conservation laws, that are first integrals of the full and truncated sets of
ODEs. For fast-slow CRNs, truncated ODEs describe the fast dynamics. We define
compatibility classes as subsets of the species concentration space, obtained
by equating the conservation laws to constants. A set of conservation laws is
complete when the corresponding compatibility classes contain a finite number
of steady states (quasi-steady states for approximate conservation laws of
fast-slow systems). Complete sets of conservation laws can be used for
reducing CRNs by variable pooling. We provide algorithmic methods for
computing linear, monomial and polynomial conservation laws of CRNs and for
testing their completeness. The resulting conservation laws and their
completeness are are either independent or dependent on the parameters. In the
latter case, we provide parametric case distinctions. In particular, we
propose a new method to compute polynomial conservation laws by comprehensive
Gröbner systems and syzygies.

Keywords: First integrals, chemical reaction networks, polynomial conservation
laws, syzygies, comprehensive Gröbner systems.

##  1 Introduction

Chemical reaction networks (CRN) have important applications in physics and
chemistry, where they are used as models for non-equilibrium thermodynamics
[54], catalytic reactions [72], combustion [60], etc. More recently, CRNs were
used to model cell and tissue physiology [67] needed for the understanding of
the fundamental mechanisms of living systems and for fighting disease.

A simple characterisation of CRNs uses stoichiometries and reaction rates [2].
The stoichiometry of a chemical reaction tells us how many molecules of each
species are consumed and produced in the reaction. For instance, the reaction
$A_{1}+A_{2}\rightarrow A_{3}$ consumes one molecule of each $A_{1}$ and
$A_{2}$ and produces one molecule of $A_{3}$. The reaction rate is a function
of species concentrations, representing the number of occurrences of the
reaction per unit time and volume. In this paper we consider that reaction
rates are monomial functions of the species concentration. Therefore, the
deterministic kinetics of CRNs is described by sets of polynomial ordinary
differential equations. For instance, if the reaction rates of the 2 reaction
CRN $A_{1}+A_{2}\rightarrow A_{3}$, $A_{1}+A_{1}\rightarrow A_{1}$ are
$R_{1}=k_{1}x_{1}x_{2}$ and $R_{2}=k_{2}x_{1}^{2}$, the CRN kinetics is
described by
$\dot{x}_{1}=-k_{1}x_{1}x_{2}-k_{2}x_{1}^{2},\,\dot{x}_{2}=-k_{1}x_{1}x_{2},\,\dot{x}_{3}=k_{1}x_{1}x_{2}$.
Our models cover in fact all systems of ODEs whose right-hand sides are
integer-coefficient polynomials in $\boldsymbol{x}$ and $\boldsymbol{k}$,
because for any set of such polynomial ODEs there is a CRN whose kinetics is
described exactly by these equations. Monomial rate CRNs are thus more general
than mass-action networks [24], which satisfy an important constraint. In mass
action networks, the probability that two species react is proportional to
their abundances, and therefore the reaction rates are monomials containing
only concentrations of reactant species with exponents equal to the number of
molecules entering the reaction. This constraint leads to algebraic properties
of mass action CRNs exploited in chemical reaction network theory (CRNT),
which has been initiated by Horn, Jackson and Feinberg [25, 24].

Conservation laws are first integrals of CRNs, i.e. quantities that are
constant on any solution of the chemical kinetics ODEs. When they exist, they
can be used for model reduction by variable elimination [8]. A reduced model
can be more easily simulated and analysed. For example, chemical kinetic ODEs
with $n$ species and $n-1$ independent conservation laws are solvable by
quadratures. Conservation laws are also used for computing steady states of
CRNs. In the case of linear conservation laws, the set of steady states can be
partitioned into stoichiometric compatibility classes [25] that are
characterized by fixed, different values of the linear conservation laws. Here
we extend this definition to any type of conservation laws. Compatibility
classes are invariant sets of the chemical kinetics ODEs and, under some
conditions that we call “completeness”, contain a finite number of steady
states.

A lot of effort has been put into computing linear conservation laws [57, 61].
Non-linear conservation laws have been addressed using Darboux polynomials
[47]. In this paper we propose new methods for computing monomial and
polynomial conservation laws. Monomial conservation laws are computed using
methods based on stoichiometry. Polynomial conservation laws are computed
using syzygies, which are relations among the polynomial r.h.s of the chemical
kinetics ODEs. We also propose algorithmic solutions for testing the
independence and the completeness of the computed conservation laws.

Conservation laws lead to exact reduced models [8]. However, many model
reduction methods for CRNs are approximate. For instance, model reduction
methods based on singular perturbations for fast-slow CRNs are based on
elimination of the fast species [68, 27, 29, 53, 42, 17]. We introduce here
approximate conservation laws which is a key concept for approximate model
reduction by variable pooling [55, 3, 17]. We discuss the application of
approximate conservation laws to model reduction in [17].

##  2 Models

In this paper we consider CRNs with $n$ species $A_{1},\ldots,A_{n}$ whose
concentrations $\boldsymbol{x}=(x_{1},\ldots,x_{n})$ follow a system of ODEs
of the form

$\dot{x}_{1}=f_{1}(\boldsymbol{k},\boldsymbol{x}),\,\dots,\,\dot{x}_{n}=f_{n}(\boldsymbol{k},\boldsymbol{x})$
(1)

where

$f_{i}(\boldsymbol{k},\boldsymbol{x})=\sum_{j=1}^{r}S_{ij}k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}\in\mathbb{Z}[\boldsymbol{k},\boldsymbol{x}]=\mathbb{Z}[k_{1},\dots,k_{r},x_{1},\dots,x_{n}].$

The monomials
$\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}=x_{1}^{\alpha_{j1}}\cdots
x_{n}^{\alpha_{j_{n}}}$ appearing in the right hand sides of (1) are defined
by $r$ multi-indices
$\boldsymbol{\alpha}_{j}=(\alpha_{j1},\ldots,\alpha_{jn})\in\mathbb{N}^{n}$
and for each monomial the variable $k_{j}$ represents a rate constant. The
variables $\boldsymbol{k}=(k_{1},\dots,k_{r})$ take values in
$\mathbb{R}_{+}^{*}$ and the integer coefficients $S_{ij}$ form a matrix
$\boldsymbol{S}=(S_{ij})\in\mathbb{Z}^{n\times r}$ which is called the
stoichiometric matrix. We denote the vector of right hand sides of (1) by

$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=(f_{1}(\boldsymbol{k},\boldsymbol{x}),f_{2}(\boldsymbol{k},\boldsymbol{x}),\ldots,f_{n}(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}}.$

Mass action networks belong to this class of models. For a mass action
reaction

$\bar{\alpha}_{j1}A_{1}+\ldots+\bar{\alpha}_{jn}A_{n}\rightarrow\beta_{j1}A_{1}+\ldots+\beta_{jn}A_{n}$

we have $S_{ij}=\beta_{ji}-\bar{\alpha}_{ji}$ and the reaction rate is
$k_{j}\boldsymbol{x}^{\boldsymbol{\bar{\alpha}}_{j}}$, i.e. the stoichiometric
coefficients $\boldsymbol{\bar{\alpha}}_{j}$ and the multi-indices
$\boldsymbol{\alpha}_{j}$ coincide. However, this may not be the case in
general, meaning that the multi-indices $\boldsymbol{\alpha}_{j}$ and
$\boldsymbol{\bar{\alpha}}_{j}$ in $S_{ij}=\beta_{ji}-\bar{\alpha}_{ji}$ and
$k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}$ are not necessarily equal. The
more general type of CRNs introduced here allows us to cover also truncated
ODEs systems which we use in our reduction procedure in Section 3.

######  Theorem 1.

Any system of polynomial ODEs

$\dot{x}_{1}=f_{1}(\boldsymbol{k},\boldsymbol{x}),\,\dots,\,\dot{x}_{n}=f_{n}(\boldsymbol{k},\boldsymbol{x})$

where
$f_{i}(\boldsymbol{k},\boldsymbol{x})\in\mathbb{Z}[\boldsymbol{k},\boldsymbol{x}]$
is a polynomial in $\boldsymbol{k}$ and $\boldsymbol{x}$, homogeneous of
degree one in $\boldsymbol{k}=(k_{1},\dots,k_{r})$, can be written as (1) for
an appropriate matrix $\boldsymbol{S}$.

###### Proof.

A constructive proof is given by Algorithm 1. ∎

0: An ODE system whose r.h.s. are sums of rational monomials in
$\boldsymbol{k}$ and $\boldsymbol{x}$, homogeneous of degree one in
$\boldsymbol{k}$, with integer coefficients:
$\dot{x}_{i}=f_{i}(\boldsymbol{k},\boldsymbol{x})=\sum_{j=1}^{r_{i}}z_{ij}k_{ij}\boldsymbol{x}^{\boldsymbol{\alpha}_{ij}}$
where $z_{ij}\in\mathbb{Z}$, $\boldsymbol{\alpha}_{ij}\in\mathbb{Z}^{n}$,
$k_{ij}\in\\{k_{1},\dots,k_{r}\\}$, for $1\leq i\leq n$, $1\leq j\leq r_{i}$.

0: An integer coefficient matrix $\boldsymbol{S}$.

1: Compute a list of distinct monomials
$\\{M_{j}=k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}},1\leq j\leq
r\\}=\bigcup_{i=1}^{n}\bigcup_{j=1}^{r_{i}}\\{k_{ij}\boldsymbol{x}^{\boldsymbol{\boldsymbol{\alpha}}_{ij}}\\}.$
Two monomials $k_{ij}\boldsymbol{x}^{\boldsymbol{\alpha}_{ij}}$,
$k_{i^{\prime}j^{\prime}}\boldsymbol{x}^{\boldsymbol{\alpha}_{i^{\prime}j^{\prime}}}$
are distinct if
$\boldsymbol{\alpha}_{ij}\neq\boldsymbol{\alpha}_{i^{\prime}j^{\prime}}$ or
$k_{ij}\neq k_{i^{\prime}j^{\prime}}$.

2: for $i:=1$ to $n$ do

3: for $j:=1$ to $r$ do

4: $S_{ij}=0$

5: for $l:=1$ to $r_{i}$ do

6: if $k_{il}\boldsymbol{x}^{\boldsymbol{\alpha}_{il}}=M_{j}$ then

7: $S_{ij}=S_{ij}+z_{il}$

8: end if

9: end for

10: end for

11: end for

Algorithm 1 $\operatorname{Smatrix}$

######  Remark 2.

In many models, parameters are given by their numerical values, sometimes
linearly dependent over the integers. In this case, the independent parameters
and the stoichiometric matrix are hidden behind these numbers. We propose in
Algorithm 2 a way to still obtain descriptions of the form (1), where the
$k_{j}$ are linearly independent over the integers and the stoichiometric
coefficients are clearly visible. For example, the ODE system
$\dot{x}_{1}=1.5x_{1}x_{2},\,\dot{x}_{2}=x_{2}+3x_{1}x_{2}$ is transformed
into $\dot{x}_{1}=k_{1}x_{1}x_{2},\,\dot{x}_{2}=k_{2}x_{2}+2k_{1}x_{1}x_{2}$
with the stoichiometric matrix

$S=\begin{pmatrix}1&0\\\ 2&1\end{pmatrix}.$

0: An ODE system of the form:
$\dot{x_{i}}=\sum_{j=1}^{m}B_{ij}\boldsymbol{x}^{\boldsymbol{\boldsymbol{\alpha}}_{j}}$,
where $\boldsymbol{\alpha}_{j}\in\mathbb{Z}^{n}$, $B_{ij}\in\mathbb{R}$, $m$
is the number of distincts $\boldsymbol{\alpha}_{j}$.

0: An ODE system whose r.h.s. are sums of rational monomials in
$\boldsymbol{k}$ and $\boldsymbol{x}$, homogeneous of degree one in
$\boldsymbol{k}$, with integer coefficients:
$\dot{x}_{i}=f_{i}(\boldsymbol{k},\boldsymbol{x})=\sum_{j=1}^{r_{i}}z_{ij}k_{ij}\boldsymbol{x}^{\boldsymbol{\alpha}_{ij}}$
where $z_{ij}\in\mathbb{Z}$, $\boldsymbol{\alpha}_{ij}\in\mathbb{Z}^{n}$,
$k_{ij}\in\\{k_{1},\dots,k_{r}\\}$, for $1\leq i\leq n$, $1\leq j\leq r_{i}$.

1: for $h:=1$ to $m$ do

2: for $i:=1$ to $n$ do

3: for $j:=1$ to $i$ do

4: if $B_{ih}/B_{jh}\in\mathbb{Z}$ then

5: $z_{ih}:=\text{sign}(B_{ih})|B_{ih}/B_{jh}|$, $k_{ih}:=B_{ih}/z_{ih}$,
$z_{jh}:=\text{sign}(B_{jh})$, $k_{jh}:=B_{jh}/z_{jh}$

6: else if $B_{jh}/B_{ih}\in\mathbb{Z}$ then

7: $z_{jh}:=\text{sign}(B_{jh})|B_{jh}/B_{ih}|$, $k_{jh}:=B_{jh}/z_{jh}$,
$z_{ih}:=\text{sign}(B_{ih})$, $k_{ih}:=B_{ih}/z_{ih}$

8: else

9: $z_{jh}:=\text{sign}(B_{jh})$, $k_{jh}:=B_{jh}/z_{jh}$,
$z_{ih}:=\text{sign}(B_{ih})$, $k_{ih}:=B_{ih}/z_{ih}$

10: end if

11: end for

12: end for

13: end for

Algorithm 2 $\operatorname{CompareAndSplitCoefficients}$

##  3 Exact Versus Approximate Conservation Laws

We start this section with a motivating example.

######  Example 3.

Let us consider the irreversible Michaelis-Menten mechanism that is
paradigmatic for enzymatic reactions. We choose rate constants corresponding
to the so-called quasi-equilibrium, studied by Michaelis and Menten. The
reaction network for this model is

$S+E\underset{k_{2}}{\overset{k_{1}}{\rightleftharpoons}}ES\overset{k_{3}\delta}{\longrightarrow}E+P,$

where $S$ is a substrate, $E$ is an enzyme, $ES$ is an enzyme-substrate
complex and $k_{1}$,$k_{2}$,$k_{3}$ are rate constants of the same order of
magnitude. Here $0<\delta<1$ is a small positive scaling parameter, indicating
that the third rate constant is small.

According to mass-action kinetics, the concentrations $x_{1}=[S]$,
$x_{2}=[SE]$ and $x_{3}=[E]$ satisfy the system of ODEs

$\displaystyle\dot{x_{1}}$ $\displaystyle=$ $\displaystyle-
k_{1}x_{1}x_{3}+k_{2}x_{2},$ (2a) $\displaystyle\dot{x_{2}}$ $\displaystyle=$
$\displaystyle k_{1}x_{1}x_{3}-k_{2}x_{2}-\delta k_{3}x_{2},$ (2b)
$\displaystyle\dot{x_{3}}$ $\displaystyle=$ $\displaystyle-
k_{1}x_{1}x_{3}+k_{2}x_{2}+\delta k_{3}x_{2}.$ (2c)

We consider the truncated system of ODEs

$\displaystyle\dot{x_{1}}$ $\displaystyle=$ $\displaystyle-
k_{1}x_{1}x_{3}+k_{2}x_{2},$ (3a) $\displaystyle\dot{x_{2}}$ $\displaystyle=$
$\displaystyle k_{1}x_{1}x_{3}-k_{2}x_{2},$ (3b) $\displaystyle\dot{x_{3}}$
$\displaystyle=$ $\displaystyle-k_{1}x_{1}x_{3}+k_{2}x_{2},$ (3c)

which is obtained by setting $\delta=0$ in (2). If $k_{i},1\leq i\leq 3$ are
of order ${O}(\delta^{0})$, then the truncated system (3) describes the
dynamics of the model on fast timescales of order ${O}(\delta^{0})$.

The steady state of the fast dynamics is obtained by equating to zero the
r.h.s. of (3). The resulting condition is called quasi-equilibrium (QE),
because it means that the complex formation rate $k_{1}x_{1}x_{3}$ is equal to
the complex dissociation rate $k_{2}x_{2}$. In other words, the reversible
reaction

$S+E\underset{k_{2}}{\overset{k_{1}}{\rightleftharpoons}}ES$

is at equilibrium. The QE condition is reached only at the end of the fast
dynamics and is satisfied with a precision of order ${O}(\delta)$ during the
slow dynamics [31].

We introduce the linear combinations of variables $x_{4}=x_{1}+x_{2}$ and
$x_{5}=x_{2}+x_{3}$, corresponding to the total substrate and total enzyme
concentrations, respectively. Addition of equation (2b) and (2c) leads to
$\dot{x_{5}}=0$ which means that for solutions of the full system, $x_{5}$ is
constant for all times. We will call such a quantity an “exact conservation
law”. The addition of equation (3a) and (3b) and the addition of equation (3b)
and (3c) lead to $\dot{x_{4}}=0$ and $\dot{x_{5}}=0$, respectively. This means
that $x_{4}$ is constant for solutions of the truncated dynamics, valid at
short times $t={O}(\delta^{0})$, but may not be constant at larger times
$t={O}(\delta^{-1})$. We call such a quantity an “approximate conservation
law”. The quantity $x_{5}$ is both an exact and approximate conservation law,
because is conserved by both (3) and (2) and therefore it is constant at all
times.

More generally, let us assume that the r.h.s. of the system in (1) can be
decomposed into a sum of two terms, the first one dominating (being much
larger than the second). This decomposition can be obtained by rescaling the
model variables and parameters by powers of a small scaling parameter
$\delta$, where $0<\delta<1$. To this aim, we define
$\boldsymbol{\bar{k}}=(\bar{k}_{1},\ldots,\bar{k}_{r})$ and
$\boldsymbol{\bar{x}}=(\bar{x}_{1},\ldots,\bar{x}_{m})$ such that
$k_{i}=\delta^{e_{i}}\bar{k}_{i}$, $x_{i}=\delta^{d_{i}}\bar{x}_{i}$ and
obtain

$f_{i}(\boldsymbol{k},\boldsymbol{x})=f_{i}^{(1)}(\boldsymbol{k},\boldsymbol{x})+f_{i}^{(2)}(\boldsymbol{k},\boldsymbol{x}),$
(4)

where

$f_{i}^{(1)}(\boldsymbol{k},\boldsymbol{x})=\delta^{b_{i}}\bar{f}_{i}^{(1)}(\boldsymbol{\bar{k}},\boldsymbol{\bar{x}})\quad\mathrm{and}\quad
f_{i}^{(2)}(\boldsymbol{k},\boldsymbol{x})=\delta^{b^{\prime}_{i}}\bar{f}_{i}^{(2)}(\boldsymbol{\bar{k}},\boldsymbol{\bar{x}},\delta)$

with $b_{i}<b^{\prime}_{i}$ (for details see [42, 17]).

For a small parameter $\delta$, the terms
$\delta^{b^{\prime}_{i}}\bar{f}_{i}^{(2)}(\boldsymbol{\bar{k}},\boldsymbol{\bar{x}},\delta)$
are dominated by the terms
$\delta^{b_{i}}\bar{f}_{i}^{(1)}(\boldsymbol{\bar{k}},$
$\boldsymbol{\bar{x}})$. This justifies to introduce the truncated system as
the system of ODEs obtained by keeping only the lowest order, dominant terms
of the system in (1), namely

$\dot{x}_{1}=f_{1}^{(1)}(\boldsymbol{k},\boldsymbol{x}),\,\ldots,\,\dot{x}_{n}=f_{n}^{(1)}(\boldsymbol{k},\boldsymbol{x}).$
(5)

We denote the vector of right hand sides of the truncated system by

$\boldsymbol{F}^{(1)}(\boldsymbol{k},\boldsymbol{x})=(f_{1}^{(1)}(\boldsymbol{k},\boldsymbol{x}),\ldots,f_{n}^{(1)}(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}}.$

######  Remark 4.

Depending on the application, the truncated system in (5) can consist of less
ODEs than the system in (1). For instance, in the case of model reduction of
fast-slow CRNs, only the ODEs corresponding to fast variables are included in
the truncated system. In this case, the truncated system describes the fast
timescale dynamics.

######  Definition 5.

  1. (a)

A function $\phi(\boldsymbol{k},\boldsymbol{x})$ is a parametric exact
conservation law if there exists a non-empty semi-algebraic set
$V_{\boldsymbol{k}}\subseteq\mathbb{R}^{r}_{+}$ such that for all
$\boldsymbol{k}\in V_{\boldsymbol{k}}$, $\phi(\boldsymbol{k},\boldsymbol{x})$
is a first integral of the full system in (1), i.e. if

$\forall\boldsymbol{k}\in
V_{\boldsymbol{k}}\quad\forall\boldsymbol{x}\in\mathbb{R}^{n}_{+}\quad(\sum_{i=1}^{n}\frac{\partial\phi}{\partial
x_{i}}(\boldsymbol{k},\boldsymbol{x})f_{i}(\boldsymbol{k},\boldsymbol{x})=0).$

  2. (b)

A function $\phi(\boldsymbol{k},\boldsymbol{x})$ is a parametric approximate
conservation law if there exists a non-empty semi-algebraic set
$V_{\boldsymbol{k}}\subseteq\mathbb{R}^{r}_{+}$ such that for all
$\boldsymbol{k}\in V_{\boldsymbol{k}}$, $\phi(\boldsymbol{k},\boldsymbol{x})$
is a first integral of the truncated system in (5), i.e. if

$\forall\boldsymbol{k}\in
V_{\boldsymbol{k}}\quad\forall\boldsymbol{x}\in\mathbb{R}^{n}_{+}\quad(\sum_{i=1}^{n}\frac{\partial\phi}{\partial
x_{i}}(\boldsymbol{k},\boldsymbol{x})f_{i}^{(1)}(\boldsymbol{k},\boldsymbol{x})=0).$

  3. (c)

A parametric exact (approximate) conservation law of the form

$\phi(\boldsymbol{k},\boldsymbol{x})=c_{1}(\boldsymbol{k})x_{1}+\dots+c_{n}(\boldsymbol{k})x_{n}$

with $c_{i}(\boldsymbol{k})\in\mathbb{R}[\boldsymbol{k}]$ is called an exact
(approximate) linear conservation law. If $c_{i}\geq 0$ for $1\leq i\leq n$,
the linear conservation law is called semi-positive.

  4. (d)

A parametric exact (approximate) conservation law of the form
$\phi(\boldsymbol{k},\boldsymbol{x})=a(\boldsymbol{k})x_{1}^{m_{1}}$ $\cdots
x_{n}^{m_{n}}$ with $m_{i}\in\mathbb{Z}$ and
$a(\boldsymbol{k})\in\mathbb{R}[\boldsymbol{k}]$, is called an exact
(approximate) rational monomial conservation law. If $m_{i}\in\mathbb{Z}_{+}$
for $1\leq i\leq n$, the conservation law is called monomial. For simplicity,
in this paper, we will call both types monomial.

  5. (e)

A parametric exact (approximate) conservation law of the form

$\sum_{i=1}^{s}a_{i}(\boldsymbol{k})x_{1}^{m_{1i}}\cdots x_{n}^{m_{ni}}$

with $m_{ji}\in\mathbb{Z}_{+}$ and
$a_{i}(\boldsymbol{k})\in\mathbb{R}[\boldsymbol{k}]$ is called a parametric
exact (approximate) polynomial conservation law.

######  Remark 6.

In some applications we are interested in conservation laws
$\phi(\boldsymbol{x})$ that do not depend on $\boldsymbol{k}$. For instance,
in model reduction of a fast-slow CRN, approximate conservation laws that are
not exact represent variables of the reduced model [55, 3, 31, 53, 17]. When
these conservation laws depend only on $\boldsymbol{x}$, one can show that
these approximate conservation laws are always slower than the variables which
compose them, a prerequisite for fast-slow reduction [17]. Moreover, in such a
case the reduced model is a CRN obtained from the initial network by graph
rewriting operations, i.e by pooling and pruning (see [53]).

######  Remark 7.

Following Remark 6, in this paper, we use the following convention. An
unconditional conservation law is a function $\phi(\boldsymbol{x})$ that does
not depend on $\boldsymbol{k}$ and satisfies the conditions of Definition 5,
while a parametric conservation law is a function
$\phi(\boldsymbol{k},\boldsymbol{x})$ that does depend on $\boldsymbol{k}$.

##  4 Irreducibility, Independence and Completeness

Some variables $x_{i}$ may not appear in a conservation law
$\phi(\boldsymbol{k},\boldsymbol{x})$. This means, in the case of a linear
conservation law, that some coefficients $c_{i}$ may be zero or, in the case
of a nonlinear conservation law, that some partial derivatives
$\frac{\partial\phi}{\partial x_{i}}(\boldsymbol{k},\boldsymbol{x})$ may
vanish. If $r$ is the number of all non-zero quantities
$\frac{\partial\phi}{\partial x_{i}}(\boldsymbol{k},\boldsymbol{x})$, then we
say that the conservation law depends on $r$ variables.

######  Definition 8.

An exact or approximate conservation law depending on $r$ variables is called
reducible if it can be split into the sum or the product of two conservation
laws such that at least one of them depends on a number of variables
$r^{\prime}$ with $1\leq r^{\prime}<r$. Otherwise it is called irreducible.

######  Definition 9.

For $\boldsymbol{k}\in\mathbb{R}^{r}_{+}$ a steady state $\boldsymbol{x}$ is a
positive solution of $\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$. We
denote the steady state variety by $\mathcal{S}_{\boldsymbol{k}}$. A steady
state is degenerate or non-degenerate if the Jacobian
$D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})$ is singular
or regular, respectively.

Degeneracy of steady states implies that $\mathcal{S}_{\boldsymbol{k}}$ is not
discrete. Reciprocally, we have the following result.

######  Theorem 10.

Assume that for $\boldsymbol{k}\in\mathbb{R}^{r}_{+}$,
$\mathcal{S}_{\boldsymbol{k}}$ is a manifold. If the local dimension at a
point $\boldsymbol{x}_{0}\in\mathcal{S}_{\boldsymbol{k}}$ is strictly
positive, then $\boldsymbol{x}_{0}$ is degenerate.

###### Proof.

There is a neighborhood $U\subset\mathbb{R}$ of zero and a smooth function
$U\rightarrow\mathcal{S}_{\boldsymbol{k}}$,
$\alpha\mapsto\boldsymbol{x}(\alpha)$ with
$\boldsymbol{x}(0)=\boldsymbol{x}_{0}$. By differentiating
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x}(\alpha))=0$ with respect to
$\alpha$ at $0$ we get
$D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x}_{0})\frac{d\boldsymbol{x}(0)}{d\alpha}=0$,
and thus $D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x}_{0})$
is singular. ∎

######  Definition 11.

A set of $s$ exact conservation laws defined by the vector

$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=(\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\ldots,\phi_{s}(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}}$

is called complete if the Jacobian matrix

$\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{k},\boldsymbol{x})=D_{\boldsymbol{x}}(\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x}),\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}}$

has rank $n$ for any
$\boldsymbol{k}\in\mathbb{R}_{+}^{r},\boldsymbol{x}\in\mathbb{R}_{+}^{n}$
satisfying $\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$. The set is
called independent if the Jacobian matrix of
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})^{\mathsf{T}}$ with respect
to $\boldsymbol{x}$ has rank $s$ for any
$\boldsymbol{k}\in\mathbb{R}_{+}^{r},\boldsymbol{x}\in\mathbb{R}_{+}^{n}$ such
that $\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$. In the case of a set
of approximate conservation laws defined by a vector
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})$, completeness is defined
with $\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})$ replaced by
$\boldsymbol{F}^{(1)}(\boldsymbol{k},\boldsymbol{x})$.

######  Proposition 12.

If a system has a complete set of conservation laws, then the intersection of
$\mathcal{S}_{\boldsymbol{k}}$ with
$\\{\boldsymbol{x}\mid\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{c}_{0}\\}\cap\mathbb{R}_{+}^{n}$,
where $\boldsymbol{c}_{0}\in\mathbb{R}^{r}$, is finite.

###### Proof.

Suppose that the intersection is non-empty and contains $\boldsymbol{x}$.
Because the rank of
$\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{k},\boldsymbol{x})$
in $\boldsymbol{x}$ is $n$, the implicit function theorem implies that, for
$\boldsymbol{k}$ and $\boldsymbol{c}_{0}$, $\boldsymbol{x}$ is isolated from
other solutions of $\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$,
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{c}_{0}$. As
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})$ is polynomial, the
intersection is finite. ∎

######  Remark 13.

When $\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})$ is linear in
$\boldsymbol{x}$ and results from a stoichiometric matrix, the set

$\mathcal{C}_{\boldsymbol{c}_{0}}=\\{\boldsymbol{x}\mid\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{c}_{0}\\}\cap\mathbb{R}_{+}^{n}$

is called stoichiometric compatibility class or reaction simplex [70, 25].
Stoichiometric compatibility classes of systems with complete sets of linear
conservation laws contain a finite number of steady states. We note that some
authors call completeness of linear conservation laws “non-degeneracy” [26].
Here, we extend these definitions and call compatibility classes all sets of
the type $\mathcal{C}_{\boldsymbol{c}_{0}}$, even if
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})$ is a non-linear
conservation law. Compatibility classes of polynomial conservation laws are
semi-algebraic sets.

######  Remark 14.

Since our concern is the number of positive solutions in $\boldsymbol{x}$ of
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$,
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{c}_{0}$, it
would be more natural to consider in Definition 11 the rank of
$\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{k},\boldsymbol{x})$
on

$\mathcal{S}_{\boldsymbol{k}}\cap\\{\boldsymbol{x}\mid\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{c}_{0}\\}\cap\mathbb{R}_{+}^{n}.$

In fact, as this rank does not depend on $\boldsymbol{c}_{0}$, it is simpler
and equivalent to impose its value on
$\mathcal{S}_{\boldsymbol{k}}\cap\mathbb{R}_{+}^{n}$.

######  Remark 15.

The independent linear conservation laws

$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=(x_{1}+x_{2},x_{2}+x_{3})$

of Example 3 are complete. More precisely, the Jacobian of

$(\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x}),\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}},$

where $\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})$ is the vector of right
hand sides of (3), has a $3\times 3$ minor

$M=\mathrm{det}(D_{\boldsymbol{x}}(-k_{1}x_{1}x_{3}+k_{2}x_{2},\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}})=-k_{2}-k_{1}x_{1}-k_{1}x_{3}.$

This minor can not be zero for positive $\boldsymbol{x}$, $\boldsymbol{k}$ on
the steady state variety which is defined by the single equation
$-k_{1}x_{1}x_{3}+k_{2}x_{2}=0$. Therefore the rank of
$\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{x},\boldsymbol{k})$
is $3$.

Furthermore, all stoichiometric compatibility classes defined by
$x_{1}+x_{2}=c_{01}$, $x_{2}+x_{3}=c_{02}$, $\boldsymbol{x}>0$ contain a
unique steady state

$\displaystyle x_{1}$ $\displaystyle=$
$\displaystyle\big{(}k_{1}(c_{01}-c_{02})-k_{2}+\sqrt{\Delta}\big{)}/(2k_{1}),$
$\displaystyle x_{2}$ $\displaystyle=$
$\displaystyle\big{(}k_{1}(c_{01}+c_{02})+k_{2}-\sqrt{\Delta}\big{)}/(2k_{1}),$
$\displaystyle x_{3}$ $\displaystyle=$
$\displaystyle\big{(}k_{1}(c_{02}-c_{01})-k_{2}+\sqrt{\Delta}\big{)}/(2k_{1}),$

where
$\Delta=(c_{01}-c_{02})^{2}k_{1}^{2}+k_{2}^{2}+2k_{1}k_{2}(c_{01}+c_{02})$.

The notions of completeness and independence in Definition 11 are effective.
Indeed, Algorithm 4 and 5 below check for completeness and independence,
respectively. Both algorithms use a parametric rank computation (see in both
cases line 1 with the call ”ParametricRank”), which we are going to present
first.

Let $P$, $N$, $M\in\mathbb{N}$ and let $\mathbb{R}[\boldsymbol{v}]$ be a
polynomial ring in $P$ indeterminates $\boldsymbol{v}=(v_{1},\dots,v_{P})$.
Furthermore, let $A\in\mathbb{R}[\boldsymbol{v}]^{M\times N}$ be a matrix
whose entries are polynomials in $\boldsymbol{v}$. We are interested in the
rank of $A$ in dependence of different values
$\bar{\boldsymbol{v}}\in\mathbb{R}^{P}$ for the variables $\boldsymbol{v}$.
Algorithm 3 (ParametricRank) provides a decomposition of the affine space
$\mathbb{R}^{P}$ into disjoint semialgebraic subsets

$S_{1},\dots,S_{I}$

such that for each subset $S_{i}$ the rank of $A(\bar{\boldsymbol{v}})$ is
constant for all points $\bar{\boldsymbol{v}}$ of $S_{i}$. The output of the
algorithm is organized as a set of pairs

$R=\\{(\Gamma_{1},r_{1}),\,\dots,\,(\Gamma_{n},r_{n})\\}.$

For each $i=1,\dots,I$, $\Gamma_{i}$ is a conjunction of polynomial
constraints in $\boldsymbol{v}$, i.e. polynomial equations and inequations in
$\boldsymbol{v}$ with coefficients in $\mathbb{R}$, and $r_{i}$ is an integer.
The semialgebraic sets $S_{1},\dots,S_{I}$ from above are defined by
$\Gamma_{1},\dots,\Gamma_{n}$ and $r_{1},\dots,r_{I}$ are the ranks of $A$ on
these sets.

In order to determine the different ranks the algorithm computes row echelon
forms by introducing case distinctions. Roughly speaking the algorithm works
in the following way. For a semialgebraic set $S$ defined by a conjunction of
polynomial constraints $\Gamma$ and for a matrix $A$ which is in row echelon
form up to the $p$-th row, the algorithm checks if there is an entry $A_{ij}$
in $A$, where $p+1\leq i\leq M$ and $p+1\leq j\leq N$, which does not vanish
for any point of $S$. If $A_{mn}$ is such an entry, then one swaps rows $p+1$
with $m$ and columns $p+1$ with $n$. One uses then the entry $A_{p+1,p+1}$ in
the new matrix as a pivot element to delete all subsequent entries and
introduces a new case consisting of the same set $S$, i.e. the same $\Gamma$,
but now with the matrix in row echelon form up to row $p+1$. If there is not
such an entry, the algorithm checks if there is an entry $A_{ij}$, where $i$,
$j$ runs as above, which does not vanish for all points of $S$. If $A_{mn}$ is
such an entry, the algorithm introduces two new cases. Both cases consist of
the old matrix $A$ in row echelon form up to row $p$. The difference is now
that one case consists of the semialgebraic set defined by $\Gamma\land
A_{mn}\neq 0$ and the other one consists of the set defined by $\Gamma\land
A_{mn}=0$. The new conditions guarantee the non-vanishing and the vanishing of
the entry $A_{mn}$ in the respective case. If there is not an entry $A_{mn}$
as described above, on $S$ the matrix $A$ is in row echelon form and has rank
$p$. The algorithm returns $\Gamma$ and $p$. The algorithm processes a stack
of cases in the way described above starting with the case which consists of
the semialgebraic set $\mathbb{R}^{P}$ defined by
$\Gamma="\operatorname{true}"$ and the matrix $A$ in row echelon form up to
row $p=0$.

Algorithm 3 ParametricRank

0: An $M\times N$-matrix $A(\boldsymbol{v})$ with polynomial entries in
parameters $\boldsymbol{v}$.

0: A set of pairs
$\\{(\Gamma_{1}(\boldsymbol{v}),r_{1}),\dots,(\Gamma_{I}(\boldsymbol{v}),r_{I})\\}$,
where each $\Gamma_{i}(\boldsymbol{v})$ is a conjunction of polynomial
equations and inequations and $r_{i}\in\\{1,\dots,N\\}$. For any real choice
$\bar{\boldsymbol{v}}$ of parameters $\boldsymbol{v}$ there is one and only
one $i\in\\{1,\dots,I\\}$ such that
$\mathbb{R}\models\Gamma_{i}(\bar{\boldsymbol{v}})$. For this $i$ we have
$\rank A(\bar{\boldsymbol{v}})=r_{i}$.

1: $I:=0$

2: create an empty stack

3: $\operatorname{push}~{}(\operatorname{true},A,0)$

4: while stack is not empty do

5: $(\Gamma,A,p):=\operatorname{pop}$

6: if $\Gamma\nvdash\operatorname{false}$ then

7: if there is $m\in\\{p+1,\dots,M\\}$, $n\in\\{p+1,\dots,N\\}$ s.t.
$\Gamma\vdash A_{mn}\neq 0$ then

8: in $A$, swap rows $p+1$ with $m$ and columns $p+1$ with $n$

9: in $A$, use row $p+1$ to obtain $A_{p+2,p+1}=\dots=A_{M,p+1}=0$

10: $\operatorname{push}~{}(\Gamma,A,p+1)$

11: else if there is $m\in\\{p+1,\dots,M\\}$, $n\in\\{p+1,\dots,N\\}$ s.t.
$\Gamma\nvdash A_{mn}=0$ then

12: $a:=A_{mn}$

13: $\operatorname{push}~{}(\Gamma\land a\neq 0,A,p)$

14: in $A$, set $A_{mn}:=0$optional optimisation

15: $\operatorname{push}~{}(\Gamma\land a=0,A,p)$

16: else $A$ is in row echelon form modulo $\Gamma$

17: $I:=I+1$

18: $(\Gamma_{I},r_{I}):=(\Gamma,p)$

19: end if

20: end if

21: end while

22: return $\\{(\Gamma_{1},r_{1}),\dots,(\Gamma_{I},r_{I})\\}$

######  Example 16.

We consider the $3\times 3$-matrix

$A^{(0)}=\begin{pmatrix}0&v_{1}&1\\\ v_{1}&v_{3}^{2}-1+v_{1}v_{2}&v_{3}\\\
v_{3}-1&v_{3}+1+v_{1}v_{3}^{3}&v_{3}^{3}\end{pmatrix}$

with polynomial entries in $\mathbb{R}[\boldsymbol{v}]$ with
$\boldsymbol{v}=(v_{1},v_{2},v_{3})$. Starting with
$(\operatorname{true},A^{(0)},0)$ we find in $A^{(0)}$ the entry
$A^{(0)}_{1,3}=1$ which clearly does not vanish on $\mathbb{R}^{3}$. We swap
the first and third column and use $1$ as a pivot to delete the entries
$v_{3}$ and $v_{3}^{3}$. We obtain

$A^{(1)}=\begin{pmatrix}1&v_{1}&0\\\ 0&v_{3}^{2}-1&v_{1}\\\
0&v_{3}+1&v_{3}-1\end{pmatrix}$

and introduce $(\operatorname{true},A^{(1)},1)$. For each entry
$A^{(1)}_{i,j}$ with $2\leq i,j\leq 3$ there is at least one point in
$\mathbb{R}^{3}$ such that $A^{(1)}_{i,j}$ vanishes. On the other hand for
each entry $A^{(1)}_{i,j}$ with $2\leq i,j\leq 3$ there is also at least one
point in $\mathbb{R}^{3}$ such that $A^{(1)}_{i,j}$ does not vanish. We choose
$A^{(1)}_{2,2}=v_{3}^{2}-1$ and introduce the cases $(v_{3}^{2}-1\neq
0,A^{(1)},1)$ and $(v_{3}^{2}-1=0,A^{(1)},1)$. Continuing with the first case,
we trivially find that the entry $A^{(1)}_{2,2}=v_{3}^{2}-1$ does not vanish
on the semialgebraic set
$\\{\boldsymbol{\overline{v}}\in\mathbb{R}^{3}\mid\overline{v}_{3}^{2}-1\neq
0\\}$ and we use it as a pivot to delete the entry $v_{3}+1$. We obtain

$A^{(2)}=\begin{pmatrix}1&v_{1}&0\\\ 0&v_{3}^{2}-1&v_{1}\\\
0&0&(v_{3}-1)^{2}-v_{1}\end{pmatrix}$

and introduce $(v_{3}^{2}-1\neq 0,A^{(2)},2)$. Since the set
$\\{\boldsymbol{\overline{v}}\in\mathbb{R}^{3}\mid\overline{v}_{3}^{2}-1\neq
0\\}$ contains points where $A^{(2)}_{3,3}$ vanishes and points where
$A^{(2)}_{3,3}$ does not vanish, we introduce the two cases

$\displaystyle(v_{3}^{2}-1\neq 0\,\land\,(v_{3}-1)^{2}\neq v_{1},A^{(2)},2),$
$\displaystyle(v_{3}^{2}-1\neq 0\,\land\,(v_{3}-1)^{2}=v_{1},A^{(2)},2).$

In the first case, the entry $A^{(2)}_{3,3}$ clearly does not vanish on the
semialgebraic set

$\\{\boldsymbol{\overline{v}}\in\mathbb{R}^{3}\mid\overline{v}_{3}^{2}-1\neq
0\,\land\,(\overline{v}_{3}-1)^{2}\neq\overline{v}_{1}\\}.$

Since there is no entry left to delete with the pivot $A^{(2)}_{3,3}$, we
simply introduce the new case

$(v_{3}^{2}-1\neq 0\,\land\,(v_{3}-1)^{2}\neq v_{1},A^{(2)},3),$

which yields the output $(v_{3}^{2}-1\neq 0\,\land\,(v_{3}-1)^{2}\neq
v_{1},3)$. In the second case we have that the entry $A^{(2)}_{3,3}$ is zero
on the constructed semialgebraic set and so the algorithm returns
$(v_{3}^{2}-1\neq 0\,\land\,(v_{3}-1)^{2}=v_{1},2)$. Finally, we are coming
back to the case $(v_{3}^{2}-1=0,A^{(1)},1)$. On the semialgebraic set
$\\{\boldsymbol{\overline{v}}\in\mathbb{R}^{3}\mid\overline{v}_{3}^{2}-1=0\\}$
the entries $v_{3}^{2}-1$, $v_{3}+1$ and $v_{3}-1$ are zero while there are
points in this set where the entry $A^{(1)}_{2,3}=v_{1}$ vanishes and points
where it does not vanish. We introduce the cases
$(v_{3}^{2}-1=0\,\land\,v_{1}\neq 0,A^{(1)},1)$ and
$(v_{3}^{2}-1=0\,\land\,v_{1}=0,A^{(1)},1)$. In the first case the entry
$A^{(1)}_{2,3}=v_{1}$ is nonzero on the semialgebraic set

$\\{\boldsymbol{\overline{v}}\in\mathbb{R}^{3}\mid\overline{v}_{3}^{2}-1=0\,\land\,\overline{v}_{1}\neq
0\\}$

and so we obtain $(v_{3}^{2}-1=0\,\land\,v_{1}\neq 0,A^{(1)},2)$, which leads
to the output $(v_{3}^{2}-1=0\,\land\,v_{1}\neq 0,2)$. In the second case all
entries $A^{(1)}_{i,j}$ with $2\leq i,j\leq 3$ are zero on the corresponding
semialgebraic set and therefore the algorithm returns
$(v_{3}^{2}-1=0\,\land\,v_{1}=0,1)$.

Let us get back to checking for completeness and independence. We start with
Algorithm 4 below which handles completeness.

0: 1\.
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=(f_{1}(\boldsymbol{k},\boldsymbol{x})),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))$;
2\.
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=(\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\dots\phi_{s}(\boldsymbol{k},\boldsymbol{x}))$;
3\. $\boldsymbol{k}=(k_{1},\dots,k_{r})$; 4\.
$\boldsymbol{x}=(x_{1},\dots,x_{n})$

0: “yes” if $\Phi$ is complete according to Definition 11, “no” otherwise

1:
$R:=\operatorname{ParametricRank}\bigl{(}\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{k},\boldsymbol{x})\bigr{)}$

2: $R_{n}:=\\{\,(\Gamma,r)\in R\mid r=n\,\\}$

3: $\varrho_{n}:=\bigvee_{(\Gamma,n)\in R_{n}}\Gamma$

4:
$\gamma:=\forall\,\boldsymbol{k}\,\forall\,\boldsymbol{x}\,(\boldsymbol{k}>0\,\land\,\boldsymbol{x}>0\,\land\,\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0\longrightarrow\varrho_{n})$

5: if $\mathbb{R}\models\gamma$ then

6: return “yes”

7: else

8: return “no”

9: end if

Algorithm 4 $\operatorname{IsComplete}$.

On the grounds of Algorithm 3, we construct in line 2–3 in Algorithm 4 an
equivalent logic condition $\varrho_{n}$ for
$\rank\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{k},\boldsymbol{x})=n$.
In line 4 we construct $\gamma$ as a direct formalization of the definition of
completeness and in line 5 we finally test validity of $\gamma$ over the
reals. Technically we use a combination of various effective quantifier
elimination procedures for the theory of real closed fields [21, 64, 65, and
the references there] combined with heuristic simplification techniques [20]
in the Redlog system [18, 40, 58]. Alternatively, one could use Satisfiability
Modulo Theories solving over the logic QF_NRA [51, 4, 1]. Algorithm 5 proceeds
analogously to Algorithm 4 and tests for independence.

0: 1\.
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=(f_{1}(\boldsymbol{k},\boldsymbol{x})),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))$;
2\.
$\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})=(\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\dots\phi_{s}(\boldsymbol{k},\boldsymbol{x}))$;
3\. $\boldsymbol{k}=(k_{1},\dots,k_{r})$; 4\.
$\boldsymbol{x}=(x_{1},\dots,x_{n})$

0: “yes” if $\Phi$ is independent according to Definition 11, “no” otherwise

1:
$R:=\operatorname{ParametricRank}(D_{\boldsymbol{x}}\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})^{T})$

2: $R_{s}:=\\{\,(\Gamma,r)\in R\mid r=s\,\\}$

3: $\varrho_{s}:=\bigvee_{(\Gamma,s)\in R_{s}}\Gamma$

4:
$\iota:=\forall\,\boldsymbol{k}\,\forall\,\boldsymbol{x}\,(\boldsymbol{k}>0\,\land\,\boldsymbol{x}>0\,\land\,\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0\longrightarrow\varrho_{s})$

5: if $\mathbb{R}\models\iota$ then

6: return “yes”

7: else

8: return “no”

9: end if

Algorithm 5 $\operatorname{IsIndependent}$

######  Example 17.

We automatically process Example 3 with our Algorithm 4

$\operatorname{IsComplete}$. Our input is

$\displaystyle\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})$ $\displaystyle=$
$\displaystyle(-k_{1}x_{1}x_{3}+k_{2}x_{2},k_{1}x_{1}x_{3}-k_{2}x_{2},-k_{1}x_{1}x_{3}+k_{2}x_{2}),$
$\displaystyle\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})$
$\displaystyle=$ $\displaystyle(x_{1}+x_{2},x_{2}+x_{3})$

and $\boldsymbol{k}=(k_{1},k_{2})$, $\boldsymbol{x}=(x_{1},x_{2},x_{3})$. We
obtain the parametric Jacobian

$\boldsymbol{J}_{\boldsymbol{F},\boldsymbol{\phi}}(\boldsymbol{k},\boldsymbol{x})=\begin{pmatrix}-k_{1}x_{3}&k_{2}&-k_{1}x_{1}\\\
k_{1}x_{3}&-k_{2}&k_{1}x_{1}\\\ -k_{1}x_{3}&k_{2}&-k_{2}x_{1}\\\ 1&1&0\\\
0&1&1\end{pmatrix}$

and compute in line 1 the parametric rank

$R=\\{(k_{1}x_{1}+k_{1}x_{3}+k_{2}\neq
0,3),(k_{1}x_{1}+k_{1}x_{3}+k_{2}=0,2)\\},$

from which we select $R_{3}=\\{(k_{1}x_{1}+k_{1}x_{3}+k_{2}\neq 0,3)\\}$ in
line 2 and

$\varrho_{3}=(k_{1}x_{1}+k_{1}x_{3}+k_{2}\neq 0)$

in line 3. Completeness is straightforwardly formalized in line 4 by

$\gamma=\forall\,\boldsymbol{k}\,\forall\,\boldsymbol{x}\,(\boldsymbol{k}>0\,\land\,\boldsymbol{x}>0\,\land\,k_{1}x_{1}x_{3}-k_{2}x_{2}=0\longrightarrow
k_{1}x_{1}+k_{1}x_{3}+k_{2}\neq 0),$

where some redundant equations are automatically removed from
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$ via the standard simplifier
described in [20]. Real quantifier elimination on $\gamma$ in line 5
equivalently yields “$\operatorname{true}$”, which confirms completeness, and
we return “yes” in line 6.

Similarly, Algorithm 5 $\operatorname{IsIndependent}$ computes

$D_{\boldsymbol{x}}\boldsymbol{\Phi}(\boldsymbol{k},\boldsymbol{x})^{\mathsf{T}}=\begin{pmatrix}1&1&0\\\
0&1&1\end{pmatrix},$

$R=\\{(\operatorname{true},2)\\}$, $R_{2}=R$ and
$\varrho_{2}=\operatorname{true}$. This yields

$\iota=\forall\,\boldsymbol{k}\,\forall\,\boldsymbol{x}\,(\boldsymbol{k}>0\,\land\,\boldsymbol{x}>0\,\land\,k_{1}x_{1}x_{3}-k_{2}x_{2}=0\longrightarrow\operatorname{true}),$

which by quantifier elimination is equivalent to “$\operatorname{true}$”, and
we return “yes”. Recall that any implication with “$\operatorname{true}$” on
the right hand side holds already due to Boolean logic [59]. In our framework
here, this corresponds to the observation that completeness and independence
hold whenever we encounter full rank of the corresponding Jacobian for all
choices of parameters. The automatic computations described here take less
than 0.01 sec altogether.

The next example shows that linear conservation laws are not always complete.

######  Example 18.

One checks easily that the system

$\dot{x}_{1}=1-x_{1}-x_{2},\quad\dot{x}_{2}=x_{1}+x_{2}-1$

has the linear conservation law
$\Phi(\boldsymbol{k},\boldsymbol{x})=x_{1}+x_{2}$ and that the Jacobian of
$(\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x}),$
$\Phi(\boldsymbol{k},\boldsymbol{x}))^{\mathsf{T}}$ is constant and has rank
$1$. Trivially the Jacobian has everywhere rank $1$ and so
$\Phi(\boldsymbol{k},\boldsymbol{x})$ is not complete. The explicit solutions
of the ODE system are

$\displaystyle x_{1}(t)$ $\displaystyle=$ $\displaystyle(1-c_{0})t+c_{1},$
$\displaystyle x_{2}(t)$ $\displaystyle=$
$\displaystyle(c_{0}-1)t+c_{0}-c_{1}.$

Thus, all invariant curves are of the form $x_{1}+x_{2}=c_{0}$. We conclude
that there are no further first integrals and so the system has no complete
set of conservation laws. We can also note that the intersection of the steady
state variety defined by $x_{1}+x_{2}=1$ with a stoichiometric compatibility
class $\\{x_{1}+x_{2}=c_{0},x_{1}>0,x_{2}>0\\}$ is either empty or a line
segment.

In the absence of the parameters $\boldsymbol{k}$ Algorithm 4
$\operatorname{IsComplete}$ computes the parametric rank
$R=\\{(\operatorname{true},1)\\}$. However, we have $n=2$, which yields
$R_{2}=\emptyset$ and $\varrho_{2}=\operatorname{false}$. This gives us

$\gamma=\forall\,\boldsymbol{k}\,\forall\,\boldsymbol{x}\,(\operatorname{true}\,\mathrel{\land}\,\boldsymbol{x}>0\,\land\,x_{1}+x_{2}-1=0\longrightarrow\operatorname{false}),$

which quantifier elimination identifies to be “$\operatorname{false}$”, and
the algorithm returns “no”. Notice that $\gamma$ is equivalent to

$\forall\,\boldsymbol{x}\,(\boldsymbol{x}>0\longrightarrow x_{1}+x_{2}-1\neq
0),$

which illustrates that with deficient rank, completeness can only hold
formally when there is no steady state in the positive orthant.

##  5 Computing Linear Conservation Laws

By definition a vector $\boldsymbol{c}=(c_{1},\ldots,c_{n})\in\mathbb{R}^{n}$
defines an exact linear conservation law

$\Phi(\boldsymbol{k},\boldsymbol{x})=\sum_{i=1}^{n}c_{i}x_{i}$

of the system in (1), if

$\sum_{i=1}^{n}c_{i}f_{i}(\boldsymbol{k},\boldsymbol{x})=\sum_{j=1}^{r}(\boldsymbol{c}\boldsymbol{S})_{j}k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}=0$

for all $\boldsymbol{k}$, $\boldsymbol{x}$. Exact linear conservation laws of
the system in (1) can be obtained from the matrix $\boldsymbol{S}$ as the
following theorem shows.

######  Theorem 19.

Assume that all monomial reaction rates with the same multi-index have
different rate constants, i.e.
$\boldsymbol{\alpha}_{j}=\boldsymbol{\alpha}_{j^{\prime}}$ with $j\neq
j^{\prime}$ implies $k_{j}\neq k_{j^{\prime}}$. This condition is always
fulfilled in our setting (see the Remark 2). Let $r_{\boldsymbol{S}}$ be the
rank of the matrix $\boldsymbol{S}$. Then there is a full rank matrix
$\boldsymbol{C}$ with $n-r_{\boldsymbol{S}}$ rows such that
$\boldsymbol{C}\boldsymbol{S}=0$. Furthermore, $\boldsymbol{C}$ can be chosen
such that its rows form a set of independent, irreducible exact linear
conservation laws unconditional on $\boldsymbol{k}$. Finally, all
$\boldsymbol{c}=(c_{1},\ldots,c_{n})$ defining unconditional, exact linear
conservation laws are linear combinations of the rows of $\boldsymbol{C}$.

###### Proof.

As a direct consequence of the rank-nullity theorem, there is a matrix
$\boldsymbol{C}$ with $n-r_{\boldsymbol{S}}$ independent rows such that
$\boldsymbol{C}\boldsymbol{S}=0$. If $\boldsymbol{c}=(c_{1},\ldots,c_{n})$
defines an exact, unconditional linear conservation law. It follows that

$\sum_{j=1}^{r}(\boldsymbol{c}\boldsymbol{S})_{j}k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}=0$

for all $\boldsymbol{k}$, $\boldsymbol{x}$. If $i\neq j$ one has
$\boldsymbol{\alpha}_{i}\neq\boldsymbol{\alpha}_{j}$ or $k_{i}\neq k_{j}$. It
follows that $k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}$ are linearly
independent in $\mathbb{R}[\boldsymbol{k},\boldsymbol{x}]$ and
$\boldsymbol{c}\boldsymbol{S}=0$. Therefore, $\boldsymbol{c}$ is a linear
combination of rows of $\boldsymbol{C}$. Those rows of $\boldsymbol{C}$ that
are not irreducible can be decomposed as a sum of irreducible linear
conservation laws. The resulting conservation laws form a generating set for
the vector space of linear conservation laws and therefore this set contains a
basis. The $n-r_{\boldsymbol{S}}$ elements of this basis can be chosen as the
irreducible rows of $\boldsymbol{C}$. ∎

######  Remark 20.

Since $S_{ij}\in\mathbb{Z}$, the coefficients $c_{i}$ can be chosen to be
integers. For conservative CRNs linear conservation laws can be chosen semi-
positive, that is all $c_{i}\geq 0$ (see [57]). Semi-positive linear
conservation laws are important, because their existence implies that
concentrations of some species are bounded at all times. Also, semi-positive
conservation laws allow easy interpretation in terms of pools of chemical
species. An algorithm to compute irreducible, semi-positive linear
conservation laws can be found in [57]. We must emphasize that it is not
always possible to transform complete systems of linear conservation laws into
semi-positive linear conservation laws.

######  Remark 21.

In Example 18 we have seen that $n-r_{\boldsymbol{S}}$ independent linear
conservation laws may not form a complete system. The system in this example
is as in (1), although it is not a mass action network.

By applying the Algorithm 1 to the truncated system
$\boldsymbol{F}^{(1)}(\boldsymbol{k},\boldsymbol{x})$ we obtain a matrix
$\boldsymbol{S}^{(1)}=(S_{ij}^{(1)})\in\mathbb{Z}^{n\times r}$ such that

$f_{i}^{(1)}(\boldsymbol{k},\boldsymbol{x})=\sum_{j=1}^{r}S^{(1)}_{ij}k_{j}\boldsymbol{x}^{\boldsymbol{\alpha}_{j}}.$

The matrix $\boldsymbol{S}^{(1)}$ is called truncated stoichiometric matrix.
Since the truncated system contains less monomials than the full system,
$\boldsymbol{S}^{(1)}$ is obtained from $\boldsymbol{S}$ by setting the
respective entries in $\boldsymbol{S}$ to zero.

######  Remark 22.

Approximate conservation laws are computed in the same way as the exact ones
after replacing the matrix $\boldsymbol{S}$ by the truncated stoichiometric
matrix $\boldsymbol{S}^{(1)}$.

######  Theorem 23.

Consider the Jacobian matrix
$\boldsymbol{J}(\boldsymbol{k},\boldsymbol{x})=D_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})$
and the matrix $\boldsymbol{C}$ introduced in Theorem 19. The rows of
$\boldsymbol{C}$ provide a complete set of linear conservation laws, if for
all $\boldsymbol{k}\in\mathbb{R}_{+}^{r}$,
$\boldsymbol{x}\in\mathbb{R}_{+}^{n}$ with
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=\mathbf{0}$ the following
conditions hold:

  1. (i)

$rk(\boldsymbol{J}(\boldsymbol{k},\boldsymbol{x}))=rk(\boldsymbol{S})=r_{\boldsymbol{S}}$,

  2. (ii)

no row of the product $J(\boldsymbol{k},\boldsymbol{x})\boldsymbol{S}$ is
zero.

###### Proof.

We check that the assumptions imply the conditions in Definition 11. From (i)
it follows that $\boldsymbol{J}(\boldsymbol{k},\boldsymbol{x})$ has
$r_{\boldsymbol{S}}$ independent rows. Moreover, from (ii) we obtain that
these rows are not in the left kernel of $\boldsymbol{S}$, which is spanned by
the $n-r_{\boldsymbol{S}}$ rows of $\boldsymbol{C}$. Thus the matrix
$(\boldsymbol{J}(\boldsymbol{k},\boldsymbol{x}),\boldsymbol{C})^{\mathsf{T}}$
has rank $n$. ∎

The CRN in Example 18 has no complete set of linear conservation laws. This
model fails to satisfy condition (ii) of the theorem, but it fulfills (i). The
network is not a mass action CRN, so one may ask whether mass action CRNs
automatically satisfy both (i) and (ii). The answer to this question is no, as
is shown by the following example.

######  Example 24.

The CRN

$\dot{x}_{1}=-k_{1}x_{1}-k_{2}x_{1}-2k_{3}x_{1}^{2},\quad\dot{x}_{2}=2k_{3}x_{1}^{2}+2k_{2}x_{1},\quad\dot{x}_{3}=k_{3}x_{1}^{2}+k_{1}x_{1}$

is a mass action network described by the reactions

$A_{1}\xrightarrow{k_{1}}A_{3},\quad A_{1}\xrightarrow{k_{2}}2A_{2},\quad
2A_{1}\xrightarrow{k_{3}}2A_{2}+A_{3}.$

Its stoichiometric matrix is

$\boldsymbol{S}=\begin{pmatrix}-1&-1&-2\\\ 0&2&2\\\ 1&0&1\end{pmatrix}$

and one easily checks that it has rank two. There is a single linear
conservation law, namely
$\phi(\boldsymbol{k},\boldsymbol{x})=2x_{1}+x_{2}+2x_{3}$. The Jacobian matrix
of the system is

$J(\boldsymbol{k},\boldsymbol{x})=\begin{pmatrix}-k_{1}-k_{2}-4k_{3}x_{1}&0&0\\\
2k_{2}+4k_{3}x_{1}&0&0\\\ k_{1}+2k_{3}x_{1}&0&0\end{pmatrix},$

which clearly has rank one. Therefore, this model does not satisfy condition
(i) of the Theorem 23, whereas (ii) is fulfilled. Furthermore, the equations
$\boldsymbol{F}(\boldsymbol{k},\boldsymbol{x})=0$,
$\phi(\boldsymbol{k},\boldsymbol{x})=c_{0}$, which define the intersection of
$\mathcal{S}_{\boldsymbol{k}}$ with the stoichiometric compatibility class,
have degenerate solutions $x_{1}=0,\,x_{2}=c_{0}-2x_{3}$ and so the
conservation law is not complete.

In [26], sufficient conditions were found for mass action CRNs to be
“injective”, meaning to have a unique steady state in any stoichiometric
compatibility class. The same conditions also imply completeness of any
$n-r_{\boldsymbol{S}}$ independent, linear conservation laws (also see [26]).

##  6 Computing Nonlinear Conservation Laws

Nonlinear conservation laws, i.e. first integrals, of polynomial systems of
ODEs have been approached by means of a variety of methods. For example,
elementary conservation laws, i.e. expressions composed of polynomials,
logarithms and exponentials, can be obtained by Darboux polynomials [52, 13].
The bottleneck of this method is the computation of all Darboux polynomials of
the ODE system. This problem has been approached in computational algebra by
using Gröbner bases [48]. For particular CRN models, roughly Lotka-Volterra
models, Darboux polynomials and conservation laws of specific kind can be
obtained with the methods described in [47].

In this section we restrict ourselves to polynomial conservation laws. We show
that polynomial conservation laws are obtained from syzygies, that are easier
to compute than Darboux polynomials. Our methods cover CRNs more generally
than those studied in [47]. However, we do not cope with parametric
exponential conservation laws existing in Lotka-Volterra models.

We also discuss monomial and rational monomial polynomials that can be
effectively computed without syzygies or Darboux polynomials. We show that the
computation of monomial conservation laws can be reduced to the computation of
linear conservation laws. Finally we use the databases Biomodels [12] and
ODEbase [45], to study the scalability of our methods and the propensity of
CRNs to have linear and nonlinear conservation laws.

###  6.1 Computing Monomial Conservation Laws

The following proposition establishes a link between linear and monomial
conservation laws.

######  Proposition 25.

Let $E$ be a system of ODEs given by

$\dot{x_{1}}=f_{1}(\boldsymbol{k},\boldsymbol{x}),\,\dots,\,\dot{x_{n}}=f_{n}(\boldsymbol{k},\boldsymbol{x})$

and let $E^{\prime}$ be the system of ODEs given by

$\dot{x_{1}}=\frac{f_{1}(\boldsymbol{k},\boldsymbol{x})}{x_{1}},\,\dots,\,\dot{x_{n}}=\frac{f_{n}(\boldsymbol{k},\boldsymbol{x})}{x_{n}}.$

Then $E$ admits a monomial conservation law if and only if $E^{\prime}$ admits
a linear conservation law with integer coefficients. Moreover, if the linear
conservation law for $E^{\prime}$ is $m_{1}x_{1}+\dots+m_{n}x_{x}$ then the
monomial conservation law for $E$ is $x_{1}^{m_{1}}\cdots x_{n}^{m_{n}}$. The
linear conservation laws of $E^{\prime}$ are irreducible if and only if the
corresponding monomial conservation laws of $E$ are irreducible.

###### Proof.

Let $M=x_{1}^{m_{1}}\cdots x_{n}^{m_{n}}$ be a monomial. Then the derivative
of $M$ is

$\dot{M}=\big{(}\sum_{i=1}^{n}m_{i}\frac{\dot{x_{i}}}{x_{i}}\big{)}M.$

The equivalences

$\dot{M}=0\Longleftrightarrow\sum_{i=1}^{n}m_{i}\frac{\dot{x_{i}}}{x_{i}}=0\Longleftrightarrow\sum_{i=1}^{n}m_{i}\frac{f_{i}(\boldsymbol{k},\boldsymbol{x})}{x_{i}}=0$

prove the first part of the proposition.

Let $\phi=m_{1}x_{1}+\dots+m_{n}x_{n}$ be the linear conservation law of
$E^{\prime}$ corresponding to the monomial conservation law $M$ of $E$. Assume
that $M$ is reducible, say $M=M_{1}M_{2}$ with $M_{k}=x_{1}^{m_{k1}}\cdots
x_{n}^{m_{kn}}$ and $\dot{M}_{k}=0$ for $k\in\\{1,2\\}$. Then by the first
part

$\dot{M}_{k}=0\Longleftrightarrow\sum_{i=1}^{n}m_{ki}\frac{\dot{x_{i}}}{x_{i}}=0$

for $k\in\\{1,2\\}$. Thus, $\phi=\phi_{1}+\phi_{2}$, where
$\phi_{k}=m_{k1}x_{1}+\dots+m_{kn}x_{n}$ are linear conservation laws of
$E^{\prime}$ for $k\in\\{1,2\\}$. This means that the linear conservation for
$E^{\prime}$ is also reducible. The reverse follows from the same reasoning. ∎

######  Remark 26.

Proposition 25 holds even when the r.h.s. of $E^{\prime}$ are sums of rational
monomials, which can be the case when $E$ is a system as in (1) and some of
the polynomials $f_{i}(\boldsymbol{k},\boldsymbol{x})$ are not a multiple of
$x_{i}$. In this case, a stoichiometric matrix $\boldsymbol{S}^{\prime}$ can
still be associated to $E^{\prime}$ using Algorithm 1 and the linear
conservation laws of $E^{\prime}$ can be computed from the matrix
$\boldsymbol{S}^{\prime}$. This procedure is summarized in Algorithm 6.

0: Polynomial ODE system $E$ written as in (1).

0: Integer coefficient matrix $\boldsymbol{C}$ whose rows define monomial
conservation laws $\boldsymbol{x}^{\boldsymbol{C}_{i}}$ for $E$.

1: Compute $E^{\prime}$ by dividing each ODE in $E$ by $x_{i}$.

2: $\boldsymbol{S}^{\prime}=\text{Smatrix}(E^{\prime})$.

3: Find $\boldsymbol{C}$ such that $\boldsymbol{C}\boldsymbol{S}^{\prime}=0$
and $\boldsymbol{C}$ has $n-r_{\boldsymbol{S}}$ irreducible, independent rows,
where $r_{\boldsymbol{S}}=rk(\boldsymbol{S}^{\prime})$.

Algorithm 6 $\operatorname{StoichiometricMonomialConservation}$

######  Example 27.

The model

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle
x_{1}(x_{2}-x_{1})-\delta x_{1},$ $\displaystyle\dot{x}_{2}$ $\displaystyle=$
$\displaystyle x_{2}(x_{1}-x_{2})$

is a mass action network and is described by the reactions

$A_{1}+A_{2}\xrightarrow{1}2A_{2},\quad
A_{1}+A_{2}\xrightarrow{1}\emptyset,\quad
A_{2}+A_{2}\xrightarrow{1}\emptyset,\quad A_{1}\xrightarrow{\delta}\emptyset.$

The truncated system

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle
x_{1}(x_{2}-x_{1}),$ $\displaystyle\dot{x}_{2}$ $\displaystyle=$
$\displaystyle x_{2}(x_{1}-x_{2})$

has the steady state variety $x_{1}=x_{2}$ and Jacobian matrix is singular on
it. After truncation and elimination of the factors $x_{1}$, $x_{2}$ we obtain
the system

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle x_{2}-x_{1},$
$\displaystyle\dot{x}_{2}$ $\displaystyle=$ $\displaystyle x_{1}-x_{2},$

which has the linear irreducible conservation law
$\phi_{l}(x_{1},x_{2})=x_{1}+x_{2}$. We conclude that the model has the
monomial irreducible approximate conservation law
$\phi_{m}(x_{1},x_{2})=x_{1}x_{2}$. The monomial conservation law is complete,
since the $2\times 2$ minor

$\mathrm{det}(D_{\boldsymbol{x}}(x_{1}(x_{2}-x_{1}),x_{1}x_{2})^{\mathsf{T}})=-2x_{1}^{2}$

of the Jacobian matrix does not vanish for $x_{1}>0$. For
$c_{0}\in\mathbb{R}_{+}$ the intersection of the steady state variety with the
compatibility class defined by $\phi_{m}(x_{1},x_{2})=c_{0}$ and
$x_{1},x_{2}>0$ is the point $x_{1}=x_{2}=\sqrt{c_{0}}$.

######  Example 28.

A model by Volpert (see [47]) is described by the mass action network

$A_{1}+A_{2}\xrightarrow{1}2A_{2},\quad A_{2}+A_{3}\xrightarrow{2}2A_{3},\quad
A_{3}+A_{1}\xrightarrow{2}2A_{1},$

and the ODEs

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle
x_{1}(x_{3}-x_{2}),$ $\displaystyle\dot{x}_{2}$ $\displaystyle=$
$\displaystyle x_{2}(x_{1}-x_{3}),$ $\displaystyle\dot{x}_{3}$
$\displaystyle=$ $\displaystyle x_{3}(x_{2}-x_{1})$

have the steady state variety $x_{1}=x_{2}=x_{3}$. Moreover, the
stoichiometric matrix

$\boldsymbol{S}=\begin{pmatrix}-1&0&1\\\ 1&-1&0\\\ 0&1&-1\end{pmatrix}$

has rank $r_{\boldsymbol{S}}=2$ and its left kernel is generated by
$\boldsymbol{c}=(1,1,1)$. The stoichiometric matrix of the system obtained
after elimination of the factors $x_{1}$, $x_{2}$, $x_{3}$ in the truncated
system is

$\boldsymbol{S}^{\prime(1)}=\begin{pmatrix}0&1&-1\\\ -1&0&1\\\
1&-1&0\end{pmatrix}$

and clearly has rank $r_{\boldsymbol{S}}=2$. Moreover, its left kernel is also
generated by $\boldsymbol{c}=(1,1,1)$. Thus the model has two approximate
conservation laws, one linear $\phi_{l}(\boldsymbol{x})=x_{1}+x_{2}+x_{3}$ and
one monomial $\phi_{m}(\boldsymbol{x})=x_{1}x_{2}x_{3}$. Both conservation
laws are complete. Indeed, the Jacobian
$D_{\boldsymbol{x}}(F_{1},F_{2},F_{3},\phi_{m})^{\mathsf{T}}$ has a non-
vanishing $3\times 3$ minor

$\mathrm{det}(D_{\boldsymbol{x}}(x_{1}(x_{3}-x_{2}),x_{2}(x_{1}-x_{3}),x_{1}x_{2}x_{3})^{\mathsf{T}})=x_{1}^{2}x_{2}x_{3}+x_{1}x_{2}^{2}x_{3}+x_{1}x_{2}x_{3}^{2}.$

Moreover, the Jacobian
$D_{\boldsymbol{x}}(F_{1},F_{2},F_{3},\phi_{l})^{\mathsf{T}}$ has a $3\times
3$ minor

$\displaystyle M$ $\displaystyle=$
$\displaystyle\mathrm{det}(D_{\boldsymbol{x}}(x_{1}(x_{3}-x_{2}),x_{2}(x_{1}-x_{3}),x_{1}+x_{2}+x_{3})^{\mathsf{T}})$
$\displaystyle=$
$\displaystyle-(x_{1}-x_{2})^{2}+2x_{1}x_{3}+2x_{2}x_{3}-x_{3}^{2}$

On $\mathcal{S}_{\boldsymbol{k}}\cap\mathbb{R}_{+}^{3}$ we have
$x_{1}=x_{2}=x_{3}$ and so $M=3x_{1}^{2}$ does not vanish there.

However, the two conservation laws are not independent on
$\mathcal{S}_{\boldsymbol{k}}\cap\mathbb{R}_{+}^{3}$. Indeed, for
$x_{1}=x_{2}=x_{3}$ we have

$D_{\boldsymbol{x}}(\phi_{l},\phi_{m})^{\mathsf{T}}=\begin{pmatrix}1&1&1\\\
x_{1}^{2}&x_{1}^{2}&x_{1}^{2}\end{pmatrix}.$

We want to point out that in this example we have used a special choice of the
rate constants. For rational rate constants, there are rational monomial
conservation laws with parameter dependent exponents [47]. Parameter dependent
monomial conservation laws can be computed with methods presented in [30].

###  6.2 Computing Polynomial Conservation Laws from Syzygies

In this subsection, we present an algorithm for computing parametric
polynomial conservation laws. Our algorithm uses syzygies, a well-known
concept in computational commutative algebra and algebraic geometry, as well
as comprehensive Gröbner systems, in order to compute parametric conservation
laws. We also explain another algorithm for computing only unconditional
conservation laws that avoids Gröbner basis computation. This algorithm is
based on an ansatz and was previously used for computing polynomial first
integrals [32]. The algorithm makes an ansatz for an unconditional polynomial
conservation law with additional indeterminates as coefficients. One will then
obtain a linear system in this indeterminates from which the coefficients of
the conservation laws can be computed.

######  Definition 29.

Let $R=\mathbb{Q}[\boldsymbol{x}]$ be the polynomial ring in indeterminates
$\boldsymbol{x}=(x_{1},\ldots,x_{n})$ over the field of rational numbers
$\mathbb{Q}$ and let $I=\langle f_{1},\ldots,f_{m}\rangle\subseteq R$ be an
ideal in $R$. A syzygy of $I$ with respect to the basis
$\\{f_{1},\ldots,f_{m}\\}$ is an element
$\boldsymbol{g}=(g_{1},\ldots,g_{m})\in R^{m}$ such that

$g_{1}f_{1}+\ldots+g_{m}f_{m}=0.$

A syzygy of the form

$\boldsymbol{g}=(0,\dots,0,f_{j},0,\dots,0,-f_{i},0,\dots,0),$

where $f_{j}$ is the $i-$th and $-f_{i}$ is the $j-$th entry of
$\boldsymbol{g}$, is called a trivial syzygy.

######  Remark 30.

In this section, for computational purposes, we consider polynomials over
rational numbers, that can be evaluated in $\mathbb{R}$.

######  Remark 31.

In the definition of a syzygy, some authors require $\\{f_{1},\ldots,f_{m}\\}$
to be a minimal generating set of the ideal $I$ [23]. However, for the purpose
of this paper, we do not require such a condition on the generators. Moreover,
it is well known that syzygies of an ideal form a submodule of $R^{m}$ as an
$R-$module.

######  Example 32.

Consider the ideal $I=\langle x,y\rangle\subseteq\mathbb{Q}[x,y]$. Since

$(-y)x+(x)y=0,$

the vector $\boldsymbol{g}=(-y,x)\in\mathbb{Q}[x,y]^{2}$ is a trivial syzygy
of the ideal $I$.

Syzygies have been extensively studied in commutative algebra and algebraic
geometry. Buchberger’s algorithm for computing Gröbner bases [11, 7] leads to
an algorithm which computes a basis of a syzygy module (see Schreyer’s results
in [6, 56]). One can also obtain a basis for the syzygy module during Gröbner
bases computations via signature-based algorithms [28, 69, 34]. Many computer
algebra systems such as Singular, CoCoA, etc. include implementations for
computing syzygies. For more on syzygies and in particular on algorithms for
their computation we refer to the books of Weispfenning [5] and Cox, Little &
O’Shea [14, Chapter 5.3].

Let us consider the system (1), i.e. the system of ODEs of the form

$\dot{x}_{1}=f_{1}(\boldsymbol{k},\boldsymbol{x}),\,\ldots,\,\dot{x}_{n}=f_{n}(\boldsymbol{k},\boldsymbol{x})\in\mathbb{Z}[\boldsymbol{k},\boldsymbol{x}],$

and assume that
$\phi(\boldsymbol{k},\boldsymbol{x})\in\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$
is a parametric exact polynomial conservation law with semi-algebraical set
$V_{\boldsymbol{k}}$ for the parameters $\boldsymbol{k}$ for (1). By
Definition 5 we have

$\sum_{i=1}^{n}\frac{\partial\phi}{\partial
x_{i}}(\boldsymbol{k},\boldsymbol{x})f_{i}(\boldsymbol{k},\boldsymbol{x})=0,\quad\forall\boldsymbol{k}\in
V_{\boldsymbol{k}},\ \forall\boldsymbol{x}\in\mathbb{R}^{n},$

which implies that

$\left(\frac{\partial\phi}{\partial
x_{1}}(\boldsymbol{k},\boldsymbol{x}),\,\dots,\,\frac{\partial\phi}{\partial
x_{n}}(\boldsymbol{k},\boldsymbol{x})\right)\in\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]^{n}$

is a syzygy of the steady state ideal $I=\langle f_{1},\dots,f_{n}\rangle$ for
all $\boldsymbol{k}\in V_{\boldsymbol{k}}$. This means that the gradient

$\gradient\phi(x)=\left(\frac{\partial\phi}{\partial
x_{1}}(\boldsymbol{k},\boldsymbol{x}),\,\dots,\,\frac{\partial\phi}{\partial
x_{n}}(\boldsymbol{k},\boldsymbol{x})\right)$

of a parametric polynomial conservation law
$\phi(\boldsymbol{k},\boldsymbol{x})$ is a syzygy of the steady state ideal
for all $\boldsymbol{k}\in V_{\boldsymbol{k}}$.

On the other hand, if we are given a syzygy $\boldsymbol{g}$ of the steady
state ideal, which is valid only for values of $\boldsymbol{k}$ coming from
the semi-algebraic set $V_{\boldsymbol{k}}$, such that $\boldsymbol{g}$ is the
gradient of a polynomial
$\phi(\boldsymbol{k},\boldsymbol{x})\in\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$,
then $\phi(\boldsymbol{k},\boldsymbol{x})$ is a parametric polynomial
conservation law. Roughly speaking, this implies that for a parametric syzygy
$\boldsymbol{g}=(g_{1},\dots,g_{n})$, one can obtain a parametric conservation
law by integrating $\boldsymbol{g}$. We are going to clarify this idea.

Consider the differential form $d\boldsymbol{g}=g_{1}dx_{1}+\dots+g_{n}dx_{n}$
and take a piecewise smooth curve connecting the origin with $\boldsymbol{x}$,
for instance

$\displaystyle\mathcal{C}=\\{(sx_{1},0,\dots,0)\mid s\in[0,1]\\}\cup$
$\displaystyle\\{(x_{1},sx_{2},0,\dots,0)\mid
s\in[0,1]\\}\cup\dots\cup\\{(x_{1},\dots,x_{n-1},sx_{n})\mid s\in[0,1]\\}.$

Integrating the differential form $d\boldsymbol{g}$ along the curve
$\mathcal{C}$ yields

$\int_{0}^{x_{1}}g_{1}(\cdot,0,...,0)dx_{1}+\int_{0}^{x_{2}}g_{2}(x_{1},\cdot,0,...,0)dx_{2}+\dots+\int_{0}^{x_{n}}g_{n}(x_{1},...,x_{n-1},\cdot)dx_{n}$

and as a result of this integration we obtain a polynomial
$\phi(\boldsymbol{k},\boldsymbol{x})\in\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$.
The rotational of the vector field generated by $\boldsymbol{g}$ defined as

$[\curl\boldsymbol{g}]_{ij}=\frac{\partial g_{i}}{\partial
x_{j}}-\frac{\partial g_{j}}{\partial x_{i}}$

is important in this procedure. Indeed, if $\curl\boldsymbol{g}\neq 0$, the
result of the integration depends on the path $\mathcal{C}$ and there is no
function $\phi(\boldsymbol{k},\boldsymbol{x})$ satisfying
$\gradient\phi(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{g}$ for all
$\boldsymbol{x}\in\mathbb{R}^{n}$. The function
$\phi(\boldsymbol{k},\boldsymbol{x})$ is well defined in a simply connected
domain $\mathcal{D}\subset\mathbb{R}^{n}$ if only if $\curl\boldsymbol{g}=0$
everywhere in $\mathcal{D}$, i.e. only when $\boldsymbol{g}$ defines a
conservative vector field in $\mathcal{D}$ [62]. If the irrotationality
condition, i.e. $\curl\boldsymbol{g}=0$, is satisfied, then one has
$\gradient\phi(\boldsymbol{k},\boldsymbol{x})=\boldsymbol{g}$ and
$\gradient\phi(\boldsymbol{k},\boldsymbol{x})\cdot
F(\boldsymbol{k},\boldsymbol{x})=0$ everywhere in $\mathcal{D}$.

Another method to obtain a parametric conservation law from a syzygy
$\boldsymbol{g}=(g_{1},\dots,g_{n})$ is symbolic integration of each
coordinate $g_{i}$. However, some considerations should be taken into account.
First of all, the integration of the polynomial in the $i-$th entry of a
syzygy vector is done with respect to the variable $x_{i}$. Thus

$\phi(\boldsymbol{k},\boldsymbol{x})=\int_{x_{i}}g_{i}dx_{i}+c_{i}$

with an integration constant
$c_{i}\in\mathbb{Q}[\boldsymbol{k},x_{1},\dots,\check{x}_{i},\dots,x_{n}]$,
that is $c_{i}$ is a polynomial of $\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$
not involving the indeterminate $x_{i}$. The polynomial $c_{i}$ can be
computed with the following heuristic. For some degree $d$ consider the
polynomial $c_{i,d}(\boldsymbol{a},\boldsymbol{k},\boldsymbol{x})$ consisting
of all monomials in $\boldsymbol{k}$, $\boldsymbol{x}$ of degree less or equal
to $d$ whose coefficients are new indeterminates
$\boldsymbol{a}=(a_{1},\dots,a_{w})$ with $w=\binom{n+r+d}{d}$. Collecting in
the equation

$\gradient(\int_{x_{i}}g_{i}dx_{i}+c_{i,d}(\boldsymbol{a},\boldsymbol{k},\boldsymbol{x}))\cdot
F(\boldsymbol{k},\boldsymbol{x})=0$

the terms in $\boldsymbol{k}$, $\boldsymbol{x}$ the coefficients define a
linear system in the indeterminates $\boldsymbol{a}$. If the system has a
solution $\boldsymbol{\overline{a}}$ in $\mathbb{R}^{w}$, then

$\int_{x_{i}}g_{i}dx_{i}+c_{i,d}(\boldsymbol{\overline{a}},\boldsymbol{k},\boldsymbol{x}))$

is a conservation law. Otherwise we need to increase the degree $d$ and repeat
the computation. We will explain below that the procedure will stop for some
degree $d$ and every $i$, if the irrotationality condition is satisfied.

As noticed above, once we obtained a syzygy, we have to check if the vector
field defined by it is conservative. It is a classical result that a vector
field is conservative over a simply connected set if and only if it is
irrotational [62]. Thus, if $\curl{\boldsymbol{g}}=0$ is not satisfied, the
linear system just described above is incompatible, as the equation
$\gradient\phi=\boldsymbol{g}$ has no solution. As an example, consider the
trivial syzygy $\boldsymbol{g}=(-y,x)$ for the ideal $\langle
x,y\rangle\subset\mathbb{Q}[x,y]$. As the vector field arising from this
particular syzygy is not irrotational, i.e., $\curl(-y,x)\neq 0$, one cannot
obtain a parametric polynomial conservation law from this syzygy. But there
exist steady state ideals with a trivial syzygy which is the gradient of a
parametric conservation law. For example, consider the ideal $\langle
x,-y\rangle\subset\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$ with the trivial
syzygy $(y,x)$. Then $\curl(y,x)=0$, which means that the trivial syzygy is
irrotational and so it is the gradient of a parametric conservation law, which
is $xy=c$, where $c$ is a constant.

In order to compute a parametric polynomial conservation law, one needs to
find a corresponding semi-algebraic set $V_{\boldsymbol{k}}$. The set
$V_{\boldsymbol{k}}$ determines the values for the parameters $\boldsymbol{k}$
for which the conservation law is valid. We obtain these semi-algebraic sets
from the branches of a comprehensive Gröbner system of the steady state ideal.
We introduce shortly some basics of comprehensive Gröbner systems.

Consider the polynomial ring $\mathbb{Q}[\boldsymbol{k}][\boldsymbol{x}]$ in
the indeterminates $\boldsymbol{x}$ over the coefficient ring
$\mathbb{Q}[\boldsymbol{k}]$. A specialisation of $\mathbb{Q}[\boldsymbol{k}]$
is a homomorphism $\sigma:\mathbb{Q}[\boldsymbol{k}]\rightarrow\mathbb{Q}$.
Every element $\boldsymbol{\alpha}\in\mathbb{Q}^{r}$ defines a specialisation
$\sigma_{\boldsymbol{\alpha}}$ by evaluating $\boldsymbol{k}$ at
$\boldsymbol{\alpha}$ and it can be canonically extended to a specialisation

$\sigma_{\boldsymbol{\alpha}}:\mathbb{Q}[\boldsymbol{k}][\boldsymbol{x}]\rightarrow\mathbb{Q}[\boldsymbol{x}].$

######  Definition 33.

Let $F$ and $G_{1},\dots,G_{l}$ be finite subsets of
$\mathbb{Q}[\boldsymbol{k}][\boldsymbol{x}]$. Further let $A_{1},\dots,A_{l}$
be semi-algebraic sets of $\mathbb{Q}^{r}$ and let $S\subseteq\mathbb{Q}^{r}$
such that $S\subset A_{1}\cup\dots\cup A_{l}$. A finite set

$\mathcal{G}=\\{(A_{1},G_{1}),\dots,(A_{l},G_{l})\\}$

is called a comprehensive Gröbner system (CGS) on $S$ for $F$ if for all
$i=1,\dots,l$ the set $\sigma_{\boldsymbol{\alpha}}(G_{i})$ is a Gröbner basis
for the ideal
$\langle\sigma_{\boldsymbol{\alpha}}(F)\rangle\subset\mathbb{Q}[\boldsymbol{x}]$
for every $\boldsymbol{\alpha}\in A_{i}$. Each $(A_{i},G_{i})$ is called a
branch of $\cal{G}$. The generic branch corresponds to the semi-algebraic set
that assumes all leading coefficients are non-zero and ignores the other
cases. If $S=\mathbb{Q}^{r}$, then $\cal{G}$ is called a comprehensive Gröbner
system for $F$.

Comprehensive Gröbner systems were introduced by Weispfenning in 1992 in [71].
Independent of each other Kapur introduced the same concept in [37], calling
it parametric Gröbner bases. Today in the literature there are many
publications dealing with the theory and algorithmics of CGSs, for example
[66, 49, 15, 43, 39, 38, 33] only mentioning a few. For a survey paper on this
topic see [44]. Moreover, several computer algebra systems provide
implementations of algorithms to compute comprehensive Gröbner systems, such
as Reduce [35, 36, 19, 63] and Singular [16].

Let $(A_{i},G_{i})$ be a branch of a CGS of a steady state ideal, i.e., the
ideal generated by $f_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,$
$f_{n}(\boldsymbol{k},\boldsymbol{x})$. One can use the algorithms described
in [14, Chapter 5.3] to compute a basis of the syzygy module of the ideal
generated by $G_{i}$. If
$\boldsymbol{g}\in\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$ is an element of
the syzygy module of the ideal generated by $G_{i}$, then for every
$\boldsymbol{\alpha}\in A_{i}$, the specialisation
$\sigma_{\boldsymbol{\alpha}}(\boldsymbol{g})$ will be a syzygy of the ideal
generated by $\sigma_{\boldsymbol{\alpha}}(G_{i})$. Following [14, Lemma 3.7],
there exist matrices $B^{1,i}$ and $B^{2,i}$ such that

$\boldsymbol{G_{i}}B^{1,i}=F\quad\mathrm{and}\quad
FB^{2,i}=\boldsymbol{G_{i}},$

where
$F=(f_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))$
and $\boldsymbol{G_{i}}$ is the matrix corresponding to $G_{i}$.

Then for every syzygy $\boldsymbol{g}$ of $G_{i}$, the product
$B^{2,i}\boldsymbol{g}=\boldsymbol{f}$ is a syzygy of $F$, and also
$B^{1,i}\boldsymbol{f}=\boldsymbol{g}$. Moreover we have

$\boldsymbol{G_{i}}\boldsymbol{g}=FB^{2,i}B^{1,i}\boldsymbol{f}=F\boldsymbol{f},$

since $F=\boldsymbol{G_{i}}B^{1,i}=FB^{2,i}B^{1,i}$. Applying
$\boldsymbol{\sigma}_{\boldsymbol{\alpha}}$ to the equation we obtain

$0=\boldsymbol{\sigma}_{\boldsymbol{\alpha}}(\boldsymbol{G_{i}})\boldsymbol{\sigma}_{\boldsymbol{\alpha}}(\boldsymbol{g})=\boldsymbol{\sigma}_{\boldsymbol{\alpha}}(F)\boldsymbol{\sigma}_{\boldsymbol{\alpha}}(\boldsymbol{f}),$

and therefore, $\boldsymbol{\sigma}_{\boldsymbol{\alpha}}(\boldsymbol{f})$ is
a syzygy of $\boldsymbol{\sigma}_{\boldsymbol{\alpha}}(F)$ for every
$\boldsymbol{\alpha}\in A_{i}$. Therefore, having computed a CGS of the steady
state ideal and then applying the conversion matrices described above, one can
obtain parametric syzygies for $F$ with the corresponding semi-algebraic sets.

Assume that the syzygy $\boldsymbol{f}$ generates a conservative vector field
and let $\phi(\boldsymbol{k},\boldsymbol{x})$ be the polynomial obtained by
integrating $\boldsymbol{f}$. Then $\phi(\boldsymbol{k},\boldsymbol{x})$ is a
parametric conservation law for all $\boldsymbol{k}\in
V_{\boldsymbol{k}}=A_{i}\cap\mathbb{Q}_{+}^{r}$. The equations and
inequalities describing $V_{\boldsymbol{k}}$ are the equations and
inequalities describing $A_{i}$ along with $k_{1}>0,\dots,k_{r}>0$, which
represent the positive orthant. By applying the cylindrical algebraic
decomposition(CAD), implemented in e.g., QEPCAD [9], to the latter equations
and inequations one may obtain another (though not necessarily better)
representation of $V_{\boldsymbol{k}}$. The program SLFQ (Simplifying Large
Formulas with QEPCAD) [10] can be used once QEPCAD is applied to obtain such a
simplified version. Our discussions is summarised in Algorithm 7.

0: A polynomial ODE system
$F(\boldsymbol{k},\boldsymbol{x})=(f_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))$
whose r.h.s. $f_{i}(\boldsymbol{k},\boldsymbol{x})$ are polynomials in
$\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$.

0: A set of pairs
$\mathcal{S}=\\{(V_{\boldsymbol{k},1},\mathcal{C}_{1}),\dots,(V_{\boldsymbol{k},l},\mathcal{C}_{l})\\},$
where $V_{\boldsymbol{k},i}\subseteq\mathbb{Q}_{+}^{r}$ is a semi-algebraic
set and $\mathcal{C}_{i}$ is a set of parametric polynomial conservation laws
valid for all $\boldsymbol{k}\in V_{\boldsymbol{k},i}$.

1: $\mathcal{S}:=\emptyset$

2:
$I:=\operatorname{Ideal}(f_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))$
in $\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$

3: Compute a CGS $\mathcal{G}:=\\{(A_{1},G_{1}),\dots,(A_{l},G_{l})\\}$ of $I$

4: for $i=1,\dots,l$ do

5: $\mathcal{C}_{i}:=\emptyset$

6: Compute a basis $\operatorname{Syz}(G_{i})$ of the syzygy module of $G_{i}$

7: Convert $\operatorname{Syz}(G_{i})$ into a basis $SF$ of the parametric
syzygies of $F$

8: for
$\boldsymbol{f}(\boldsymbol{k},\boldsymbol{x})=(f_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))\in
SF$ do

9: if $\curl{\boldsymbol{f}}=0$ then

10:
$\phi(\boldsymbol{k},\boldsymbol{x}):=\int_{0}^{x_{1}}f_{1}dx_{1}+\int_{0}^{x_{2}}f_{2}dx_{2}+\dots+\int_{0}^{x_{n}}f_{n}dx_{n}$

11:
$\mathcal{C}_{i}:=\mathcal{C}_{i}\cup\\{\phi(\boldsymbol{k},\boldsymbol{x})\\}$

12: end if

13: end for

14: Compute $V_{\boldsymbol{k},i}:=A_{i}\cap\mathbb{Q}^{r}_{+}$

15: $\mathcal{S}:=\mathcal{S}\cup\\{(V_{\boldsymbol{k},i},\mathcal{C}_{i})\\}$

16: end for

17: return $\mathcal{S}$.

Algorithm 7 $\operatorname{ParametricPolynomialConservationLawsViaSyzygies}$

######  Remark 34.

By Definition 11 the set of $s$ parametric polynomial conservation laws
$\\{\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,\phi_{s}(\boldsymbol{k},\boldsymbol{x})\\}$
is independent if their Jacobian matrix

$\left(\begin{matrix}\frac{\partial\phi_{1}(\boldsymbol{k},\boldsymbol{x})}{\partial
x_{1}}&\dots&\frac{\partial\phi_{1}(\boldsymbol{k},\boldsymbol{x})}{\partial
x_{n}}\\\ \vdots&&\vdots\\\
\frac{\partial\phi_{s}(\boldsymbol{k},\boldsymbol{x})}{\partial
x_{1}}&\dots&\frac{\partial\phi_{s}(\boldsymbol{k},\boldsymbol{x})}{\partial
x_{n}}\end{matrix}\right)$

has rank $s$. So if $\boldsymbol{g}_{1},\dots,\boldsymbol{g}_{s}$ are syzygies
behind
$\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,\phi_{s}(\boldsymbol{k},\boldsymbol{x})$,
then independence of
$\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,\phi_{s}(\boldsymbol{k},\boldsymbol{x})$
implies linear independence of $\boldsymbol{g}_{1},\dots,\boldsymbol{g}_{s}$.
Conversely, if linearly independent syzygies
$\boldsymbol{g}_{1},\dots,\boldsymbol{g}_{s}$ are the gradients of parametric
conservation laws
$\phi_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,\phi_{s}(\boldsymbol{k},\boldsymbol{x})$,
then these conservation laws form an independent set. Since Algorithm 7 uses
for each branch a basis $\operatorname{Syz}(G_{i})$ of the syzygy module of
$G_{i}$, each pair of its output $(V_{\boldsymbol{k},i},\mathcal{C}_{i})$
defines a set of independent parametric polynomial conservation laws
$\mathcal{C}_{i}$ valid for all $\boldsymbol{k}\in V_{\boldsymbol{k},i}$.
Moreover, if there exists a complete set of parametric polynomial conservation
laws, then the Algorithm 7 returns such a set.

######  Remark 35.

Let $(A_{1},G_{1}),\dots,(A_{l},G_{l})$ be a CGS of the steady state ideal, as
in line 3 of Algorithm 7. By Definition 5 parametric conservation laws depend
on the rate constants $\boldsymbol{k}$, whereas unconditional conservation
laws are independent of $\boldsymbol{k}$. If $\phi$ is an unconditional
conservation law, then it is a parametric conservation law as well, hence,
$\phi$ must be obtained from a syzygy that does not contain variables
$\boldsymbol{k}$. This means that such a syzygy is a syzygy of every $G_{i}$.
Therefore, in order to compute unconditional conservation laws, we need to
find those syzygies which only depend on the variables $\boldsymbol{x}$ and
not on $\boldsymbol{k}$, for every branch of a CGS. In order to do so, we use
the elimination property of Gröbner bases. Fix a block term order on the
variables $\\{\boldsymbol{k},\boldsymbol{x}\\}$ with
$\boldsymbol{k}>\boldsymbol{x}$. let $\mathcal{M}$ be the module of syzygies
of the steady state ideal, and let $\mathcal{M_{\boldsymbol{x}}}$ be the
submodule of $\mathcal{M}$ that contains those elements that only include
variables $\boldsymbol{x}$. If $\mathcal{G_{M}}$ is a Gröbner basis for
$\mathcal{M}$ with respect to the block order, then
$G_{\mathcal{M}}\cap\mathbb{Q}[\boldsymbol{x}]^{m}$ is a Gröbner basis for
$\mathcal{M_{\boldsymbol{x}}}$, see [41, Thm 3.4.5]. This means that in order
to compute unconditional conservation laws, one can modify Algorithm 7 by
adding a computation of $G_{\mathcal{M}}\cap\mathbb{Q}[\boldsymbol{x}]^{m}$
after line 6.

######  Remark 36.

In the context of CRNs, some authors work with the ring
$\mathbb{Q}(\boldsymbol{k})[\boldsymbol{x}]$ without specialising
$\boldsymbol{k}$, while in algorithm 7, we work with the ring
$\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$, considering case distinctions
done for CGS computation, as $\boldsymbol{k}$ can be specialised in
$\mathbb{Q}_{+}^{r}$. One can see that a Gröbner bases of the steady state
ideal, as an ideal in $\mathbb{Q}(\boldsymbol{k})[\boldsymbol{x}]$,
corresponds to the Gröbner basis in the generic branch of a comprehensive
Gröbner system in the ring $\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$, with
the corresponding semi-algebraic set being the set generated by the leading
coefficients not equal to zero, saturated with respect to the product of the
leading terms of $F$. This fact about the generic branch has been considered
in [50], and a new algorithm for computing CGS has been presented based on
this. Among computer algebra systems, Singular has an implementation for
computing Gröbner bases over $\mathbb{Q}(\boldsymbol{k})[\boldsymbol{x}]$.

We implemented Algorithm 7 in Singular version 4.3.1, in a library called
polconslaw.lib. Our implementation uses grobcov library for computing
comprehensive Gröbner systems. We also implemented the approach in Remark 36
for computing polynomial conservation laws corresponding to the generic branch
of a CGS in a Singular library, called genpolconslaw.lib, which avoids step 3
in Algorithm 7 (for computing a comprehensive Gröbner system), and just
computes the conservation laws that correspond to the generic branch.
Therefore, this approach is much faster than computing the whole comprehensive
Gröbner system and picking up the generic branch. We have benchmarked the two
libraries using Biomodels [12], a comprehensive database of real life
biochemical CRN models. The polynomial CRN biomodels were reformatted as input
data for symbolic computation software and stored in the database ODEbase
[45].

SI Table 1 (see Supplementary Information) shows the output of computations in
generic branch obtained using our Singular library genpolconslaw.lib. SI Table
2 (see Supplementary Information) shows the result of parametric polynomial
conservation law computations using our Singular library polconslaw.lib.
Interestingly, the generic branch contains mostly unconditional linear
conservation laws, whereas the non-generic branches contain mostly parametric
polynomial, nonlinear conservation laws. The large propensity of real life
biochemical models for unconditional linear conservation laws may be explained
by Eigen’s idea that the hypercycle is a general organization principle in
biology [22]. A consequence of this principle is that biochemical CRNs can be
decomposed into cycles in which each species is a product of at least one
reaction and reactant of at least another one. The total amounts of different
moieties are unconditionally conserved in cycles.

Figure 1 shows that the number of syzygies in the generic branch scales
quadratically with the number $n$ of variables of the model, whereas the total
number of conservation laws (linear and polynomial) is sublinear in $n$.

We conclude this section with a series of examples clarifying our
computational methods.

Figure 1: Computation of polynomial conservation laws arising from the generic
branch of comprehensive Gröbner system for CRN models from the database
Biomodels. $n$ is the number of variables and ODEs in the CRN model. The
number of syzygies (blue dots) scales like $n^{2}$, which is mainly due to the
$n(n-1)/2$ trivial syzygies. The number of conservation laws (red dots)
depends sublinearly on $n$. The maximum running time for this computation was
3600 seconds (on a single core of an Intel i5-9400, 2.9 GHz processor).

######  Example 37.

We apply Algorithm 7 to the system

$\displaystyle\dot{x_{1}}=$ $\displaystyle-k_{1}x_{1},$
$\displaystyle\dot{x_{2}}=$ $\displaystyle-k_{2}x_{2},$
$\displaystyle\dot{x_{3}}=$ $\displaystyle(k_{1}+k_{2})x_{1}x_{2}.$

In line 6 of Algorithm 7 applied to the generic branch of a comprehensive
Gröbner basis, which is given by the semi-algebraic set defined by $k_{1}\neq
0$ and $k_{2}\neq 0$, the following basis for the syzygy module is computed.

$\operatorname{Syz}(F)=\\{(x_{2},x_{1},1),(-k_{2}x_{2},k_{1}x_{1},0)\\}.$

The first element of the syzygy basis is conservative and integrates to the
unconditional conservation law $\phi_{1}=x_{1}x_{2}+x_{3}$. The second element
is not conservative, because $\curl(-k_{2}x_{2},k_{1}x_{1},0)$ has the nonzero
entry

$\frac{\partial}{\partial x_{2}}(-k_{2}x_{2})-\frac{\partial}{\partial
x_{1}}(k_{1}x_{1})=-(k_{1}+k_{2})\neq 0.$

Other branches of the comprehensive Gröbner basis are not interesting, since
they lead to non-conservative syzygies.

######  Example 38.

For the system

$\displaystyle\dot{x}_{1}=$ $\displaystyle\ x_{2},$
$\displaystyle\dot{x}_{2}=$ $\displaystyle\ x_{1},$
$\displaystyle\dot{x}_{3}=$ $\displaystyle\ x_{1}x_{2},$

one can compute the following basis for the syzygy module:

$\\{(0,x_{2},-1),(x_{1},-x_{2},0)\\}.$

Both elements are conservative and integrate to the independent unconditional
conservation laws

$\frac{1}{2}x_{2}^{2}-x_{3}\quad\mathrm{and}\quad\frac{1}{2}x_{1}^{2}-\frac{1}{2}x_{2}^{2}.$

The two conservation laws are polynomial and one can easily check that they
form a complete set. Also, since the generating set of the steady state ideal
contains only polynomials depending on the variables $\boldsymbol{x}$ and not
on the variables $\boldsymbol{k}$, it is clear that the obtained conservation
laws are unconditional.

######  Example 39 (BIOMD0000000629, [46]).

The ODEs corresponding to Biomodel 629 in the ODEbase repository [46] are the
following:

$\displaystyle\dot{x_{1}}=$ $\displaystyle-k_{2}x_{1}x_{3}+k_{3}x_{2},$
$\displaystyle\dot{x_{2}}=$ $\displaystyle
k_{2}x_{1}x_{3}-k_{3}x_{2}-k_{4}x_{2}x_{4}+k_{5}x_{5},$
$\displaystyle\dot{x_{3}}=$ $\displaystyle-k_{2}x_{1}x_{3}+k_{3}x_{2},$
$\displaystyle\dot{x_{4}}=$ $\displaystyle-k_{4}x_{2}x_{4}+k_{5}x_{5},$
$\displaystyle\dot{x_{5}}=$ $\displaystyle k_{4}x_{2}x_{4}-k_{5}x_{5}.$

A comprehensive Gröbner system for the steady state ideal consists of 10
branches. We consider first the generic branch and after that one non generic
branch. One easily checks that the generic branch of the comprehensive Gröbner
basis consists of $\\{k_{2}\neq 0\ ,k_{4}\neq 0\\}$ and the set formed by the
five right hand sides of the system. Following the Remark 36, we carry out
computations over $\mathbb{Q}(\boldsymbol{k})[\boldsymbol{x}]$ in order to
find a basis of the syzygy module. We find

$\displaystyle\\{\boldsymbol{g}_{1}=(-1,0,1,0,0),\
\boldsymbol{g}_{2}=(0,-1,-1,1,0),\ \boldsymbol{g}_{3}=(0,0,0,1,1),$
$\displaystyle\boldsymbol{g}_{4}=(k_{2}x_{1}x_{3}-k_{4}x_{2}x_{4}-k_{3}x_{2}+k_{5}x_{5},k_{2}x_{1}x_{3}-k_{3}x_{2},0,0,0)\\}.$

Integrating the elements of the syzygy basis leads to the linear independent
unconditional conservation laws

$\phi_{1}(\boldsymbol{k},\boldsymbol{x})=-x_{1}+x_{3},\
\phi_{2}(\boldsymbol{k},\boldsymbol{x})=-x_{2}-x_{3}+x_{4},\
\phi_{3}(\boldsymbol{k},\boldsymbol{x})=x_{4}+x_{5}.$

One should note that $\boldsymbol{g}_{4}(\boldsymbol{k},\boldsymbol{x})$ is
not conservative as $\curl\boldsymbol{g}_{4}\neq 0$. Syzygies of the
generators of all other branches of comprehensive Gröbner system are not
conservative, and therefore, there is no polynomial conservation law. We want
to point out that the above Gröbner basis of the generic branch is not
reduced. If one would continue the computations with the reduced Gröbner basis

$\left(\\{k_{2}k_{4}\neq
0\\},\\{k_{4}x_{2}x_{4}-k_{5}x_{5},k_{2}x_{1}x_{3}-k_{3}x_{2}\\}\right),$

one would end up with a basis of the syzygy module formed by the trivial
syzygy

$(k_{2}x_{1}x_{3}-k_{3}x_{2},-k_{4}x_{2}x_{4}+k_{5}x_{5}).$

This syzygy is not irrotational and so no parametric conservation law would be
found.

As discussed earlier in Remark 35 one can compute unconditional conservation
laws from the syzygy module by an additional Gröbner basis computation. On the
other hand, unconditional conservation laws can be computed via ansatz, i.e.,
guessing the conservation laws with indeterminate coefficients, and then
obtaining coefficients via definition of first integrability. This method has
been used for computing polynomial first integrals of dynamical systems [32].
Computations based on an ansatz avoids Gröbner basis (and syzygy) computations
in the ring $\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$, as well as the
computation of additional unnecessary parametric conservation laws. However,
the matrices obtained via this method can be huge and the computations can be
hard to carry on. In the following, we explain this method for unconditional
conservation law.

Given a degree bound $d$, consider

$\phi(\boldsymbol{x})=\sum_{0\neq|\boldsymbol{\alpha}|\leq
d}c_{\boldsymbol{\alpha}}\boldsymbol{x}^{\boldsymbol{\alpha}}$

where $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}_{0}^{n}$
and $c_{\boldsymbol{\alpha}}$ are indeterminates for all $\boldsymbol{\alpha}$
with $0\neq|\boldsymbol{\alpha}|\leq d$. Then the condition

$\gradient{\phi(\boldsymbol{x})}\cdot F=0$

leads to a linear system of equations over $\mathbb{Q}$ in the
$c_{\boldsymbol{\alpha}}$ obtained by equating the coefficients of the
monomials with respect to the variables $\boldsymbol{x}$ and $\boldsymbol{k}$
in $\gradient{\phi(\boldsymbol{x})}\cdot F$ to zero. The solution space of the
linear system defines the $\mathbb{Q}$-vector space of unconditional
conservation laws up to order $d$. Algorithm 8 below uses this idea. If the
solution space is empty, then there are no unconditional conservation laws up
to order $d$. In this case one might increase the degree bound $d$ and run the
algorithm again.

0: 1\. Degree bound $d\in\mathbb{N}$2\. A polynomial ODE system
$F(\boldsymbol{k},\boldsymbol{x})=(f_{1}(\boldsymbol{k},\boldsymbol{x}),\dots,f_{n}(\boldsymbol{k},\boldsymbol{x}))$
whose r.h.s. $f_{i}(\boldsymbol{k},\boldsymbol{x})$ are polynomials in
$\mathbb{Q}[\boldsymbol{k},\boldsymbol{x}]$.

0: A set $\mathcal{S}$ of unconditional polynomial conservation laws for
$F(\boldsymbol{k},\boldsymbol{x})$ up to degree $d$.

1: (Ansatz)Let  $\phi(\boldsymbol{x})=\sum_{0\neq|\boldsymbol{\alpha}|\leq
d}c_{\boldsymbol{\alpha}}\boldsymbol{x}^{\boldsymbol{\alpha}},$ where
coefficients $c_{\boldsymbol{\alpha}}$ are indeterminate.

2: Obtain linear system of equations over $\mathbb{Q}$ in the
$c_{\boldsymbol{\alpha}}$ obtained by equating the coefficients of the
monomials with respect to the variables $\boldsymbol{x}$ and $\boldsymbol{k}$
in $\gradient{\phi(\boldsymbol{x})}\cdot F$ to zero.

3: Find the set of real solutions of the linear equations and obtain
$c_{\boldsymbol{\alpha}}$ and consequently
$\phi(\boldsymbol{k},\boldsymbol{x})$.

4: return $\mathcal{S}$, The set of conservation laws $\phi(\boldsymbol{x})$.

Algorithm 8 $\operatorname{PolynomialConservationLawsViaAnsatz}$

######  Example 40.

We consider the system of Example 38 and compute the set of unconditional
conservation laws up to order $d=2$ using Algorithm 8. With

$\displaystyle\phi$ $\displaystyle=$ $\displaystyle
c_{002}x_{3}^{2}+c_{011}x_{2}x_{3}+c_{020}x_{2}^{2}+c_{101}x_{1}x_{3}+c_{110}x_{1}x_{2}+c_{200}x_{1}^{2}$
$\displaystyle+c_{001}x_{3}+c_{010}x_{2}+c_{100}x_{1}$

we compute

$\displaystyle\gradient{\phi(\boldsymbol{x})}\cdot F$ $\displaystyle=$
$\displaystyle
c_{101}x_{2}x_{1}^{2}+c_{110}x_{1}^{2}+c_{011}x_{1}x_{2}^{2}+2c_{002}x_{3}x_{1}x_{2}+c_{110}x_{2}^{2}+c_{100}x_{2}$
$\displaystyle+(c_{001}+2c_{020}+2c_{200})x_{1}x_{2}+c_{011}x_{3}x_{1}+c_{010}x_{1}+c_{101}x_{3}x_{2}.$

The linear system formed by the coefficients of the monomials in
$\boldsymbol{x}$ in the last equation has the two dimensional solution space
with basis $c_{200}=1$, $c_{001}=-2$, all other variables zero and
$c_{020}=1$, $c_{001}=-2$, all other variables zero. Therefore the
$\mathbb{Q}$-vector space of unconditional conservation laws up to order two
has basis

$\phi_{1}=x_{1}^{2}-2x_{3}\quad\mathrm{and}\quad\phi_{2}=x_{2}^{2}-2x_{3}.$

It is easily checked that the two conservation laws of Example 38 can be
linearly represented by this basis.

A parametric version of Algorithm 8 can be easily adapted. In order to do so,
one can ansatz

$\phi(\boldsymbol{k},\boldsymbol{x}):=\phi(\boldsymbol{k},\boldsymbol{x})=\sum_{0\neq|\boldsymbol{\beta}+\boldsymbol{\alpha}|\leq
d}c_{\boldsymbol{\alpha},\boldsymbol{\beta}}\boldsymbol{k}^{\boldsymbol{\beta}}\boldsymbol{x}^{\boldsymbol{\alpha}},$

and then carry on parametric linear algebra computations as in steps 2 and 3
of the algorithm. This is a future research direction.

## Availability

ODEbase is available at https://www.odebase.org/. The Singular libraries
polconslaw.lib and genpolconslaw.lib are available at Zenodo under the record
https://zenodo.org/record/7474523 and through Github at the following url
https://github.com/oradules/polynomial_conservation_laws.

## Supplementary Information

The Supplementary Information containing the SI Tables 1 and 2 is available at
Zenodo under the record https://doi.org/10.5281/zenodo.7494868.

## Acknowledgement

This work has been supported by the interdisciplinary bilateral project
ANR-17-CE40-0036 and DFG-391322026 SYMBIONT. Ovidiu Radulescu thanks Werner
Seiler for very insightful discussions and for his hospitality in the
University of Kassel where the ideas about approximate conservation laws were
initiated. Hamid Rahkooy would like to thank Deepak Kapur and Amir Hashemi for
the discussions about comprehensive Gröbner bases and syzygy computations. The
project SYMBIONT owes a lot to Andreas Weber who sadly left us in 2020, but
who is still present in our memories.

## References

  * [1] E. Ábrahám, J. Abbott, B. Becker, A. M. Bigatti, M. Brain, B. Buchberger, A. Cimatti, J. H. Davenport, M. England, P. Fontaine, S. Forrest, A. Griggio, D. Kroening, W. M. Seiler, and T. Sturm, $SC^{2}$: Satisfiability checking meets symbolic computation, in Proc. CICM 2016, M. Kohlhase, M. Johansson, B. Miller, L. de Moura, and F. Tompa, eds., vol. 9791 of LNCS, Springer, 2016, pp. 28–43.
  * [2] R. Aris and R. Mah, Independence of chemical reactions, Industrial & Engineering Chemistry Fundamentals, 2 (1963), pp. 90–94.
  * [3] P. Auger, R. B. de La Parra, J.-C. Poggiale, E. Sánchez, and L. Sanz, Aggregation methods in dynamical systems and applications in population and community dynamics, Physics of Life Reviews, 5 (2008), pp. 79–105.
  * [4] C. Barrett, P. Fontaine, and C. Tinelli, The SMT-LIB standard: Version 2.6, tech. report, Department of Computer Science, The University of Iowa, 2017.
  * [5] T. Becker, V. Weispfenning, and H. Kredel, Gröbner Bases, a Computational Approach to Commutative Algebra, vol. 141 of Graduate Texts in Mathematics, Springer, 1993.
  * [6] C. Berkesch and F.-O. Schreyer, Syzygies, finite length modules, and random curves, Commutative algebra and noncommutative algebraic geometry, 1 (2015), pp. 25–52.
  * [7] N. K. Bose, Gröbner Bases: An Algorithmic Method in Polynomial Ideal Theory, Springer Netherlands, Dordrecht, 1995, pp. 89–127, https://doi.org/10.1007/978-94-017-0275-1_4.
  * [8] F. Boulier, M. Lefranc, F. Lemaire, and P.-E. Morant, Model reduction of chemical reaction systems using elimination, Mathematics in Computer Science, 5 (2011), pp. 289–301.
  * [9] C. W. Brown, QEPCAD B: A program for computing with semi-algebraic sets using CADs, ACM SIGSAM Bulletin, 37 (2003), pp. 97–108, https://doi.org/10.1145/968708.968710.
  * [10] C. W. Brown and C. Gross, Efficient preprocessing methods for quantifier elimination, in Computer Algebra in Scientific Computing: 9th International Workshop, CASC 2006, Chisinau, Moldova, September 11–15, 2006, V. G. Ganzha, E. W. Mayr, and E. V. Vorozhtsov, eds., vol. 4194 of Lecture Notes in Computer Science, Springer, 2006, pp. 89–100.
  * [11] B. Buchberger, Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal, doctoral dissertation, Mathematical Institute, University of Innsbruck, Innsbruck, Austria, 1965.
  * [12] V. Chelliah, N. Juty, I. Ajmera, R. Ali, M. Dumousseau, M. Glont, M. Hucka, G. Jalowicki, S. Keating, V. Knight-Schrijver, A. Lloret-Villas, K. Nath Natarajan, J.-B. Pettit, N. Rodriguez, M. Schubert, S. M. Wimalaratne, Y. Zhao, H. Hermjakob, N. Le Novère, and C. Laibe, BioModels: Ten-year anniversary, Nucl. Acids Res., 43 (2015), p. D542–D548, https://doi.org/10.1093/nar/gku1181.
  * [13] G. Chèze, Computation of darboux polynomials and rational first integrals with bounded degree in polynomial time, Journal of Complexity, 27 (2011), pp. 246–262.
  * [14] D. A. Cox, J. Little, and D. O’shea, Using algebraic geometry, vol. 185, Springer Science & Business Media, 2006.
  * [15] M. D. Darmian, A. Hashemi, and A. Montes, Erratum to ”a new algorithm for discussing Gröbner bases with parameters” [J. symbolic comput. 33(1-2) (2002) 183-208], J. Symb. Comput., 46 (2011), pp. 1187–1188, https://doi.org/10.1016/j.jsc.2011.05.002.
  * [16] W. Decker, G.-M. Greuel, G. Pfister, and H. Schönemann, Singular 4-2-0 — A computer algebra system for polynomial computations.  http://www.singular.uni-kl.de, 2020.
  * [17] A. Desoeuvres, A. Iosif, C. Lüders, O. Radulescu, H. Rahkooy, M. Seiß, and T. Sturm, Reduction of chemical reaction networks with approximate conservation laws, 2022, https://doi.org/10.48550/ARXIV.2212.13474, https://arxiv.org/abs/2212.13474.
  * [18] A. Dolzmann and T. Sturm, Redlog: Computer algebra meets computer logic, ACM SIGSAM Bulletin, 31 (1997), pp. 2–9.
  * [19] A. Dolzmann and T. Sturm, Redlog: Computer algebra meets computer logic, ACM SIGSAM Bulletin, 31 (1997), pp. 2–9, https://doi.org/10.1145/261320.261324.
  * [20] A. Dolzmann and T. Sturm, Simplification of quantifier-free formulae over ordered fields, J. Symb. Comput., 24 (1997), pp. 209–231.
  * [21] A. Dolzmann, T. Sturm, and V. Weispfenning, Real quantifier elimination in practice, in Algorithmic Algebra and Number Theory, B. H. Matzat, G.-M. Greuel, and G. Hiss, eds., Springer, 1998, pp. 221–247, https://doi.org/10.1007/978-3-642-59932-3_11.
  * [22] M. Eigen and P. Schuster, A principle of natural self-organization, Naturwissenschaften, 64 (1977), pp. 541–565.
  * [23] D. Eisenbud, The geometry of syzygies: A second course in commutative algebra and algebraic geometry, vol. 229, Springer, 2005.
  * [24] M. Feinberg, Foundations of Chemical Reaction Network Theory, vol. 202 of Applied Mathematical Sciences, Springer, 2019, https://doi.org/10.1007/978-3-030-03858-8.
  * [25] M. Feinberg and F. J. Horn, Dynamics of open chemical systems and the algebraic structure of the underlying reaction network, Chemical Engineering Science, 29 (1974), pp. 775–787.
  * [26] E. Feliu and C. Wiuf, Preclusion of switch behavior in networks with mass-action kinetics, Applied Mathematics and Computation, 219 (2012), pp. 1449–1467.
  * [27] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, Journal of differential equations, 31 (1979), pp. 53–98.
  * [28] S. Gao, F. V. IV, and M. Wang, A new framework for computing Gröbner bases, Math. Comput., 85 (2016), pp. 449–465, https://doi.org/10.1090/mcom/2969.
  * [29] A. Goeke, S. Walcher, and E. Zerz, Classical quasi-steady state reduction—a mathematical characterization, Physica D: Nonlinear Phenomena, 345 (2017), pp. 11–26, https://doi.org/10.1016/j.physd.2016.12.002.
  * [30] L. Goldman, Integrals of multinomial systems of ordinary differential equations, Journal of pure and applied algebra, 45 (1987), pp. 225–240.
  * [31] A. N. Gorban, O. Radulescu, and A. Y. Zinovyev, Asymptotology of chemical reaction networks, Chemical Engineering Science, 65 (2010), pp. 2310–2324.
  * [32] A. Goriely, Integrability and nonintegrability of dynamical systems, vol. 19, World Scientific, 2001.
  * [33] A. Hashemi, M. D. Darmian, and M. Barkhordar, Gröbner systems conversion, Math. Comput. Sci., 11 (2017), pp. 61–77, https://doi.org/10.1007/s11786-017-0295-3.
  * [34] M. Hauke, Signature Gröbner bases A comprehensive survey and a new algorithmic approach, master thesis, Graz University of Technology, Graz, Austria, 2020.
  * [35] A. C. Hearn, Reduce: A user-oriented interactive system for algebraic simplification, in Proc. Symposium on Interactive Systems for Experimental Applied Mathematics, ACM, 1967, https://doi.org/10.1145/2402536.2402544.
  * [36] A. C. Hearn, Reduce: The first forty years, in Algorithmic Algebra and Logic: Proceedings of the A3L 2005, BOD, Norderstedt, Germany, 2005, pp. 19–24.
  * [37] D. kapur, An approach for solving systems of parametric polynomial equations, Principles and Practice of Constraint Programming, (1995), p. 217–224.
  * [38] D. Kapur, Comprehensive Gröbner basis theory for a parametric polynomial ideal and the associated completion algorithm, J. Syst. Sci. Complex., 30 (2017), pp. 196–233, https://doi.org/10.1007/s11424-017-6337-8.
  * [39] D. Kapur, Y. Sun, and D. Wang, An efficient method for computing comprehensive Gröbner bases, J. Symb. Comput., 52 (2013), pp. 124–142, https://doi.org/10.1016/j.jsc.2012.05.015.
  * [40] M. Košta, New Concepts for Real Quantifier Elimination by Virtual Substitution, doctoral dissertation, Saarland University, Germany, Dec. 2016, https://doi.org/10.22028/D291-26679.
  * [41] M. Kreuzer and L. Robbiano, Computational Commutative Algebra 1, Springer Publishing Company, Incorporated, 2008.
  * [42] N. Kruff, C. Lüders, O. Radulescu, T. Sturm, and S. Walcher, Algorithmic reduction of biological networks with multiple time scales, Mathematics in Computer Science, in press, (2020), https://arxiv.org/abs/2010.10129.
  * [43] A. M. Lozano, The Gröbner Cover, vol. 27 of Algorithms and Computation in Mathematics, Springer, 2018, https://doi.org/10.1007/978-3-030-03904-2.
  * [44] D. Lu, Y. Sun, and D. Wang, A survey on algorithms for computing comprehensive Gröbner systems and comprehensive Gröbner bases, J. Syst. Sci. Complex., 32 (2019), pp. 234–255, https://doi.org/10.1007/s11424-019-8357-z.
  * [45] C. Lüders, T. Sturm, and O. Radulescu, Odebase: a repository of ode systems for systems biology, Bioinformatics Advances, 2 (2022), p. vbac027.
  * [46] C. Lüders, T. Sturm, and O. Radulescu, ODEbase: A repository of ODE systems for systems biology, CoRR, abs/2201.08980 (2022), http://arxiv.org/abs/2201.08980, https://arxiv.org/abs/2201.08980.
  * [47] A. Mahdi, A. Ferragut, C. Valls, and C. Wiuf, Conservation laws in biochemical reaction networks, SIAM Journal on Applied Dynamical Systems, 16 (2017), pp. 2213–2232.
  * [48] Y.-K. Man, Computing closed form solutions of first order odes using the prelle-singer procedure, Journal of Symbolic Computation, 16 (1993), pp. 423–443.
  * [49] A. Montes, A new algorithm for discussing Gröbner bases with parameters, J. Symb. Comput., 33 (2002), pp. 183–208, https://doi.org/10.1006/jsco.2001.0504.
  * [50] K. Nabeshima, Generic Gröbner basis of a parametric ideal and its application to a comprehensive Gröbner system.  Applications of Computer Algebra Conference, https://www.rs.tus.ac.jp/nabeshima/ACA2022/2022aca/nabeshima22.pdf, 2022\.
  * [51] R. Nieuwenhuis, A. Oliveras, and C. Tinelli, Solving SAT and SAT modulo theories: From an abstract Davis–Putnam–Logemann–Loveland procedure to DPLL(T), J. ACM, 53 (2006), pp. 937–977, https://doi.org/10.1145/1217856.1217859.
  * [52] M. J. Prelle and M. F. Singer, Elementary first integrals of differential equations, Transactions of the American Mathematical Society, 279 (1983), pp. 215–229.
  * [53] O. Radulescu, A. N. Gorban, A. Zinovyev, and V. Noel, Reduction of dynamical biochemical reactions networks in computational biology, Frontiers in genetics, 3 (2012), p. 131.
  * [54] F. Schlögl, Chemical reaction models for non-equilibrium phase transitions, Zeitschrift für physik, 253 (1972), pp. 147–161.
  * [55] K. R. Schneider and T. Wilhelm, Model reduction by extended quasi-steady-state approximation, Journal of mathematical biology, 40 (2000), pp. 443–450.
  * [56] F. O. Schreyer, A standard basis approach to syzygies of canonical curves., Journal für die reine und angewandte Mathematik, 421 (1991), pp. 83–124.
  * [57] S. Schuster and T. Höfer, Determining all extreme semi-positive conservation relations in chemical reaction systems: a test criterion for conservativity, Journal of the Chemical Society, Faraday Transactions, 87 (1991), pp. 2561–2566.
  * [58] A. Seidl, Cylindrical Decomposition Under Application-Oriented Paradigms, doctoral dissertation, University of Passau, Germany, Mar. 2006.
  * [59] A. M. Seidl and T. Sturm, Boolean quantification in a first-order context, in Computer Algebra in Scientific Computing. Proceedings of the CASC 2003, V. G. Ganzha, E. W. Mayr, and E. V. Vorozhtsov, eds., Institut für Informatik, Technische Universität München, München, Germany, 2003, pp. 329–345.
  * [60] N. N. Semenov, Some problems relating to chain reactions and to the theory of combustion, Nobel Lecture, (1956), pp. 487–514.
  * [61] S. Soliman, Invariants and other structural properties of biochemical models as a constraint satisfaction problem, Algorithms for Molecular Biology, 7 (2012), pp. 1–9.
  * [62] M. Spivak, Calculus on manifolds: a modern approach to classical theorems of advanced calculus, CRC press, 2018.
  * [63] T. Sturm, Redlog online resources for applied quantifier elimination, Acta Acad. Ab., 67 (2007), pp. 177–191.
  * [64] T. Sturm, A survey of some methods for real quantifier elimination, decision, and satisfiability and their applications, Math. Comput. Sci., 11 (2017), pp. 483–502, https://doi.org/10.1007/s11786-017-0319-z.
  * [65] T. Sturm, Thirty years of virtual substitution, in Proc. ISSAC 2018, ACM, 2018, pp. 11–16.
  * [66] A. Suzuki and Y. Sato, An alternative approach to comprehensive Gröbner bases, J. Symb. Comput., 36 (2003), pp. 649–667, https://doi.org/10.1016/S0747-7171(03)00098-1.
  * [67] N. Theret, J. Feret, A. Hodgkinson, P. Boutillier, P. Vignet, and O. Radulescu, Integrative models for tgf-$\beta$ signaling and extracellular matrix, in Extracellular Matrix Omics, Springer, 2020, pp. 209–225.
  * [68] A. N. Tikhonov, Systems of differential equations containing small parameters in the derivatives, Mat. Sb. (N. S.), 73 (1952), pp. 575–586.
  * [69] F. Volny, New Algorithms for Computing Groebner Bases, doctoral dissertation, Clemson University, Clemson, USA, 2011.
  * [70] J. Wei and C. D. Prater, The structure and analysis of complex reaction systems, in Advances in Catalysis, vol. 13, Elsevier, 1962, pp. 203–392.
  * [71] V. Weispfenning, Comprehensive Gröbner bases, J. Symb. Comput., 14 (1992), pp. 1–30, https://doi.org/10.1016/0747-7171(92)90023-W.
  * [72] G. Yablonskii, V. Bykov, V. Elokhin, and A. Gorban, Kinetic models of catalytic reactions, Elsevier, 1991.