Realizations of kinetic differential equations

The induced kinetic differential equation of a reaction network endowed with mass action type kinetics is a system of polynomial differential equations. The problem studied here is: Given a polynomial differential equation, is it possible to find a network which induces the equation? If yes, can we find a network with some chemically relevant properties (implying also important dynamic consequences), such as reversibility, weak reversibility, zero deficiency, detailed balancing, complex balancing, mass conservation, etc. The constructive answers presented to a series of questions of the above type are useful when fitting a differential equation to measurements, or when trying to find out the dynamic behavior of the solutions of a differential equation. It turns out that some of the results can be applied when trying to solve purely mathematical problems, like the existence of positive solutions to polynomial equation.


Introduction
Our goal here is to find reaction networks inducing a given polynomial differential equation (or classes of equations with coefficients as symbolic parameters) with as many good properties (e.g. weak reversibility, reversibility, small deficiency, small number of linkage classes or reaction steps, mass conservation, etc.) as possible.
It has been known that under mass action kinetics, the necessary and sufficient condition for realizability is a sign structure in the system of polynomial differential equations [18]. (See Lemma 1 below.) Once we know a system can arise from a mass action system, the question still lies: can the system arise from a mass action system with good properties? This is what this present work looks at.
There are several sources of motivation for this problem.
1. The most practical being that having fitted a system of differential equations to data, you may wonder whether the obtained equations can be interpreted as the mass action type deterministic model of an appropriate reaction network. 2. There is an internal requirement within this branch of science: One should like to know as much as possible about the structure of differential equations that arises from modelling a chemical system. 3. Given a system of polynomial differential equations in any field of pure or applied mathematics, you may wish to have statements on stability or oscillation, similar to those offered by the Zero Deficiency Theorem [16], Volpert's theorem [43], or the Global Attractor Theorem [5]. Then it comes in handy to see that your system of differential equations belongs to a well behaving class. 4. Last but not least: results of formal reaction kinetics (to use the expression of [1,2]) on the stationary point may offer solutions to problems in algebra, e.g. showing the existence of positive roots of a polynomial easier than by classical methods of algebraic geometry [12,26] (e.g. based on the result by [3]) if the polynomial is known to be the right hand side of the induced kinetic differential equation of a reversible or weakly reversible reaction network.
The structure of our paper is as follows. Section 2 introduces the essential concepts of reaction networks and mass action systems. Section 3 formulates the problem we are interested in, that of realizability of kinetic differential equations. Section 4 treats two special cases: finding realizations for compartmental models (defined later) and weakly reversible networks. Section 5 focuses on the general problem of realizability. We first review existing algorithms available. Then we outline several procedures on a reaction network that preserves the system of differential equations, including adding and removing vertices from the reaction graph. Section 6 explores the relation between weakly reversible and complex balanced realizations. Here we work with families of symbolic kinetic differential equations. In Section 7, we prove the geometric meaning of zero deficiency, and discuss when is a realization unique. We also pose several conjectures that may be of interest to the reader. An Appendix contains a very large number of enlightening examples. The present paper intends to be a review of the realizability problem, while containing some new results. A few remarks on notation and the use of words are in order. 1. The set of vectors in R M with strictly positive coordinates is denoted by R M >0 . Also x y := x y 1 1 x y 2 2 · · · x y M M for any x i > 0 and y i ≥ 0. If Y = (y 1 , y 2 , . . . , y N ), then x Y = (x y 1 , x y 2 , . . . , x y N ) T .
2. We use differential equation to mean a system of ordinary differential equations. In the same vein, we use polynomial to usually mean a system of polynomial equations, i.e. a vector-valued function on R M >0 . By monomial we mean a scalar-valued function on R M >0 . 3. We say a reaction network induces a differential equation, if we formulate the mass action type differential equation with given reaction rate coefficients, and a differential equation is realized by a reaction network, if the network induces the given differential equation.

Chemical reaction networks and mass action systems
Here we recapitulate shortly the necessary concepts of the mathematical description of reaction networks; for more details, see e.g. [16,42].
A reaction network consists of a set of R reaction steps among the chemical species X(1), X(2), . . . , X(M): The positive number k r is the rate coefficient * for the r-th reaction.
On the two sides of arrows in (2.1), which represent a reaction step, one has formal linear combinations of the species, called complexes, the left one being the reactant complex, and the right one being the product complex. Associated to each complex is the vector of the coefficients of these linear combinations. By an abuse of notation, we will refer to this vector as a complex. The difference between the product complex and the reactant complex is the reaction vector. These vectors form the columns of the stoichiometric matrix. The stoichiometric space S is the linear span of the reaction vectors, where S = dim S. (This number can also be interpreted as the number of independent reaction steps.) The number of complexes is usually denoted by N.
With the reaction sets and complexes as described above, each reaction network is naturally associated to a directed graph. Definition 1. The reaction graph (or Feinberg-Horn-Jackson graph) of a given reaction network is a finite directed graph with no self-loops, where the vertices are the complexes of the reaction network, and the edges are precisely the reactions.
The connected components of a reaction graph are also called linkage classes, and L denotes the number of linkage classes in a reaction network. The deficiency of the reaction network is δ := N−L−S , where N is the number of complexes (i.e., number of nodes in the reaction graph), and S = dim S. The strong components of the reaction graph are also called strong linkage classes. An ergodic component (or strong terminal class) is a strong component with no reaction step such that the reactant complex is in the component while the product complex is outside of it, i.e., an ergodic component is a strongly connected component of the raction graph. The number of ergodic components will be denoted by T .
A reaction network is reversible if its reaction graph as a relation is symmetric, i.e., i → j is a reaction if and only if j → i is a reaction. It is weakly reversible, if the transitive closure of its reaction graph is symmetric; or if all the reaction arrows are edges of a directed cycle in the graph; or all of its strong components are ergodic. Thus, weak reversibility implies T = L, but the converse is false [17].
The reversible reaction network (with rate coefficients k r+ , k r− for the r-th pair of reversible reactions) is detailed balanced at the positive stationary concentration x * if all the steps proceed with the same rate in both directions. In other words, x * is detailed balanced if dx dt = 0 implies k r+ (x * ) α r = k r− (x * ) β r , (r = 1, 2, . . . , P).
A detailed balanced stationary point exists if and only if for all u ∈ Ker(γ) one has κ u = 1, where The induced kinetic differential equation can be written in many different forms, see e.g. [42,Sec. 6.3]. We introduce the form that will be used frequently below.
For a reaction graph with N complexes, let their associated vectors be columns of the matrix Y := (y 1 , y 2 , . . . , y N ). Implicitly we are imposing an ordering on the set of complexes. For a reaction i → j (from the i-th complex to the jth complex), assume that the rate coefficient is k ji > 0. Define the matrix K to be so that K is a matrix with nonnegative entries, and 0 along the diagonal. Then the column-conserving matrix A k := K − diag(K 1) is the (weighted) Laplacian of the reaction graph. With these matrices defined, the induced kinetic differential equation (2.2) of the reaction network can be rewritten in the following form:ẋ where x Y is the vector of monomials (x y 1 , x y 2 , . . . , x y N ) T . Note that the system of differential equations (2.5) is completely determined by the complexes (columns of Y) and the connectivity and rate coefficients of the reaction graph (given by A k ). Thus, we call (Y, A k ) a realization of the differential equation (2.5).
The reaction system (2.5) is complex balanced at the positive stationary concentration x * , if all the complexes are destroyed and formed with the same rate at this concentration: If one remains within the framework of mass action type kinetics -as we do here -then the deficiency zero theorem provides a simple structural necessary and sufficient condition for complex balancing [14,19,20]. Theorem 1. A mass action system is complex balanced for all choices of the rate coefficients if and only if it is weakly reversible and its deficiency is zero.
One might notice that in the deficiency zero theorem, it is the mass action system that is described as complex balanced, not one of its stationary concentrations. This is because complex balancing implies a host of good structural and dynamical properties, including: 1. If one of the stationary concentration is complex balanced, then all stationary concentrations of the system are complex balanced. 2. There is exactly one stationary concentration within every stoichiometric compatibility class. 3. Let x * be any complex balanced stationary concentration. The function defined on R n >0 is a Lyapunov function, with global minimum at x = x * . 4. Every positive stationary concentration is locally asymptotically stable within its stoichiometric compatibility class. 5. The set of complex balanced stationary concentrations Z k satisfies the equation ln Z k = ln x * + S ⊥ .

Formulation of problems
In defining the system of differential equations for a mass action system, notice that given a set of reactions and rate coefficients (or equivalently, given a reaction graph and positive weights for each edge), it is straightforward to write down the differential equation. Unfortunately, given a system of differential equations (satisfying a mild sign condition introduced below), there is in general no unique reaction graph that is associated to it [6].
A natural question arises: given a differential equation, what are the possible reaction graphs that could generate it if we assume mass action kinetics. Since all the information about a reaction graph is contained inside (Y, A k ), in other words, we are asking what are the possible realizations (Y, A k ) of a given system of differential equations. Different realizations of the same differential equation are said to be dynamical equivalent. Dynamical equivalence is written as linear relations in [7].
Before searching for dynamically equivalent realizations, we must first ask whether or not a given system of differential equations can arise from a mass action system. In other words, is it even possible to find one realization. We first provide a necessary and sufficient condition for realizability.
Assume that we are given a system of differential equations of the forṁ where Z := (z 1 , z 2 , . . . , z N ), and Y := (y 1 , y 2 , . . . , y N ), and x y := x y 1 1 x y 2 2 · · · x y M M for any x i > 0. Definition 2. We say that a system of differential equations (3.1) is kinetic if A system being kinetic reflects the physical assumption that in order to lose a chemical species via a reaction, that species must participate as a reactant. For example, the systeṁ is not kinetic because of the term −xz inẏ does not involve the variable y.
It turns out that a differential equation is kinetic if and only if the differential equation (3.1) comes from a mass action system. This has been shown in [18] using a constructive method: Lemma 1 (Hungarian Lemma). A polynomial system of differential equations of the form (3.1) can be realized by a mass action system if and only if the system of differential equations is kinetic.
One realization (the canonical realization) is the collection of reactions Example 1. If the differential equationẋ is given (with k 1 , k 2 > 0), then it is possible to find an inducing reaction network with only two complexes defined by the exponents of the monomials present in the right hand side: The right hand side contains two monomials, and there exists a realization with two complexes, which is also reversible, and is deficiency zero. This also happens to be the canonical realization.
However, the canonical realization is usually quite complicated, contains too many complexes and reaction steps (see Example 13 in the Appendix). Its main advantage is that it can be automatically constructed and used as a starting point. Now we present the set of problems we are interested in (also posed earlier in [13,Section 4.7]).

Problem.
1. Given a kinetic differential equationẋ = Z · x Y , find A k such that Y · A k = Z. 2. Find reversible or weakly reversible; complex balanced or detailed balanced realizations; mass conserving realizations. Decide if such a realization exists. 3. Find realizations with minimal deficiency, with minimal number of reaction steps, or with minimal number of complexes. 4. Without further restrictions, the question about uniqueness arises at all levels. 5. In all of the above problems, we may add more columns to Y [34]. For example, if the first problem has no solution, we may consider the following: find Y whose first N columns form precisely Y and find A k such that Y · A k · x Y = Z · x Y . In this case we can also ask what is the minimum number of columns that must be added.
We will answer partially the questions above, especially pertaining to reversible, weakly reversible, complex balanced and detailed balanced realizations. We also consider the uniqueness in Section 7 and mass conserving realizations in the Appendix.
In considering the uniqueness question, the concepts of sparse and dense realizations have been useful [35]. The aim is to find an M × N matrix of nonnegative integer components Y and the Kirchoff matrix A k giving the product Z = Y · A k . The pair (Y, A k ) is a realization of the matrix Z. Since (Y, A k ) uniquely identifies the reaction graph, when we refer to a subgraph of the realization, we mean a subgraph of the reaction graph corresponding to (Y, A k ).
Definition 3. If the number of zeros in the above setting in A k is maximal, then a sparse realization has been obtained. If the number of zeros in A k is minimal then a dense realization has been obtained.
The following properties of dense and sparse realizations of a given chemical reaction network were proved in [38]. Suppose Z is given and fix the matrix Y, i.e., fix the set of complexes from which we will search for realizations.
1. The reaction graph of the dense realizations is unique. 2. The reaction graph of any dynamically equivalent realization is a subgraph of the dense realization. 3. The reaction graph of a system of differential equations is unique if and only if the reaction graph of the dense and sparse realizations are identical.
The above statements hold when the word "realization" is replaced by weakly reversible realization [8]. The above results were extended to the case of constrained realizations, where the rate coefficients fulfil arbitrary polytopic constraints [9]. An important special case of this is when a subset of possible reactions is excluded from the network.
We have seen that the canonical representation may be the sparse realization. In some cases, based on Theorem 3 below we shall be able to find the dense realization. Later, we shall give a complete solution to a generalization of part 1 of Problem 3 for the case when the numerical values of the components of the matrices Z and Y are known. We also review results for the solution of the other related problems.
To date, the problem of finding a realization has used tools and techniques from different mathematical disciplines, ranging from graph theory, combinatorics, linear algebra, to optimization.

Realizations of compartmental models and reversible systems
In this section, we consider two class of examples: compartmental models and weakly reversible networks. Because of the simpler mathematical structure of compartmental models, the realization problem can be solved even if the rate coefficients are unspecified parameters, i.e., symbolic. Then we will look at an ad hoc method to find reversible realizations.
A compartmental system are built using the reaction steps First we give a complete, although ineffective solution of the first task of Problem 3 in a slightly modified form. Modification means that the empty complex will be dealt with separately, because in some calculations this implies simplification. First, we need two definitions. assuming that k nq > 0 (1 ≤ n q ≤ N) for all the reaction steps present. Let us introduce the notation (slightly different from the usual) so that there exists a compartmental matrix B ∈ R N×N and a nonnegative vector v ∈ (R + 0 ) M so that hold. If u = 0, then B is a Kirchhoff matrix. (One may call (4.4) outflow, (4.5) inflow, and they are missing if u = 0 or v = 0, respectively.) A partial answer to the realization problem follows.
Theorem 3. Suppose that the differential equatioṅ has the property that there exist a compartmental matrix B and a nonnegative vector v ∈ (R + 0 ) M so that (4.8) holds. Then the equation has a realization of the form (4.3)-(4.5).
Proof. First we prove that the conditions of the theorem imply the differential equation is kinetic. Suppose Z m n < 0 for some indices m, n. Then, y m q B q n = N q n y m q B q n + y m n B n n , and the last sum being nonnegative, y m n B n n should be negative, thus y m n > 0.
The components of Y and v are nonnegative, thus the components of b = Yv are nonnegative, as well.
Next we construct the inducing realization in the following way. Let the reaction rate coefficients of the reaction network (4.3)-(4.5) be defined as follows: Then simple verification proves the statement.
It may happen that a given kinetic differential equation has zero, one or an infinite number of realizations with complexes, or columns of Y, determined by the exponents of its monomials, see Examples 11, 12 and 15 and 2, respectively. For some nonlinear differential equations one can provide a zero deficiency reaction network. Additional conditions may allow to have a realization which is also reversible or weakly reversible. One of the early approaches to this problem was [41,Theorem 9], see also [42,Theorem 11.12] providing a necessary and sufficient condition that a generalized compartmental system induces a given differential equation. If the coefficients of the kinetic polynomial ODEs are known, then a deficiency zero realization can be determined via mixed integer linear programming, if it exists [27]. Moreover, weakly reversible deficiency zero realizations can be determined using simple linear programming [36].
is the induced kinetic differential equation of a Proof. Necessity is obvious. Instead of direct calculations sufficiency can be learned from the application of Theorem 3, as in this case Y -being diagonal with positive components -is invertible, and its inverse can easily be calculated.
1. A necessary and sufficient condition of reversibility in the closed case is sign symmetry of the coefficient matrix. Similar conditions can be formulated for the other two cases, as well. 2. Once we have a zero deficiency realization, it is easy to check detailed balancing: only the circuit conditions [15] are to be checked.
Let us study reversibility, another relevant property, in more detail, applying different methodswithout taking care of the value of deficiency.
An equivalent and useful characterization of the induced kinetic differential equation of reversible reaction networks seems to be missing. What we only have is a heuristic algorithm without the proof that it leads to a reversible realization if there exists one, and the realization is independent from the order of the steps.
1. Consider the polynomial equationẋ with b, g, s > 0 (a transform of the first repressilator model in [12]). The canonical representation of this equation is shown in Fig. 1. It is a reversible reaction network with exactly those seven complexes which are given as exponents of the monomials. It consists of a single linkage class and has a three-dimensional stoichiometric space, and it is also a sparse realization. Reversibility allows us to deduce the existence of a positive stationary point by the theorems of [3,28]. In the third part below a more systematic approach to this example will show an infinite number of realizations (not necessarily reversible), actually we shall receive the dense realization. 2. To find a reversible realization one can also use a heuristic: one tries to collect pairs of reactant and product complex vectors. Suppose our goal is to find a reversible reaction network inducing (4.11)-(4.13). Let us rewrite it in the following form Now we consider the exponent of the monomial as a reactant vector α and the vector with which this monomial has been multiplied, as the reaction vector β − α, then one should have β = , thus we need a term where the exponent of the monomial is β, and it is Repeating steps of this kind we arrive at the decomposition which again leads to the network in Fig. 1.

Conjecture 1.
If there exists a reversible realization then the above process (independently from the order of actions) ends and necessarily provides one of them which may or may not depend on the order of steps.
3. Finally, one can apply Theorem 3 to get all the realizations (with the exponents as complexes only) and choose from them the reversible one(s). Suppose again that the Eqs. (4.11)-(4.13) are given, and our goal is to find a reaction network inducing it using Theorem 3. In this case we have and we have to find a compartmental solution in B of the equation Z = YB, or which for the coordinates has the following consequences.
and then the equations reduce to having the solution implying that u 4 = u 5 = u 6 = 0.
As to the inflow: the components of the solution v to Yv = b should fulfil which means that any v defines a solution for which Summing up, Fig. 1 shows the sparse realization while if all the reaction rate coefficients in Fig.  (4.17) are positive, one arrives at the dense realization. Fig. (4.17) also shows that the condition of reversibility is this: α = β = γ = s and v 4 = v 5 = v 6 = 0, and that it is exactly the sparse realization which is reversible.
In this example we were lucky to receive a canonical realization with so many nice properties. The reader can easily verify that this is the case with the second three dimensional repressilator model and with the six dimensional one of [12]: the canonic realization is again reversible thus assures the existence of a positive stationary point.
Let us emphasize that via a known result in reaction kinetics we were able to assert the existence of a positive solution to a second degree polynomial (in three variables), or to put it another way, we could apply formal reaction kinetics in mathematics.
As to uniqueness of the stationary point: it does not follow from the mentioned theorems. One may try several methods, see e.g. [24]. However, uniqueness of the reversible realization itself does follow (in this special case) from the method based on Theorem 3.

On computing realizations of kinetic equations
In the previous section, we looked at two classes of systems with certain network properties, namely that of (generalized) compartmental models (i.e., complexes consisting of a single species) and that of reversible networks. In this section, we will provide an algorithmic approach to the realizability problem when the coefficients in the ODE are given. To date, there are algorithms based on linear programming (LP) and mixed integer linear programming (MILP) techniques. These are reviewed in Section 5.1.
In the compartmental models and in the method for finding reversible realizations, we only used the complexes that showed up as exponents of monomials in the kinetic differential equation. This may not necessarily be the case for more general examples [34]. However, once a set of complexes has been decided on (and must include the exponents of monomials), the algorithms presented below can be used to compute realizations with desirable structural properties.
In Section 5.2, we consider ways to modify an existing realization while preserving dynamical equivalence, including adding to/subtracting from the set of complexes. Finally, we consider cases when a realization with certain properties (e.g., complex-balanced) is independent of the choice of the set of complexes.

Review: algorithms to find realizations
Qualitative properties can be rewritten in terms of a mixed integer linear programming (MILP) problem [38]. The works [39,21] determine weakly reversible realizations, while [36] provides reversible and detailed balanced realizations. The paper [33] finds detailed balanced subnetworks in reactions.
Recall the discussion about possibly enlarging the set of complexes (Problem 5) from Y to Y. We write the dynamical equivalence conditions for finding a realization (using linear programming method) with this predetermined complex set Y. (Note that the columns of Y must be columns of Y.) We want to find a realization of the kinetic differential equationẋ = Z · x Y with the complex set Y. In other words, we want to find Kirchoff matrix A k such that Add columns of 0 to the matrix Z and obtain Z, where the number of columns added is the number of additional columns in Y compared to Y. Then the equation (5.1) becomes Therefore, kinetic realizability using Y is equivalent to the feasibility of the following convex constraint set Next we review algorithms for finding realizations with certain structural properties, namely, reversible, weakly reversible, or complex balancing realizations. In equations (5.3)-(5.5) above, the matrices Y and Z do not show up explicitly. Hence, to simplify notation, in what follows below, we let Y denote the complexes to be considered, some of which may not appear explicitly in the differential equation, and Z defines a kinetic differential equation (with possibly columns of 0 so that Z is the same size as Y).
To impose reversibility is to consider the following feasibility problem [22]: The statement (5.9) is the condition that the reaction step i → j is present if and only if reaction j → i is present. This feasibility problem can be solved in the framework of mixed integer linear programming (MILP).
For weak reversibility [23]: 14) The inequalities (5.13) ensures that the structure of A k and A k are the same, i.e., a reaction is present in one if and only if it is present in the other. Since every weakly reversible network is complex balancing for some choice of parameter, the algorithm searches for a complex balanced A k (in equation (5.14)), which shares the network structure we want but there is no relations in terms of the rate coefficients we seek.
In this MILP problem the decision variables are Z, off-diagonals of A k , off-diagonals of another Kirchoff matrix A k , while Y is a given matrix, and ε 1 is a fixed small constant. (We will explore the relationship between finding weakly reversible and complex balancing realizations in Section 6.) More efficient methods to compute weakly reversible realizations are available; for example, a polynomial time algorithm can be found in [31].
For finding detailed balanced realizations, it is useful to have the numerical value of one stationary concentration x * [36]: Finally, we consider complex balanced realizations. For this purpose, it is useful to have the numerical value of one stationary concentration x * [36]:

Modifying network while preserving dynamical equivalence
In this section, we introduce several modifications to an existing realization while preserving dynamical equivalence. These include: 1. positive linear combinations of realizations, 2. adding new reactions, and 3. eliminating complexes.
This list is non-exhaustive, but we hope that some of the ideas presented in this section will motivate other types of modifications that preserve dynamical equivalences.
, to the kinetic differential equationẋ = Z · x Y . Then we can easily produce other realizations which may be simpler, or fulfil some further properties [11].
Theorem 5. Suppose that A 1 k ,. . . ,A I k are realizations of a kinetic differential equation, with complex set Y, and let A k : k defines a reversible realization for all i = 1, . . . , I, then A k defines a reversible realization. 4. If A i k defines a weakly reversible realization for all i = 1, . . . , I, then A k defines a weakly reversible realization.

Proof.
1. Obviously, the off-diagonal terms of A k are nonnegative. Also the column sums of A k are conserved. Thus A k is Kirchoff. We see that 2. This trivially follows from the calculation above.
k defines a weakly reversible realization, the reaction q → p must be part of a cycle, i.e., the off-diagonal terms in A i k corresponding to the edges in this cycle are positive, hence those positions in A k are also positive. In other words, this particular cycle appears in the graph defined by A k as well.
Now we look at modifying an existing realization by adding new reactions. Suppose that A k defines a realization with complexes y 1 , y 2 , . . . , y N , and suppose that y is a convex combination of the remaining complexes, i.e., y = n v n y i where v n ≥ 0 and n v n = 1.
Then we can add the following reactions: for every n such that v n 0, add a reaction from y to the complex y n , with rate coefficient v n . These reactions together contributes 0 to the differential equation.
We wonder if we can add reactions in the opposite direction, i.e., reactions of the form y n → y. There are cases when it is possible and cases when it is not. While we do not have an answer to this question, we illustrate the complexity using two examples.
Example 3. First we show an example where adding a complex situated in the convex cone of the original complexes seems to be superfluous: its effect on the induced kinetic differential equation can be expressed only using reaction steps among the original complexes. Consider the system 0 2X 2X + 2Y 2Y To add the reaction step 2X → X + Y to get the network 0 2X 2X + 2Y 2Y X + Y one must borrow from the reactions 2X → 0 and 2X → 2X + 2Y in the sense that The example above may lead us to believe that it is always possible to add a reaction from some y n on the boundary of the convex hall to some complex y in the interior of the convex hull of other complexes. The example below illustrates that this is not always the case, and further work is needed to explore when it is possible to add such a reaction. 2X + 2Y 2Y and tries to add the reaction step 2X → X + Y (note that the new complex X + Y is a convex combination of the old ones), there is no way of choosing the rate coefficients so as to get the induced kinetic differential equation of the extended one. Perhaps instead of adding or removing reactions from an existing realization, one may want to add or remove complexes from an existing realization. As noted in [34], it is always possible to add a complex that has not appeared in the current realization, as long as the weighted sum of reaction vectors is 0. For example, one can always add the reactions 0 k ← − X k − → 2X without changing the induced kinetic differential equation. The question lies, then, on when it is possible to eliminate a complex from your realization.
Theorem 6. Suppose that (Y, A k ) is a realization with N + 1 complexes y 1 , . . . , y N , y N+1 , where the complex y N+1 participates in at least one reaction. Suppose that the monomial x y N+1 does not appear explicitly in the kinetic differential equation. Then the differential equation can also be realized using only the complexes y 1 , . . . , y N .
Proof. We can write the complex matrix and Laplacian matrix as where u, v ∈ R N ≥0 . Then the system of ODEs can be written as Since the complex y N+1 admits some outcoming reactions, we have that N n=1 v n 0, and, using the fact that the coefficient of x y N+1 is zero, we find We conclude that y N+1 is a convex combination of the remaining complexes. More importantly, the vu gives a realization with only the complexes y 1 , . . . , y N .

Remark 2.
If x y appears explicitly as a term in the kinetic differential equation, then that complex cannot be removed.
After seeing the procedure of eliminating extraneous complexes, one may be tempted to remove all complexes that do not appear as exponents of monomials without ruining desirable properties (e.g., weak reversibility, complex balanced).
Indeed, if one is interested in finding a reversible, or weakly reversible, or detailed-balanced, or complex-balanced realization, any complex that does not appear in the differential equation (after simplification) can be eliminated from the network [7]. Theorem 7. A kinetic differential equationẋ = f(x) has a complex-balanced realization if and only if it has a complex-balanced realization where the complexes (the columns of Y) are precisely the monomials in f(x). Analogous results hold for detailed-balanced, weakly reversible, and reversible realizations.
As a consequence, the algorithms described in Section 5.1 can be used to determine whether a given system of kinetic differential equations admits a complex-balanced (or detailed-balanced, weakly reversible, or reversible) realization. This is a finite calculation as there is no need to compute using different sets of complexes.

Relation between weakly reversible and complex balanced realizations for symbolic kinetic equations
When we discussed the algorithm for computing weakly reversible realizations, the reader may have noticed that we introduced an additional Kirchoff matrix A k that does not appear in the realization (Y, A k ). This reflects the fact that there is a subtle algebraic relation between weak reversibility and complex balancing. In this section, we explore this relationship.
In order to study the possible network structures that can arise from kinetic differential equation to those with uncertain coefficient, we make the following definition: where sign(·) is the componentwise sign function.
According to [18], the system of differential equationsẋ = Z · x Y is kinetic if Z ∈ Z(Y, Z σ ) for any kinetically compatible Z σ . Moreover, for any Y ∈ Z M×N and kinetically compatible Z σ , the set Z(Y, Z σ ) is an unbounded convex cone in R M×N . where clearly Z σ is kinetically compatible with Y. Note that the kinetic matrix corresponding to the Lotka-Volterra equations is Z = α −β 0 0 δ −γ . In particular Z ∈ Z(Y, Z σ ).
We defined Z(Y, Z σ ) in order to study possible behaviours a kinetic differential equation may exhibit for some choice of rate coefficients. Definition 6. Let Z σ be kinetically compatible with the matrix Y. We say that Z(Y, Z σ ) has the capacity to be complex balanced (or capacity to be weakly reversible) if there exists Z ∈ Z(Y, Z σ ) such thatẋ = Z · x Y has a complex balanced (or weakly reversible) realization.
Theorem 8. Let Z σ be kinetically compatible with the matrix Y. The set Z(Y, Z σ ) has the capacity to be complex balanced if and only if it has the capacity to be weakly reversible.
Proof. It is known that complex balancing implies weak reversibility. Suppose there is some Z ∈ Z(Y, Z σ ) with a weakly reversible realization (Y, A k ). Then there is a positive vector p in the kernel of the Laplacian A k . Fix an arbitrary positive state x * . Consider the following scaled Laplacian matrix where the division is elementwise. Then i.e., A k defines a complex balanced realization with stationary concentration x * .

On deficiency zero realizations and uniqueness
Recall that the deficiency of a reaction network is the nonnegative integer δ = N − L − S , where N is the number of complexes, L is the number of linkage classes, and S is the dimension of the stoichiometric subspace. A network with lower deficiency tends to have nicer dynamical properties; for example, by the Deficiency Zero Theorem or Deficiency One Theorem. Let us formulate a necessary and sufficient condition on network structure for zero deficiency [16].
We say that the vectors y 0 , y 1 , . . . , y N ∈ R M are affinely independent if y 1 −y 0 , . . . , y N −y 0 are linearly independent. Note that in this case y 0 − y n , . . . , y N − y n are linearly independent for any n = 1, 2, . . . , N. Also, linearly independent vectors are affinely independent. Proof. Consider the -th linkage class. Let N be the number of complexes in this linkage class, and let S be the stoichiometric space belonging to this linkage class. Let δ := N − 1 − dim S be the deficiency of the -th linkage class. Finally, let S be the full stoichiometric space.
The deficiency of the entire network can be related to the deficiencies of the linkage classes: For weakly reversible realizations, the deficiency δ is minimized if the complexes are precisely the exponents of monomials [7] because each addition complex (not appearing as exponents of monomials) increases δ by 1. It has been shown that any deficiency zero realization with one terminal strong linkage class is unique if the set of complexes is fixed [11]. Combining these two facts, we have the following Proposition:  Nonetheless, Proposition 1 naturally leads us to the following conjecture: Conjecture 2. Any weakly reversible, deficiency zero reaction network has a unique realization.
Throughout this work, we have seen over and over again that realizability (more precisely, identifying the network structure) is an underdetermined problem. Is it possible to use this flexibility in network structure to get more out of a kinetic differential equation?
An idea comes to mind; we will formulate a conjecture following an example.
Example 9. Consider the mass action system for any rate coefficients k 1 , k 2 , k 3 , k 4 > 0. This system has a weakly reversible, deficiency zero realization: The induced kinetic differential equation of both systems iṡ In particular, the latter mass action system is complex-balanced for all choice of rate coefficients. This implies that the kinetic differential equation enjoys all the dynamical properties of complex balanced systems, although we could have been looking at the network that is not weakly reversible.
Conjecture 3. Consider a reaction network N 1 , and assume that there exists a weakly reversible reaction network N 2 of deficiency δ such that any differential equation induced by N 1 and positive k values can also be induced by N 2 and positive k values. Then, the set of rate coefficient values of N 1 for which N 1 induces systems that have complex balanced realizations contains an algebraic variety of codimension not larger than δ.

Discussion
We have been talking about finding realizations (with certain network properties) of kinetic differential equations. Starting with a system of polynomial differential equations (possibly with symbolic parameters), one could try to find a reaction graph (Feinberg-Horn-Jackson reaction graph) that induces the differential equation under mass action kinetics. After fixing a set of complexes, algorithms exist to find realizations (Section 5.1). Furthermore, conditions on the resulting reaction graph (e.g., weakly reversible, complex balancing, omit certain reactions) can be incorporated into the algorithms. † We also considered operations allowed on a reaction graph that will preserve the kinetic differential equation. For example, in certain cases reactions can be added, and complexes can be removed. Our discussion in Section 5.2 is far from complete; it would be interesting to further explore the allowed operations while preserving dynamical equivalence.
In Section 6, we explored the relationship between weakly reversible and complex balancing realizations. Hidden inside the proof of the main theorem in that section is the scaling equation (6.2). We believe there is information hidden about the manifold in parameter space that decides whether a reaction network is complex balanced or not.
Finally, in our discussion of deficiency zero realizations and uniqueness, we posed two conjectures. The first is that weakly reversible, deficiency zero networks have unique realization. The latter is more complicated; it asks for what rate coefficients is the reaction network dynamically equivalent to complex balancing.
As the reader can see, a theoretical understanding of the relationships between network structure, the sign structure of the differential equation, and dynamical property is lacking. We hope to see further studies exploring the realizability problem. Development, and Innovation Office, Hungary (SNN 125739). Gheorghe Craciun and Polly Y. Yu have been supported by National Science Foundation grant DMS-1816238. is given with k 1 , k 2 > 0, then at least four complexes are needed to realize the equation. How do we prove this?
Two complexes are not enough With the notations of Theorem 3 we have As Y is invertible, v = Y −1 0 = 0 should hold. However, YB = B = Z cannot hold, because Z is not a compartmental matrix.
Three complexes are not enough either Include ( α 1 α 2 ) among the complex vectors with unknown nonnegative integers α 1 and α 2 so that this vector is different from the other two we have as yet included. Then our task is to find nonnegative rate coefficients k 21 , k 31 , k 12 , k 32 , k 13 , k 23 (at least two of them positive) so that First of all, one should have zero coefficients before x α 1 y α 2 , thus α 1 = k 31 k 31 +k 32 α 2 = k 32 k 31 +k 32 , should hold, therefore 1 − α 1 = α 2 1 − α 2 = α 1 . Equating the coefficients of x and y here with those of the right hand side gives an equation without solution for k 12 , . . . , α 2 so that α 1 and α 2 are nonnegative integers, and among the nonnegative numbers k nq at least two are positive. The equations are as follows.
Note that this realization has the minimal number of reaction steps, (one reaction step is not enough), and is of deficiency zero, as well.
Summing up: here the right hand side of (8.1) contains two monomials, still there exists no realization with two complexes, more are needed. Furthermore, three complexes are enough if we allow complexes in (R + 0 ) M instead of N M 0 . In particular, α 1 = 4 5 , α 2 = 3 5 is an appropriate choice. In the pioneering paper [20] this generalization to non-negative real coefficients for complexes is allowed, although here we restrict ourselves to nonnegative integer coefficients.
Example 12. Invertibility of the complex matrix Y is not sufficient to have a realization with the exponents of the monomials as complex vectors, as the following example shows.
x = x 3 − xy 2 ,ẏ = 2x 3 implies that Y = 3 1 0 2 , Z = 1 −1 2 0 and the (unique) solution to with some negative numbers. We are forced to formulate a problem: Characterize those matrices with nonnegative (integer) components the inverse of which multiplied with a matrix without negative cross effect gives a compartmental matrix. As this problem again requires finding the feasible solutions of a linear programming (more precisely, MIL) problem, there is not too much hope to find an explicit symbolic solution, but at he same time one can have numerical solution(s) using LP-solvers.
Example 13. Consider the polynomial equatioṅ x = k 1 x 2 y 5 z − k −1 x 3 y 4 zẏ = −k 1 x 2 y 5 z + k −1 x 3 y 4 zż = 0 (8.6) with k 1 , k −1 > 0. The canonical representation of this equation is a weakly reversible reaction network with N = 4 complexes, L = 1 linkage class and a S = 2dimensional stoichiometric space, thus it has a deficiency 1 = N − L − S . Direct calculation shows that the only positive stationary point is x * = x 0 + y 0 1 + K y * = K(x 0 + y 0 ) 1 + K z * = z 0 , however, its existence directly follows from weak reversibility [3]. Example 14. The canonical realization of the differential equatioṅ x = −4x + 2y + 3zẏ = x − 5y + z + 1ż = 2x + 3y − 4z (8.9) is the irreversible, deficiency three reaction network consisting of a single linkage class of Fig. 2. However, it is "obviously" the induced kinetic differential equation of the zero deficiency weakly reversible reaction network in Fig. 3. (Here again, the unique positive stationary point can explicitely be determined.) Let us see a nonlinear example.
Example 15. The canonical realization of the differential equatioṅ x = −12x 3 + 6y 2 ,ẏ = 2x 3 − 4y 2 + 2 (8.10) is the irreversible, deficiency 3 reaction network consisting of two linkage classes However, it is (much less) "obviously" the induced kinetic differential equation of the weakly reversible, zero deficiency reaction network in Fig. 4 consisting of a single linkage class. To show this Contrary to what the title suggests, here we show a network with thee complexes inducing a differential equation having two terms on the right hand side.
Example 17. The induced kinetic differential equation of the reaction network is (8.6) which we repeat here: x = k 1 x 2 y 5 z − k −1 x 3 y 4 zẏ = −k 1 x 2 y 5 z + k −1 x 3 y 4 zż = 0 (8.13) This is of the formẋ therefore an inducing reaction network is just the one in Eq. (8.12).

Mass conserving realizations
The reaction (2.1) is mass conserving, if there exists a vector ∈ (R + ) M such that for all r = 1, 2, . . . , R M m=1 α(m, r) m = M m=1 β(m, r) m (8.15) holds. Note that to ensure mass conservation the existence of a positive vector in the left kernel of the span of the range of the right hand side is necessary but not sufficient, [17, p. 89]. Consider again the differential equationẋ = Zx Y . (8.16) According to the definition (8.15) the existence of a vector ∈ (R + ) M such that is a necessary condition of mass conservation. However, it is only sufficient with some restrictions.