Simplification of Reaction Networks , Confluence and Elementary Modes †

Reaction networks can be simplified by eliminating linear intermediate species in partial steady states. In this paper, we study the question whether this rewrite procedure is confluent, so that for any given reaction network with kinetic constraints, a unique normal form will be obtained independently of the elimination order. We first show that confluence fails for the elimination of intermediates even without kinetics, if “dependent reactions” introduced by the simplification are not removed. This leads us to revising the simplification algorithm into a variant of the double description method for computing elementary modes, so that it keeps track of kinetic information. Folklore results on elementary modes then imply the confluence of the revised simplification algorithm with respect to the network structure, i.e., the structure of fully simplified networks is unique. We show however that the kinetic rates assigned to the reactions may not be unique, and provide a biological example where two different simplified networks can be obtained. Finally, we give a criterion on the structure of the initial network that is sufficient to guarantee the confluence of both the structure and the kinetic rates.


Introduction
Chemical reaction networks are widely used in systems biology for modeling the dynamics of biochemical molecular systems [1][2][3][4].A chemical reaction network has a graph structure that can be identified with an (unmarked) Petri net [5].Beside of this, it assigns to each of its reactions a kinetic rate that models the reaction's speed.Chemical reaction networks can either be given a deterministic semantics in terms of ordinary differential equations (ODEs), which describes the evolution of the average concentrations of the species of the network over time, or a stochastic semantics in terms of continuous time Markov chains, which defines the evolution of molecule distributions of the different species over time.In this paper, we focus on the deterministic semantics.
Reaction networks modeling molecular biological systems-see, e.g., the examples in the BioModels database [6]-may become very large if modeling sufficient details.Therefore, biologists like to abstract whole subnetworks into single black-box reactions, usually in an adhoc manner that ignores kinetic information [7,8].The absence or loss of kinetic information, however, limits the applicability of formal analysis techniques.Therefore, much effort has been spent on simplification methods for reaction networks that preserve the kinetic information (see [9] for an overview).
The classical example for a structural simplification method is Michaelis-Menten's reduction of enzymatic networks with mass-action kinetics [10].It removes the intermediate species-the complex C and enzyme E-under the assumption that their concentrations C(t) and E(t) are quasi steady, i.e., approximately constant for all time points t after a short initial phase.Segel [11] shows how to infer Michaelis-Menten's simplification from the assumptions that C(t) is constant and that the conservation law C(t) + E(t) = E(0) holds.This is equivalent to our exact steadiness assumption for both C(t) and E(t).
The ODEs for C inferred from this network jointly with exact steady state assumptions for C and E entail that the concentration of substrate S must be constant too, even if the network is used in a bigger context where the intermediate C is neither produced nor consumed.In the literature, this consequence is usually mentioned but ignored when considering the production rate of product P as a function of the concentration of S for the enzymatic network in isolation (see e.g., [12]).This oversimplification can be avoided when studying the enzymatic network in the context of a larger network.For instance, the steady state assumptions for C, E, and thus S can be satisfied in the context of the reaction network with the reaction ∅ k 4 − → S which produces S with constant speed k 4 , and the reaction P k 5 P − − → ∅ which degrades P with mass-action kinetics with rate constant k 5 .In this context, the concentration of P will saturate quickly under exact steady state assumptions for C, E, and thus S, as illustrated in Figure 1, while in other contexts it may grow without bound or even oscillate.The Michaelis-Menten simplification of the enzymatic network indeed preserves the dynamics of a network in any context which does not produce nor consume the intermediates E and C, under the assumption that E and C are exactly in steady state with respect to the network in the context.Whether exact steady state assumptions are realistic is an interesting question since the concentrations may be at most close to steady in practice.In the literature it has been argued that the Michaelis-Menten simplification yields a good approximation under appropriate conditions [11,13,14], which typically depend on the context.Whether such properties can be extended to more general simplification methods as developed in the present article is an interesting question but out of the scope of the paper.
Alternatively, much work was spent on simplifying the ODEs inferred from a given reaction network [15,16], rather than the reaction network by itself.Indeed, any structural simplification method on the network level, that preserves the kinetic information with respect to the deterministic semantics, must induce a reduction method on the ODE level.The opposite must not be true, since some ODEs may not be derivable from any reaction network or may be inferred from many different ones [17].Furthermore, it is not clear what it could mean for an ODE simplification method to be contextual.Therefore, ODE simplification alone cannot be understood as a simplification of biological systems.
A general structural simplification algorithm for reaction networks with deterministic semantics was first presented by Radulescu et al.They proposed yet another method [18] for simplifying reaction networks with kinetic expressions in partial steady states.Their method assumes the same linearity restriction considered in this paper, preserves exactly the deterministic semantics, but uses different algorithmic techniques.Their simplification algorithm is based on a graph of intermediate species.It computes cycles for simplifying the network structure rather than on elementary modes, and spanning trees for simplifying the kinetic expressions.A set of intermediate species is eliminated in one step, leading to a unique result, that is included in the results found with the algorithm of the present paper.We have not understood yet what distinguishes this result from the others obtained with our algorithm; a clarification of this point might shed light on the relationship between the two methods.In the same paper, the authors also observe that applying the method iteratively to intermediates one by one leads to different results even with different structure.The reason is that dependency elimination is lost in this manner.
A purely structural simplification algorithm method for reaction networks without kinetic rates was proposed in [19].The method allows to remove some intermediate species by combining the reaction producing and consuming them.For instance, one can simplify the network with the following two reactions on the left into the single reaction on the right, by removing the intermediate species B: Since no partial steady state assumptions can be imposed in a kinetics free framework, the intermediate elimination rules need some further restrictions.Given these, the simplification steps were shown correct with respect to the attractor semantics.contextual equivalence relation was obtained by instantiating the general framework for observational program semantics from [20].Rather than being based on termination as observable for concurrent programs, it relies on the asymptotic behaviours of the networks represented by the terminal connected components, which are often called attractors.

Outline
We first recall some basic notions on confluence, multisets, and commutative semigroups in Section 2. In Section 3 we recall the basics on reaction networks without kinetics and elementary flux modes.In Section 4, we present the rewrite rules for intermediate elimination, illustrate the failure of confluence, and propose a rewrite rule for eliminating dependent reactions, which however turns out to be non-confluent on its own.In Section 5 we present the refined algorithm in the case without kinetics based on the notion of flux networks for representing reaction networks, and prove its confluence by reduction to a folklore result on elementary flux modes.In Section 6, we introduce reaction networks with kinetic expressions, and extend them with kinetic constraints.In Section 7, we lift the revised algorithm to constrained flux networks with kinetics.In Section 8 we present a linearity restriction, that is preserved by reductions, and thus structurally confluent.In Section 9, we present a counter example that shows that full confluence is still not achieved, and present a further syntactic restriction based on elementary modes avoiding this problem.Section 10 provides a biological example of non-confluence with kinetics.Section 11 studies the relation between the simplification and the underlying ODEs simplification.Finally, we conclude in Section 12.

Preliminaries
We recall basic notions on confluence of binary relations, on multisets, and more general commutative semigroups.We will denote the set of all natural numbers including 0 by N and the set of integers by Z.

Confluence Notions
We recall the main confluence notions and their relationships from the literature.Let (S, ∼) be a set with an equivalence relation and → ⊆ S × S a binary relation.In most cases, ∼ will be chosen as the equality relation of the set S, which is The confluence notions are illustrated by the diagrams in Figure 2. Clearly, a confluence of relation → is confluent if its reflexive transitive closure → * commutes with itself.It is also obvious that local confluence implies confluence, and well known that the converse does not hold.In this paper, we will always use binary relations that are terminating, i.e., for any s ∈ S there exists a k ∈ N such that {s | s → k s } = ∅, i.e., the length k of sequences of reduction steps starting with s is bounded.It is well known that locally confluent and terminating relations are confluent (Newman's lemma).The conditions that have to be satisfied by simulations are illustrated by the diagrams in Figure 3.

Multisets
Let R be a finite set.A multiset M with elements in R is a function M : R → N.For any r ∈ R we call M(r) the number of occurrences of r in M. We say that r is a member of multiset M and write r ∈ M if M(r) = 0. We denote by M R the set of all multisets (over R), and will simply write M if the set R is clear from the context.
Given numbers k, n 1 , . . ., n k ∈ N and a subset {r 1 , . . ., r k } ⊆ R with k different elements, we denote by i=1 n i r i the multiset that for any 1 ≤ i ≤ k contains M(r i ) = n i occurrences of r i and M(r) = 0 occurrences of all other elements in R.
The sum of two multisets M 1 + M M 2 is the multiset M that satisfies M(r) = M 1 (r) + N M 2 (r) for all r ∈ R. The empty multiset 0 M is the function that maps all elements of R to 0. The algebra of multisets (M, + M , 0 M ) over a given set R is a commutative semigroup with a neutral element.
It should be noticed that our notation may give rise to some ambiguities, since we will also write + for the addition of natural numbers instead of + N .This may be problematics if R = N.In this case, the notation introduced below we will permit us to write (n for sums of natural numbers.

Commutative Semigroups
Let (G, + G , 0 G ) and (F , + F , 0 F ) be two semigroups with neutral element.Beside of the algebras of multisets (depending on the choice of R) we are interested in the algebra of vectors of naturals determined by the values of h on singleton multisets in M R via the equation: Given a homomorphism h : M R → F , we define the interpretation M F = h(M) for all multisets M ∈ M R .Clearly, the interpretation depends on the homomorphism h, even though only its co-domain F appears in our notation.This works smoothly since there will never be any ambiguity about the homomophism that is chosen.If R = F , then we use the homomorphism eval F : M F → F with eval F (1 f ) = f for all elements f ∈ F .In this case, any multiset and thus by the above equation: We will also use this notation in order to distinguish the operator + of multisets in M N from the operator + of natural numbers in N which we overloaded (as stated earlier).For instance, if n, m ∈ N, then (2n + 5m) M N is a multiset of natural numbers while (2n + 5m) N is a natural number.Note also that different multisets may have the same interpretation.For instance if n = 3 and m = 4, then (2n + 5m) N = 26 = (n2 + m5) N where we use eval N as homomorphism while (2n + 5m) M N = (n2 + m5) M N where we use id M N as homomorphism.
For any subset G ⊆ G of a semigroup, we can define the (positive integer convex) cone of G, as the set of all positive integer linear combinations of elements of G: Here we use eval G as homomorphism.

Reaction Networks without Kinetics
Let Spec be a finite set of species that is totally ordered.A (chemical) solution with species in Spec is a multiset of species s : Spec → N. A (chemical) reaction with species in Spec is a function r : Spec → Z, which assigns to each species A the stoichiometry of A in r.A chemical reaction r consumes the chemical solution Cons r = −r |{A∈Spec|r(A)<0} and produces the chemical solution Prod r = r |{A∈Spec|r(A)>0} .Clearly r(A) = Prod r (A) − Cons r (A) for all species A, while Cons r and Prod r are disjoint multisets in chemical reactions r (since their definition is based on stoichiometries).
We will freely identify a reaction r with the pair of chemical solutions consumed and produced by r.We will denote such pairs as Cons r AProd r .For instance, B + 2C AA is the chemical reaction r with r(A) = 1, r(B) = −1, and r(C) = −2.Note also that we do not consider 2A + B A3A + 2C as a chemical reaction, since the species A belongs to the chemical solutions on both sides.When removing 2A on both sides, we obtain a chemical reaction B AA + 2C.The rewrite relation of a chemical reaction r contains all pairs of chemical solutions (s, s ) such that s (A) = s(A) + N r(A) for all species A. Definition 3. A reaction network (without kinetics) over Spec is a finite set of chemical reactions over Spec, with a total order.
To any reaction network N with total order < we assign a unique vector of reactions r = (r 1 , . . ., r n ) such that N = {r 1 , . . ., r n } and r 1 < . . .< r n .Conversely, for any tuple of distinct reactions r = (r 1 , . . ., r n ), we write N r for the reaction network {r 1 , . . ., r n } with the total order r 1 < . . .< r n .
Any reaction network can be represented by a bipartite graph as for a a Petri net, with a node for each species and a node of a different type for each reaction.We will draw species nodes with ovals and reaction nodes with squares.An arrow labeled by k from the node of a species A to the node of a reaction r means that A is consumed k times by r, i.e., r(A) = −k.Conversely, an arrow with label k from the node of a reaction r to the node of a species A means that A is produced k times by r, i.e., r(A) = k.We will freely omit the labels k = 1.
Example 1.Consider the reaction network presented in Figure 4.It has m = 2 species Spec = {X, Y} and n = 4 reactions {r 1 , . . ., r 4 } in that order.Reaction r 1 produces two molecules of species X out of nothing, reaction r 2 transforms an X into a molecule Y, while r 3 transforms a molecule Y back into a molecule X. Reaction r 4 degrades a molecule X.
The set of chemical reactions defines an algebra (R, + R , 0 R ) where 0 R is the empty reaction A, and + R is the addition of integer valued functions on Spec.Note that s As + R s As = 0 R for any two disjoint chemical solutions s and s .By interpretation in this algebra (that is using the identity homomorphism), we can evaluate each multiset of chemical reactions M as a chemical reaction M R itself, as shown in Section 2.2.Definition 4.An invariant of a reaction network N without kinetics is a multiset M of reactions of N such that M R = 0 R .We denote the set of all invariants of N by inv(N).
The reaction network in Figure 4 has the set of invariants {(n We next relate the notion of invariants of a reaction network to the kernel of its stoichiometry matrix.
A reaction network and the associated graph and stoichiometry matrix.

Stoichiometry Matrices
The stoichiometry information of a reaction network is usually collected in its stoichiometry matrix.For this we consider a set of species Spec = {A 1 , . . ., A m } and a reaction network N = {r 1 , . . ., r n }, such that both sets are totally ordered by the indices of their elements.
The stoichiometry matrix S of N is the m × n matrix of integers, such that the entry of S at row i and column j is equal to r j (A i ) for all 1 ≤ i ≤ m and 1 ≤ j ≤ n.Note that reaction r j contributes in the j th column, while species A j contributes the j's row of S. For instance, the stoichiometry matrix of the reaction network in Figure 4 is given on the right.
It can now be noticed that, for any vector v = (n 1 , . . ., n n ) of natural numbers, the multiset n 1 r 1 + ... + n n r n is an invariant of reaction network N if and only if its stoichiometry matrix satisfies Sv = 0, i.e., if v belongs to the kernel of the stoichiometry matrix.Therefore, we define the (positive integer) kernel of a matrix S by:

Elementary Modes
The support of a vector v = (n 1 , . . ., n k ) is the subset of indices i such that n i is non-null, i.e., supp(v) = {i ∈ {1, . . ., k} | n i = 0}.Definition 5.An elementary mode of an m × n matrix S over Z is a vector v ∈ ker + (S) \ {0 N n } such that: v is on an extreme ray: there exists no v ∈ ker + (S) \ {0 N n } such that supp(v ) supp(v), and v is factorised: there exists no v ∈ ker + (S) such that v = kv for some natural number k ≥ 2.
The condition v ∈ ker + (S) means that an elementary mode must be a (positive integer) steady state of S. Geometrically, the set of all positive integer steady states forms a pointed cone, that is generated by convex combinations of its extreme rays.The first condition states that any elementary flux mode v must belong to some extreme ray of the cone.The second condition requires that an elementary mode is maximally factorised, i.e., it is the vector on the extreme ray with the smallest norm.
Theorem 1 (Folklore [21]).Let S be an m × n matrix of integers.Then the set E of all elementary modes of S has finite cardinality and satisfies ker + (S) = cone(E).
The intuition is ker + (S) is a cone with a finite number of extreme rays, so that these extreme ray generated the cone.The set of elementary modes E contains exactly one point on each of the extreme rays of ker + (S).Therefore, ker + (S) = cone(E), i.e., the set of elementary modes is a finite generator of ker + (S).
Let us point out two differences between the definition of elementary mode considered here and in [21].First, we added condition 2. Without this condition, any multiple of an elementary mode would be an elementary mode, so that there would be infinitely many.The double-description method as recalled there, however, computes the set of elementary modes in the above sense, so this difference is minor.Second, note that [21] considers a slightly more general problem, where some of the coordinates of v may be negative.This corresponds to the addition of reversible reactions that we do not consider in the present paper.

Elementary Flux Modes
We next lift the concept of elementary modes from matrices to reaction networks, via the stoichiometry matrix.Given a vector of reactions r = (r 1 , . . ., r n ) and a vector v = (n 1 , . . ., n n ) of natural numbers we define the multiset of reactions vr and the corresponding reaction r v as follows: Definition 6.An elementary flux mode of a reaction network N = N r is a multiset of reactions vr such that the vector v is an elementary mode of the stoichiometry matrix of N.
The kernel condition v ∈ ker + (S) of elementary modes v yields that any elementary flux mode vr satisfies r v = 0 R , i.e., the reaction defined by the elementary flux mode must be empty.For instance, reconsider the reaction network in Example 1 with m = 2 species and n = 4 reactions r = (r 1 , r 2 , r 3 , r 4 ) in that order.Its stoichiometry matrix has two elementary modes: the vectors v 1 = (1, 0, 0, 2) and v 2 = (0, 1, 1, 0).The corresponding elementary flux modes are the multisets of reactions v 1 r = r 1 + 2r 4 and v 2 r = r 2 + r 3 illustrated in Figure 5 by the arrows coloured in apricot and aquamarine respectively.First consider the multiset r 1 + 2r 4 : the first reaction r 1 produces 2X which are then degraded by 2r 4 .So the reaction r v 1 = (r 1 + 2r 4 ) R = 0 R is indeed empty.Consider now the multiset of reactions r 2 + r 3 : its first reaction r 2 transforms X to Y and its second reaction r 3 does the inverse.Thus, r v 2 = (r 2 + r 3 ) R = 0 R is the empty reaction too.The intuition is that applying to a chemical solution at the same time all reactions of an elementary flux mode with their multiplicities does not have any effect.It should be noticed that the vector v = (1, 1, 1, 2) is also a solution of the steady state equation Sv = 0, and thus the multiset of reactions r 1 + r 2 + r 3 + 2r 4 is also an invariant of the example network.It is the multiset sum of two elementary flux modes v 1 r + M v 2 r which is also equal to (v 1 + N 4 v 2 )r.

Simplifying Reaction Networks without Kinetics
We study the question whether the step-by-step intermediate elimination relation proposed in [22] is confluent in the case of reaction networks without kinetics.We present a counter example against the confluence and illustrate the reason for this problem.

Intermediate Elimination
Let I ⊆ Spec be a finite set of species that we will call intermediate species or intermediates for short.The simplification procedure will remove all intermediates from a given reaction network, step-by-step and in arbitrary order.
Our objective is to remove an intermediate X ∈ I from a network N by merging any pair of reactions of N, a reaction r that produces X and another reaction r that consumes it.This is done by the (INTER) rule in Figure 6, and is based on the merge operation r X r which returns a linear combination of r and r and thus of the reactions in the initial network: Since r produces r(X) molecules X while r consumes r (X) molecules X, we have (r X r )(X) = 0. Therefore, X is not present in the solutions consumed and produced by reaction r X r .
(INTER) In Example 2 below, we will denote vectors (n 1 , . . ., n n ) of natural numbers by 1 n 1 . . .n n n , while freely omitting components i n i with n i = 0 and simplifying component j 1 to j.For instance if r = (r 1 , . . ., r 4 ), we can write r 14 2 instead of r (1,0,0,2) .Example 2. We consider the network N in Figure 7 with species Spec = {A, B, X, Y} and reaction vector r = (r 1 , . . ., r 4 ).We consider the elimination of the intermediates in I = {X, Y} in both possible orders.On the top, we first eliminate the intermediate species X from N, obtaining network N X .We have to combine reaction r 1 producing 2 X molecules with reaction r 2 which consumes 1 X molecule.We obtain the reaction r 1 X r 2 = r 12 2 = (r 1 + 2r 2 ) R that transforms one A molecule into 2 Y molecules.We proceed in the same way with the other 3 pairs of reactions that produce and consume X.Then, we can remove the intermediate species Y from network N X and obtain the network N XY in the top right.Note that we keep empty reactions such as r 23  It turns out that N XY and N YX differ in that the former contains the reaction r 12 2 3 2 4 2 in addition to the reactions r 14 2 and r 23 shared by both networks.

Eliminating Dependent Reactions
Example 2 shows that intermediate elimination with the (INTER) rule alone is not confluent, given that it may produce two different networks that cannot be simplified any further, N XY and N YX , depending on whether we first eliminate the intermediate X or the intermediate Y.The reaction network N XY contains an additional reaction, which is a linear combination of two other reactions: In order to solve this non-confluence problem, we propose the new simplification rule (DEP) in Figure 6.It eliminates a reaction that is a positive linear combination of other reactions of the network, i.e., some reaction Unfortunately, the simplification relation with rules (INTER) and (DEP) is still not confluent.The problem is that even applying rule (DEP) alone fails to be confluent as shown by the following counter example.
Example 3. Consider the network N in Figure 8 in the absence of intermediates, i.e., where I = ∅.There are two ways of applying rule (DEP) to this network, since r 4 = (r 1 + 2r 2 ) R and r 2 = (r 3 + r 4 ) R .We can thus either eliminate r 4 leading to N r 4 or r 2 leading to N r 2 .The two results are different even though they contain no more dependencies.
Network N r 4 This example shows that general dependency elimination cannot be done in a confluent manner.On the other hand, what we need in order to solve the confluence problem for intermediate elimination as illustrated in Figure 7, is a little more restricted: it is sufficient to remove those dependent reactions that were introduced by intermediate elimination.Such dependencies can be identified from the vectors of natural numbers that we used to name the reactions.In the example, we have r 12 2 3 2 4 2 = (r 14 2 + 2r 23 ) R , so the dependency of this reaction follows from the dependency of the vectors 12 2 3 2 4 2 = (( 142 ) + 2( 23)) N 4 .

Simplifying Flux Networks
We next introduce vector representations of reaction networks without kinetics, called flux networks, and show that the simplification of such representations can indeed be done in a confluent manner.
For the reminder of this section, we fix an n-tuple r of distinct reactions and a subset of species I ⊆ Spec.

Vector Representations of Reaction Networks
The objective is to simplify the initial reaction network N r by removing the intermediates from I. The iterative elimination of intermediate species generates a sequence of networks with reactions in The idea is now to use the vectors v ∈ N n as representations of reactions r v .These vectors will tell us about the provenance of the reaction obtained when simplifying the network N r .
The mapping of vectors v ∈ N n to reactions r v ∈ R is a homomorphism between commutative semigroups, whose image is cone(N r ).It should be noticed, however, that it is not an isomorphism since any element of ker + (S) will be mapped to 0 R , where S is the stoichiometry matrix of N r .Therefore, it makes a difference whether we will work with vectors in N n representing a reaction or with the reactions itself.Intuitively, the difference is that we know where the reaction does come from.Definition 7.An n-ary flux network V is a finite subset of vectors in N n that is totally ordered.
Any n-ary flux network V defines a reaction network r V = {r v | v ∈ V}, that we call the reaction network represented by V.The total order of the reactions in network r V is the one induced by the total order of V.

Simplification Rules
Let I ⊆ Spec be a finite set of species that we call intermediates.In Figure 9, we rewrite the simplification rules (F-INTER) and (F-DEP) so that they apply to flux networks.For this we have to lift the merge operation from reactions to vectors that represent them.For any v 1 , v 2 ∈ N n we define: In the rule for the dependency elimination, we now use a notation for linear combinations of vectors in N n .Given a vector v = (v 1 , . . ., v k ) of vectors in N n and a vector v = (n 1 , . . ., n k ) of natural numbers we define: The counter example for the non-confluence of dependency elimination can no more be applied in this way, since rule (F-DEP) is not based on the dependency of the reactions as with (DEP) but on the dependencies of the vectors that define the reactions. (F-INTER)

Factorization
The simplification relation with axioms (F-INTER) and (F-DEP) is still not confluent, as shown in Example 4.
Example 4. We consider the vector of initial reactions r = (r 1 , r 2 , r 3 ) of the network N in Figure 10.Let I = {X, Y, Z} be the set of intermediate species.Note that N = r V 3 where V 3 = {(1, 0, 0), (0, 1, 0), (0, 0, 1)} ⊆ N 3 is the flux network to which we apply the simplification algorithm.If we remove the species X first from V 3 , we obtain a flux network representing the reaction network N X , and from that we get a flux network representing N XYZ by eliminating Y and Z (in any order).This flux network has only one flux vector which is 1 2 2 2 3 2 = (2(123)) N 3 .If we remove Y first we obtain a flux network representing N Y , and from that a flux network representing N YXZ by removing X and Z.The latter flux network is the singleton with the flux vector 123.
Network What is needed is a rule for the factorization of scalar multiples (kv) N n of vector v.This is done by the rule (F-FACT) in Figure 9, which, in the previous example, allows to simplify N XYZ into N YXZ .
We first note a consequence of the Folklore Theorem 1 and the following Lemma that is equally well known.Lemma 2. Let v, v ∈ N n be two elementary modes of the same matrix S.
be such that n i /n i is maximal.Without loss of generality we can assume that n i /n i ≥ 1 since otherwise, we can exchange v and v .Consider the vector of integers w = n i v − n i v.For any j ∈ supp(v) we have: Therefore, w ∈ N n , and thus w ∈ ker + (S).Furthermore, i / ∈ supp(w) so that supp(w) supp(v).Since v is an elementary mode of S this implies that w = 0, and thus n i v = n i v. Without loss of generality we can assume that n i and n i have no common prime factors.If n i = n i = 1 we are done.Otherwise n i ≥ 2 since n i /n i ≥ 1.Thus v can be factorized by n i , contradiction.

Corollary 1.
Let S be an m × n matrix of integers and E ⊆ N n the set of all elementary modes of S. For any set E ⊆ N n such that cone(E ) = ker + (S), if E is irreducible by F-FACT and F-DEP then E = E.

Proof. By Theorem 1, we have cone(E) = ker + (S) = cone(E ).
We first show that and v i are positive, it follows that supp(v i ) ⊆ supp(v) for all 1 ≤ i ≤ k.Consider i = 1.Since v is an elementary mode and v 1 ∈ ker + (S), this implies that supp(v 1 ) = supp(v).Since v 1 is factorized, and a member of ker + (S) with minimal support, it is also an elementary mode of S. Lemma 2 thus implies that v 1 = v, and so v ∈ E .(It also follows that k = n 1 = 1).
We next show that

Proving Confluence via Elementary Modes
Given a tuple of initial reactions r of size n and a set of intermediates I ⊆ Spec as parameters, we obtain a simplification relation on flux networks: We now show that this relation is confluent for all possible choices of the parameters.The proof is by reduction to the Corollary 1 of the folklore Theorem 1 on elementary modes.We start with an fundamental property of the diamond operator X , that we formulate in a sufficiently general manner so that is can be reused later on.
) be a commutative semi-ring and h : N n → G a semi-group homomorphism with respect to addition.Given a tuple (v 1 , . . ., v k ) of vectors in N n , a tuple (g 1 , . . ., g k ) of elements of G, and a species X ∈ Spec, we define: It then holds that: Proof.We use some elementary rules of commutative semi-rings to distribute and factorize the sums contained in the definition of the diamond: Our next objective is to show that the simplification preserves the invariants, when relativised to r.For any flux network V, we therefore define the set of relatived invariants of V as follows: For V n = {(1, 0, . . ., 0), . . ., (0, . . ., 0, 1)} ⊆ N n with the vectors ordered in the way they are enumerated, note that we have r V n = N r and inv r (V n ) = inv(N r ).We next show that such relativised invariants are preserved by the simplification of flux networks.
Proof.We assume V F V and first show the inclusion inv r (V) ⊆ inv r (V ). Let Otherwise, we can assume without loss of generality that v 1 = v v with v and v as in rule (F-DEP).Suppose that these have the forms v = (m 1 , . . ., m l ) and v = (w 1 , . . ., w l ).Since This yields (n We can assume without loss of generality that n i = 0 for all 1 ≤ i ≤ k.Let P, C, prod, and cons be as introduced in the Diamond Lemma 3, where G = N n , homomorphism h the identity on N n , and g i = n i for all 1 ≤ i ≤ k.The lemma then yields: This multiset is an invariant, since (∑ k i=1 n i r v i ) R = 0 R .It follows that: This implies (∑ k i=1 prod n i v i ) N n r ∈ inv r (V ).Since prod = 0 and since inv r (V ) is closed by factorization with nonzero factors, it follows that ∑ k i=1 n i r v i ∈ inv r (V ) as required.The proof of the inverse inclusion inv r (V) ⊇ inv r (V ) differs in that the Diamond Lemma is not needed.Let We distinguish three cases depending on which rule was applied: And thus, (n Otherwise, we can assume without loss of generality that v 1 = v v with v and v as in the rule.Suppose that these have the forms v = (m 1 , . . ., m l ) and v = (w 1 , . . ., w l ).Since r v v = (m 1 r w 1 + ... + m l r w l ) R , it follows that: This yields (n Case V F-INTER V .Suppose that the intermediate species X ∈ I was eliminated thereby.We recall that ∑ k i=1 n i r v i ∈ inv(r V ).Without loss of generality, we can assume that all elements of V occur exactly once in this sum.Let V = {v 1 , . . ., v l }, By the rule (F-INTER) we have: We start the reminder of the proof with the case where all species are intermediates so that I = Spec.
Proof.Given that V is irreducible by F-INTER , all intermediates species must be eliminated in all reactions of r V .Since I = Spec this implies that all species are eliminated in all reactions of r V , so for all v ∈ r V it follows that , where S is the stoichiometry matrix of N r .
Proof.By Lemmas 4 and 5 we have: Theorem 2. Consider the simplification relation for flux networks F that is parametrised by I = Spec and a tuple of initial reactions r.If V n * F V for some flux network V that is irreducible for F , then V = E, where E is the set of elementary modes of the stoichiometry matrix of N r .

Proof.
From Proposition 1 it follows that cone(V) = ker + (S) where S is the stoichiometry matrix of N r .Furthermore, V is irreducible with respect to F-FACT ∪ F-DEP , so that Corollary 1 implies V = E. Proof.We notice that F is terminating, since (F-INTER) reduces the number of intermediate species X ∈ I for which there exists a vector v such that r v (X) = 0, (F-DEP) reduces the number of vectors in the set, and (F-FACT) reduces the norm of one of the vectors.We first consider the case Spec = I.Let V be such that V n * F V, where F is parametrised by I and a tuple r of initial reactions.Suppose that , where E is the set of elementary modes of the stoichiometry matrix of r V n = N r .
We next reduce the general case where Spec ⊆ I to the case Spec = I.We define r |I by restricting all reactions in the tuple r to I, i.e., if r = (r 1 , . . ., r n ) then r |I = (r 1|I , . . ., r n|I ).We then observe that the relation F with respect to r coincides with the relation F with respect to r |I .Hence the confluence result from the case I = Spec can be applied.
As shown by Theorem 2, the exhaustive simplification of flux networks V with F can be used to compute the set of elementary modes of the stoichiometry matrix of the reaction network r V .
Interestingly, this algorithm is essentially the same as the double description method, as recalled for instance in [21].The correspondence comes from the fact that any reaction network can be identified with its stoichiometry matrix, so that the algorithm can be formulated either for the one or the other representation.Still there is a minor difference between this algorithm and the one in [21].The algorithm presented here is slightly more flexible, in that the rule (F-DEP) can be applied at any stage of the simplification while in the double description method as described in [21], the rule (F-DEP) is applied at the same time as the rule (F-INTER).However, as we have shown with the confluence Theorem 2, this additional freedom in the application order of the rules does not affect the final result.

Reaction Networks with Deterministic Semantics
We now consider reactions with kinetic expressions, and recall some basic definitions.We first define expressions and networks with kinetics.Then we recall how to associate a system of equations to a reaction network.Finally we use this system of equations to define the deterministic semantics of reaction networks.

Kinetic Expressions
We now define a class of kinetic expressions.Their syntax is the same as that of arithmetic expressions, by their semantics is by interpretation as functions of type R + → R, Let Param be a set of parameters of type R + .As set of variables of type R + → R + , we will use the set Spec.A variable A ∈ Spec is intended to represent the temporal evolution of the concentration of A over time.
We define the set of expressions Expr by the terms with the abstract syntax in Figure 11.Expressions describe functions of type R + to R. They are built from species A of type R + to R + , and constant functions defined by parameters k ∈ R + , constants c ∈ R, and expressions e(0), standing for the value of e at time 0. Beside of these, expressions can be constructed by addition, subtraction, multiplication, and division.For convenience, we will use parenthesis (e) whenever the priority of the operators might not be clear.For any species e we denote by pSpec(e) the subset of species that occur properly in e, that is outside of a sub-expression e(0).So for instance pSpec(B = A(0)) = {B}. .
Semantics The semantics of expressions is parametrised by a function β : Param → R + that interprets all parameters as positive real numbers.In order to simplify the notation, we assume that β is fixed, but notice that our simplification algorithms will be correct for any interpretation β.
The value of an expression e α ∈ (R + → R + ) ∪ {⊥} is specified in Figure 11 for any variable assignment α : Vars → (R + → R + ).It may either be a function of type R + → R or undefined ⊥.The latter is necessary for the interpretation of 1/e α , which is defined only if e α (t) = 0 for any time point t.We call an expression e nonnegative if e ≥ 0 is valid, i.e., if for all non-negative assignment α and all time points t ∈ R + , we have e α (t) ≥ 0. Definition 8.A kinetic expression is a nonnegative expression e ∈ Expr.

Constrained Flux Networks
The next objective is to add kinetic expressions to reactions and flux networks.Furthermore, we need to be able to express constraints about these kinetic expressions in order to express partial steady state hypothesis and conservation laws.This will lead us to the notion of constrained flux networks.A reaction with kinetics expressions is a pair r; e where r is a reaction without kinetics and e is a kinetic expression.As before we now use flux reaction to represent a reaction but now with a kinetic expression.Definition 9.An n-ary flux reaction with kinetic expression is a pair v; e composed of a vector v ∈ N n and a kinetic expression e ∈ Expr.Given a tuple of reactions r = (r 1 , . . ., r n ), the flux reaction v; e represents the reaction r v ; e.
The set C of constraints on kinetic functions is defined in Figure 12.A constraint C ∈ C is a conjunction of atomic constraints.The first kind is an equation e = e stating that the expressions e and e must have the same value but different from ⊥.The atomic constraint cst(e) requires that e is a constant function, e = 0 that e may never becomes equal to zero, and e ≥ 0 that e is always non-negative.More formally, we define in Figure 12   Let W = V&C be a constrained flux network.We denote by Expr(W) the set of kinetic expressions e such that v; e ∈ V or such that e occurs in the constraint C. We set:

Systems of Constrained Equations with ODEs
We now recall how to assign systems of equations to constrained flux networks.Note that systems constrained equations may contain both constraints and ordinary differential equations (ODEs) in particular.
The set of systems of constrained equations is defined in Figure 13.They are conjunctions of constraints C and ODEs Ȧ = e where A ∈ Spec and e ∈ Expr.Note that the constraints may subsume the non-differential arithmetic equations e = e .We denote by Spec(E) the set of (free) variables occurring in E, and by Expr(E) the set of expressions contained in E. The denotation of a system of constrained equations E is a value in E α ∈ B ∪ {⊥} as defined in Figure 13.The set of solutions of sol(E) is the set of assignments α : Spec → R + that make E true, i.e., We say that a constrained equation E logically implies another E and write Clearly, E |=| E if and only if E and E logically imply each other, i.e., E |= E and E |= E.

Deterministic Semantics
We assign to any constrained flux network W = V&C a system of constrained equations E(W) in Figure 14.Note that E(W) does depend on the tuple of initial reactions r and the set of intermediates I.The system contains an ODE for any species A stating that the change of the concentration of A is equal to the sum of the rates r v (A)e of the flux reactions v; e ∈ V.The factor r v (A) makes the rate negative if A is consumed and positive if A is produced.It also takes care of the multiplicities of consumption and production.Finally, the constraint C of the constrained flux network is added to the system of constrained equations.Example 5. We consider the flux network for the classical Michaelis-Menten example [10].Its system of constrained equations is then represented in Figure 15.It contains four ODEs and two constant constraints.

Contextual Equivalence
Two constrained flux networks W and W are non-contextually equivalent, denoted W W , if their systems of constrained equations are logically equivalent: We now extend the definition to a contextual equivalence.The idea is that networks can be exchanged with equivalent networks in any context, without affecting the semantics.As contexts, we use flux networks themselves W = V &C .We define the combination of a network W = V&C and the context W as follows: We now assume a set of intermediate species I ⊆ Spec and call a context W compatible if pSpecs(W ) ∩ I = ∅.Definition 12. Two constrained flux networks W and W are (contextually) equivalent if they have the same solutions in any compatible context, that is: We note that the definition of equivalence of constrained flux networks has two parameters: r and I.The equivalence ∼ depends on the tuple of initial reactions r, since the non-contextual equivalence relation relies on the deterministic semantics of constrained flux networks, which in turn depends on r.The equivalence relation ∼ also depends on the set of intermediates I since the notion of compatibility depends on it.
Our simplification algorithm will rewrite constrained flux networks up to logical equivalence of constraints.Therefore, we can hope for confluence only up to logical equivalence.More formally, we defined the similarity relation ∼ = as the least equivalence relation on constrained flux networks that satisfies the following two inference rules for all C, C , V, e, e : The first rule states that expressions that are logically equivalent under the constraints of the constrained flux network can be replaced by each other.The second rule allows to exchange logically equivalent constraints by each other.Similar networks are trivially equivalent: Lemma 6. Similarity W ∼ = W implies contextual equivalence W ∼ W .
Proof.Straightforward from the definitions.

Simplification of Constrained Flux Networks
Our next objective is to simplify constrained flux networks by lifting the confluent simplification algorithm for flux networks to the case with kinetic expressions.This will require to impose partial steady state and linearity restrictions on the constrained flux networks, since otherwise, we would not know how to remove intermediates from the constrained equations assigned to the constrained flux network.

Linear Steadiness of Intermediate Species
The following restriction will allow us to eliminate an intermediate species from the constraint equations of a constrained flux network.Definition 13.We say that a species X ∈ I is linearly steady in a constrained flux network V&C if it satisfies the following four conditions: Partial steady state: the concentration of X is steady, i.e., C |= cst(X).Linear consumption: if a reaction in V consumes X then its kinetic expression is linear in X, that is: if v; e ∈ V such that r v (X) < 0 then C |= e = Xe for some expression e such that X / ∈ pSpecs(e ).Independent production: if a reaction in V produces X then its kinetic expression does not contain X except for subexpressions X(0): for any v; e ∈ V, if r v (X) > 0 then X / ∈ pSpecs(e).Nonzero consumption: the consumption of X is nonzero: Suppose that X is linearly steady in W = V&C.Since X is in partial steady state, we have C |= cst(X) and hence C |= Ẋ = 0.The constrained equations of W thus imply that the production and consumption of X are equal: The linear consumption of X imposes that C |= cons = Xe for some expression e such that X ∈ pSpecs(e).The independent production of X imposes that X ∈ pSpecs(prod).Because of nonzero consumption, we have: where the expression prod e does not contain the species X properly.Therefore, we can eliminate the variable X from the constrained equation E(W) by substituting X by prod e .This give us hope that we can also eliminate linearly steady intermediate species from the constrained flux networks too by adapting the rule (F-INTER) to kinetic expressions.

Simplification
We now lift the simplification rules for flux networks to constrained flux networks.The lifted rules are presented in Figure 16.They define the simplification relation for constrained flux networks: The first rule (C-INTER) eliminates a linearly steady intermediate species X, by merging any pair of reactions, so that the one produces and another consumes X.The rule can be applied only under the hypothesis that the constraints of the network imply that X is linearly steady, and so that X is in partial steady state in particular.It should also be noticed that the conditions on the initial value of X are preserved by the constraint X(0) = X prod cons .As argued above, the linear steadiness of X implies that the latter is equivalent to some other expression that does not contain species X properly.So except for constraints on the initial value X(0), the species X got removed from the constrained flux network.The rule also replaces X by X(0) in all the kinetic expressions of reactions, in which X is used as a modifier, i.e., v; e ∈ V such that r v = 0 and X ∈ pSpecs(e).Furthermore, the same substitution is applied to the constraints of the flux network.The rule (C-MOD) removes an intermediate that is never a reactant or a product of a reaction, and replaces X with its initial value X(0).Then the rule (C-DEP) removes a dependent reaction.In contrast to the case without kinetics, the kinetic expressions of the remaining reactions need to be modified.The last rule (C-SIM) states that simplification is applied modulo similarity of constraint flux reaction networks.
The simplification defined here is sound for the contextual equivalence relation of constrained flux networks: The proof is given in Appendix A. The arguments are direct from the definitions, except that the Diamond Lemma 3 is needed in for C-INTER .

Michaelis-Menten
We illustrate the simplification on the classical Michaelis-Menten example [10].We consider the simplification of a three-step enzymatic scheme with mass-action kinetics into a single reaction with Michaelis-Menten kinetics.In the initial network, MMnet depicted in Figure 17, a substrate S can bind to an enzyme E and form a complex C. The complex can either dissociate back to S and E, or produce a product P, while releasing E. We assume here that the enzyme E and the complex C are intermediate species, i.e., they are at steady-state and cannot interact with the context.Therefore, the intermediate species E and C are linearly steady in this network.
We first look at the elimination of the intermediate C with (C-INTER).To this end, we merge each reaction that produces C (that is, reaction r 1 ) with each reaction that consumes C (reactions r 2 and r 3 ) and obtain the network MMnet C .Thus, merging reactions r 1 and r 2 (resp.r 1 and r 3 ) of MMnet (Figure 17) results in the reaction r 12 (resp.r 13 ) of MMnet C .The simplification also replaces the atomic constraint cst(C) with cst(k 1 SE/(k 2 + k 3 )).Since we also have the constraint S, and the parameters are constant too, we can rewrite cst(k 1 SE/(k 2 + k 3 )) into the similar cst(S).We also add the constraint To remove E before the elimination of C, one would merge r 3 with r 1 , r 2 with r 1 , and obtain the network MMnet E .At this point, in both networks we have an intermediate species that is neither a product nor a reactant of any reaction, but is used as a modifier.We can then remove it with (C-MOD), replacing E with E(0) (resp.C with C(0)).We obtain respectively the networks MMnet CE and MMnet EC .Note that these networks are similar.We can rewrite ), and use this equation to rewrite the kinetic rate.Additionally, we can also use it to transform the kinetic expression of r 12 into the usual one for Michaelis-Menten, using the following transformation.First, we have: We can then rewrite it into: By replacing E(0) with this expression in the kinetic expressions of r 13 in MMnet CE , we obtain after basic rewriting the classical rate: The following diagram illustrates the confluence of the simplifications on these networks.

Preservation of Linear Steadiness
As we have seen in the previous section, to remove an intermediate species X, we need to impose that X is linearly steady.When removing a set of intermediate species I, we then need that any X ∈ I is linearly steady, and moreover than when we remove one intermediate species, the other species in I remain linearly steady.We therefore introduce the following additional conditions, and denote by LinNets the set of networks that satisfy these conditions.We then prove that the set LinNets is stable under the simplification, i.e., that the simplification of a network in LinNets is still a network in LinNets.

LinNets
We first define the new conditions, and present some examples to motivate them.For any flux v; ex, we note Definition 14.We denote by LinNets the set of constrained flux networks W such that W is similar to a constrained flux network V&C and that for all intermediates X ∈ I: 1. Either X is linearly steady in V&C, or X is only a modifier, that is for any v; e ∈ V we have r v (X) = 0. 2. No intermediate species different from X occurs in the kinetic expression of a reaction that consumes X: for any v; e ∈ V, if X ∈ Cons I (r v ), then Specs(e) ∩ I ⊆ {X}. 3. The rate of a reaction that produces X but does not consume an intermediate species does not depend on the concentration of any intermediate species: for any v; e ∈ V, if X ∈ Prod I (r v ) and Cons I (r v ) = ∅, then I ∩ Specs(e) = ∅.4. The total stoichiometry of the intermediate species in the reactant (resp.product) of a reaction is never greater than one: for any v; e ∈ V, |Cons I (r v )| ≤ 1 and |Prod I (r v )| ≤ 1.
Note that, as a consequence of the stoichiometry condition, the sets Cons I (r v ) and Prod I (r v ) are either empty or consist of a single intermediate species.
We illustrate the motivations for these new conditions on the following examples.Let us first consider the case where a reaction consuming X has a kinetic rate that depends on another intermediate (here Y in the kinetic rate of r 2 ), so that Condition 2 is not satisfied.It is illustrated in Figure 18.If we remove Y first by merging r 3 and r 4 , then we compute the expression Y = k 3 k 4 X.We replace Y with this expression in the kinetic rate of r 2 , obtaining a reaction with a non-linear kinetic Similarly, consider a reaction producing X with a kinetic rate that depends on Y, i.e., a network where Condition 3 is not satisfied (Figure 19).If we remove the intermediate Y, the kinetic expression of r 1 becomes We obtain a reaction producing X, with a kinetic expression depending on X.Therefore, the differential equation for X will not have the required form: ), and we cannot compute an expression for X.This kind of situation may also appear as the result of the simplification of reactions where one intermediate has a stoichiometry greater than one, i.e., Condition 4 is not satisfied (Figure 20).In this network, reaction r 3 produces two molecules of Y.If we remove Y, the merging of r 3 and r 4 is a reaction that produces one molecule of X (and one of C), with kinetic expression k 3 X.
Finally, if we have two intermediate species that are both reactants (or both products) in the same reactions (Condition 4 again), then the stoichiometry of one intermediate can become greater than one as a result of the elimination of the other intermediate (Figure 21).If we remove Y, the merging of r 3 and r 4 is a reaction with two molecules of X as reactants.

Stability of LinNets
Now, we prove that LinNets is stable for our simplification.We first consider the following proposition.Proposition 3. Let W, W 0 be reaction networks such that W 0 ∈ LinNets and W 0 * C W. Let v; e ∈ W be a flux that depends on v 1 ; e 1 , . . ., v k ; e k ∈ W. Then:

•
there exists an index i such that Cons I (r The proof of this proposition is quite long and requires some new notions and definitions, and is given in Appendix B.
We now prove that LinNets is stable for the simplification.
Proposition 4. The set of networks LinNets is stable for the simplification, that is if W ∈ LinNets and W W , then W ∈ LinNets.
Proof.If the simplification is done with (C-MOD) or (C-SIM), then the conditions of LinNets are trivially preserved.
Let us assume that the simplification is done with the rule (C-DEP), removing a flux v; e that depends on v i ; e i with coefficient a i .So the simplified network contains the fluxes v i ; e i + a i e.By Proposition 3, there is an i such that Cons I (r v i ) = Cons I (r v ) and Prod I (r v i ) = Prod I (r v ), and for any other The fourth condition on the stoichiometry is trivially preserved by the simplification.For j = i, Cons I (r v j ) = Prod I (r v j ) = ∅ implies that the conditions on the kinetic expressions are directly satisfied.If v i ; e i + a i e consumes X, then since Cons I (r v i ) = Cons I (r v ), the flux v; e consumes X too.Then by induction, e and e i are linear in X, and no other intermediate species occurs in them.Therefore, this is also the case for e i + a i e.
If v i ; e i + a i e produces X without consuming any other intermediate, then it is also the case for v; e.Then by induction, e, e i , and e i + a i e do not depend on the concentration of any intermediate species.Therefore, the kinetic conditions are satisfied, and W ∈ LinNets.
Finally, consider the case of a rule (C-INTER) applied on a species X.We denote by prod and cons the expressions defined in the rule.Since W ∈ LinNets, note that for any Z ∈ I, we have Z / ∈ Vars(cons).
Let v; e ∈ W be a reaction such that Cons I (r v ) = {Y}.We consider the second condition, on the linearity of Y in e.If v; e ∈ W (that is the flux has not been changed at all by the simplification rule), then the linearity condition is trivially preserved by induction.v; e cannot be the simplification of a flux v; e with X as modifier, since that would contradict the linearity condition in W.  Then the set LinNets is stable for the simplification.This directly implies that the simplification can remove every intermediate species in I.

Confluence of the Simplification Relation
We now study the confluence of the simplification relation.We first show that the structural confluence, that is the confluence of the fluxes without kinetics, is a direct consequence of the previous results.We next present an example that illustrates that, however, the distribution of the kinetics between the fluxes can be different.Finally, we give a criterion on the modes of a network, that guarantees the full confluence, that is confluence of the structure and the rates.
In the following, we only consider networks in LinNets.

Structural Confluence
We say that two constrained flux networks W = V&C and W = V &C are structurally similar, denoted W ∼ = struc W , if they have the same structure, that is the same fluxes when neglecting the kinetic expressions: Proof.This is a direct consequence of the stability of the LinNets (Proposition 4) and of the confluence of the simplification without kinetics (Theorem 2).

Non-Confluence of the Kinetic Rates
Let us consider the reaction network W, depicted on Figure 22, with 7 species and 6 fluxes.The intermediate species are X, Y and Z.Initially, the kinetics are all mass-action.
We can remove the intermediate species in different orders with (C-INTER).If we start by eliminating X, followed by Y and finally Z, we obtain the network W XYZ , while if we first eliminate X, then Z and Y, we obtain W XZY .These networks are different, since W XYZ has one additional flux.This illustrates the necessity of the rule (C-DEP) to obtain the same network structure.
Indeed, the additional flux v 123456 in W XYZ is dependent on v 123 and v 456 and can therefore be removed with (C-DEP), while updating the kinetic expressions.We then obtain exactly the network W XZY .However, v 123456 also depends on v 25 and v 1346 .Therefore, we could as well remove it while updating these reactions.In that case, we obtain the different network W XYZd .This network has the same structure as W XZY , but not the same distribution of rates between the fluxes.

Criterion for the Full Confluence
We now give a criterion that guaranties the full confluence of the simplification.Definition 15.A vector of reactions r = (r 1 , . . ., r n ) is uniquely decomposable if any mode v ∈ ker + (S) has an unique decomposition in elementary modes, where S is the stoichiometric matrix of N r |I .
Then the mode v = (1, 1, 1, 1, 1, 1) can be decompose in either v 1 + v 4 , or in v 2 + v 3 .The two decompositions are illustrated in Figure 23.r is not uniquely decomposable, and the simplification is not confluent for the kinetic rates, as we have seen before.Theorem 4 (Confluence).If the initial vector of reactions r is uniquely decomposable, then the relation C on (LinNets, ∼ =) is confluent, for both the structure and the kinetic rates.
The theorem is the consequence of the following lemmas, that analyze the different critical pairs.That is, if a network W can be simplified in two different manners into W 1 and W 2 , then these two networks can be simplified into W 1 and W 2 such that W 1 ∼ = W 2 .
Lemma 7. Assume r is uniquely decomposable.Let W be a network such that W C-DEP W i for i ∈ {1, 2}.
Then ∃W i such that W i * Proof.Let v i e i be the dependent flux removed when simplifying W into W i , for i ∈ {1, 2}, with v i ; e i dependent on v 1 i ; e 1 i , . . ., v So the simplified networks are trivially the same, that is Assume v 1 = v 2 , and that for any i ∈ {1, 2}, for any j, v i = v j 3−i , that is v 1 does not depend on v 2 and reciprocally.Then we can still remove v 1 in W 2 , and v 2 in W 1 , and we find the same network modulo similarity: If v 1 depends on v 2 , and v 2 depends on v 1 , then since dependencies are positive linear combinations, that directly implies Proof.This case is trivial, since the substitutions commute.
Lemma 11.Let W be a network such that W Proof.Once again, since the two removed species cannot be the same, this case is trivial.
Proof.The full proof is given in Appendix C. The idea is that after removing one intermediate species, we can still remove the other one, either with (C-MOD) or with (C-INTER).In the second case, some dependent fluxes are generated, that we can eliminate to find the same simplified network, whatever the order of elimination of the intermediate species.

An Example from the BioModels Database
We have shown that the simplification system that we presented can exhibit non-confluence of the rates, even in a simple scenario with a small number of intermediates.To find if such a situation occurs in practice, we investigated the SBML models in the curated BioModels database [6].We were thus able to find a network, for the model BIOMD0000000173, that does not verify the confluence criterion, and such that two different simplified networks can be identified.Note that this was the only model not satisfying the criterion, when considering every model of the BioModels database with mass-action kinetics and with three or four linear intermediate species.
The network identified is a model of the Smad-based signal transduction mechanisms from the cell membrane to the nucleus, presented in [23].We only consider here a sub-network W of this model, sufficient to illustrate the non-confluence.It is represented in Figure 24.
In this network, a molecule of S4 c , that represents the species Smad4 in the cytoplasm, can bind with either a molecule of Smad2 in a phosphorylated form (pS2 c ) and form the complex S24 c (reaction r 5 ), or with a molecule of G in a phosphorylated form (pG c ), and form the complex G4 c (reaction r 22 ).These two reactions are reversible (r 5 and r 22 ).The same transformations can occur in the nucleus (reactions r 6 , r 6 , r 23 and r 23 ).The species Smad4 can also move from the cytoplasm to the nucleus, or reciprocally (r 1 and r 1 ).Finally, the complex of Smad2 and Smad4 can move from the cytoplasm to the nucleus (r 7 ).
Ck off G4c   We assume that I = {S4 c , S24 c , S4 n , S24 n }.The network is in LinNets.Therefore, we can consider the elimination of the four intermediate species.According to the order of the simplification, we can then obtain two different networks, with the same structure, but with different kinetic expressions.They are represented in Figure 25.The network W 1 is obtained by removing S4 n first, then S24 n , then S24 c , and finally S4 c and the dependent fluxes.The network W 2 is obtained by removing S4 n , then S24 n , S4 c , S24 c and the dependent fluxes.In We now show that the criterion is not satisfied in the initial network, i.e., that W is not uniquely decomposable.We consider the following mode: This mode has two possible decompositions into elementary modes, {v red , v blue } and {v green , v magenta }, with: We represent these fluxes in Figure 26, where we omit the non-intermediate species and the kinetic expressions for the sake of simplicity.

Simplification of Systems of Equations
In this section, we study the relation between the simplification C on reaction networks and a simplification ⇒ on systems.We show that the assignment of a system E(W) to a network W is a simulation for the simplifications.

Simplification of Systems of Equations
The simplification of systems is illustrated in Figure 27.The first rule replaces a constant variable x with its initial value x(0), and ẋ with 0. The second rule extends the simplification to similar systems.We define the simplification: Lemma 13.The simplification is correct for the equivalence, that is: Proof.It is trivial, since the substitutions commute.

Simulation
The assignment of a system E(W) to a network W is a simulation from (LinNets, C ) to (Systems, ⇒ * ).This is the equation for A in the system E(W ) for the simplified network.Therefore, E(W) ∼ = E(W ).

Conclusions
We have first shown that when neglecting the kinetic expressions, the elimination of linear intermediate species and dependent reactions is a reformulation of the double description method, that computes the elementary modes, and therefore that the network structure of simplified networks is unique.In a second time, when considering kinetic expressions, we provided a biological example illustrating that the simplification can produce two networks with the same structure but different kinetics.We then gave a sufficient criterion on the network structure of the initial network that guarantees the confluence of both the structure and the rates.
Note that the criterion seems to be satisfied in most cases in practice.When looking at the networks with mass-action kinetics from the BioModels database [6], and considering at most four intermediate species, only the Smad-model BIOMD0000000173 was identified as not satisfying the criterion.On the other hand, the linearly steadiness as well as the conditions required for a network to be in LinNets (such as the stoichiometry conditions, etc.) are not always satisfied in real biological networks, and these are therefore a real restriction on our simplification approach.
• we denote the vector of a path by ∑ ṽ = ∑ 1≤i≤k v i ; • a path is circular if Prod I (r v k ) = Cons I (r v 1 ) = ∅, and non-circular otherwise; • for a circular path ṽ, we denote the number of intermediate species occurring in the path by ṽ(X) = |{1 ≤ i ≤ k | Prod I (r v i ) = {X}}|; • for a non-circular path ṽ, we denote its beginning and its end by Cons I (r ṽ) = Cons I (r v 1 ) and Prod I (r ṽ) = Prod I (r v k ).In addition, we define the multiset ṽ(X) = |{1 ≤ i < k | Prod I (r v i ) = Cons I (r v i+1 ) = {X}}|.Note that we do not count Cons I (r ṽ) and Prod I (r ṽ) in this multiset.
Example A1.For instance, consider the initial network W 0 , with I = {X, Y, Z}, in Figure A1 (left).We denote by v i the unit vector for reaction r i .
The path v 3 v 4 is circular.It has for vector (0, 0, 1, 1, 0).The multiset is defined by v 3 v 4 (X) = 0 and We will see later that if a path ṽ satisfies some conditions w.r.t.some reaction network W, then there is a corresponding flux v in W such that ∑ ṽ = v.The reciprocal property also holds.Such a particular path, called a flux-path, is formally defined as follows.
Definition A2.Let W, W 0 be reaction networks such that W 0 * C W,

•
a non-circular flux-path ṽ is a non-circular path in W such that for any intermediate species X with ṽ(X) > 0, we have X ∈ Spec(W 0 )\Spec(W) (meaning that one of the simplification steps W 0 * C W removes X from W 0 ), and such that Cons I (r ṽ) and Prod I (r ṽ) are either the empty solution ∅, or the intermediate species that are still in W, • a circular flux-path ṽ is a circular path in W if there is at most one intermediate species X such that ṽ(X) > 0 and X ∈ Spec(W) (i.e., X is not yet simplified), • we call flux-path a path that is either a circular or a non-circular flux-path.• a flux-path ṽ is said to correspond to a flux v in W if ∑ ṽ = v.
Example A2.Let us consider the simplified network W in Figure A1 (right).The non-circular path v 1 v 2 v 3 is not a non-circular flux-path, since v 1 v 2 v 3 (X) > 0 and X has not been removed.The path v 2 v 3 is a non-circular flux-path, and corresponds to the flux v 23 .
The path v 3 v 4 is a circular flux-path, since there is a unique intermediate species, Z, that has not been removed and such that v 3 v 4 (Z) > 0. It corresponds to v 34 .
We first prove the following lemma on the dependent flux.Lemma A1.Let v; e be a flux that depends on v 1 ; e 1 , . . ., v k ; e k with coefficients a 1 , . . ., a k .For any i, let ṽi be a corresponding flux-path for v i .Then: v = ∑ 1≤i≤k a i ∑ ṽi .
Proof.We directly have v = ∑ 1≤i≤k a i v i and for any i, v i = ∑ ṽi .
We can now prove the key lemma of this section.Lemma A2.Let W, W 0 be reaction networks such that W 0 ∈ LinNets and W 0 * C W. Then the following properties hold: 1. for any flux v; e ∈ W, there is a corresponding flux-path ṽ for W such that ∑ ṽ = v.Moreover, if v is not dependent then, for any intermediate species X, we have ṽ(X) ≤ 1; 2. for any flux-path ṽ for W such that for any X, ṽ(X) ≤ 1, there is a corresponding flux v; e ∈ W, that is ∑ ṽ = v; 3. if v; e ∈ W depends on v 1 ; e 1 , . . ., v k ; e k ∈ W, then

•
there exists an index i such that Cons I (r v i ) = Cons I (r v ), Prod I (r v i ) = Prod I (r v ), and • for any j = i, Prod I (r v j ) = Cons I (r v j ) = ∅ and any flux-path ṽj that corresponds to v is circular.
Proof.We proceed by induction on the simplification steps.We start by proving each conclusion of the Lemma for the base case, that is W = W 0 .
(1) for flux v; e in the initial network W 0 , v is necessarily a unary vector v i for some i.So we can directly associate the flux-path ṽ = v i that trivially corresponds to v. Since it is a unary vector, the flux v is also necessarily not dependent.Because a flux-path of size 1 is always non-circular, we also have, for any intermediate species X, ṽ(X) = 0, and thus ṽ(X) ≤ 1 as required.
(2) any flux-path ṽ for W 0 is necessarily of size 1.Otherwise, for ṽ being a non-circular flux-path, there would exist some X ∈ Specs(W 0 )\Specs(W 0 ) = ∅.And for ṽ being a circular flux-path, there would exist at least two species X and Y such that ṽ(X) > 0 and ṽ(Y) > 0, which contradicts the definition of circular flux-path.Then there exists v i ; e ∈ W 0 such that v i = ṽ.(3) as said above, a flux v; e ∈ W 0 v can not be dependent.Now, considering the inductive case, we assume that the Lemma is true for a network W such that W 0 (1) Let v; e ∈ W, then it is the case that v; e ∈ W for some expression e because the rule (C-DEP) only removes a dependent flux and modifies some kinetic expressions.By induction hypothesis, there is a flux-path ṽ for W such that ∑ ṽ = v.Also, because Specs(W) = Specs(W ), any flux-path for W is also a flux-path for W, which proves that ṽ is a flux-path for v. Finally, if v is dependent in W, it is necessarily dependent in W and satisfies ∀X ∈ I. ṽ(X) ≤ 1 by induction hypothesis.(2) Let ṽ be a flux-path for W such that, for any intermediate species X, ṽ(X) ≤ 1.Again, because Specs(W) = Specs(W ), ṽ is also a flux-path for W .By induction hypothesis, there is a corresponding flux v; e ∈ W .If v; e is not the flux that is removed by the application of (C-DEP), then this flux still occurs in W (possibly with an updated kinetic) and we conclude directly.We now show that it can not actually be otherwise, and more precisely, that assuming v removed by (C-DEP) contradicts ṽ(X) ≤ 1.

Figure 1 . 4 −
Figure 1.Evolution of the concentration of S, E, C and P in enzymatic network with mass-action kinetics with the parameters k 1 = k 2 = k 3 = 1, the initial concentrations E(0) = 1, C(0) = 2, S(0) = 4, P(0) = 0, and in the context of the network with a reaction ∅ k 4 − → S which produces S with constant

Figure 5 .
Figure 5.The elementary modes of the reaction network in Figure 4.

Figure 6 .
Figure 6.Simplification of reaction networks without kinetics with respect to a set I of intermediate species.

Figure 7 .
Figure 7. Elimination of intermediates X and Y in reaction network N in both possible orders, leading to two different final results N XY and N YX .
. At the bottom, we show network N Y , obtained by eliminating the intermediate species Y first.The only reaction producing Y in N is r 2 and the only reaction consuming Y is r 3 .Merging them produces reaction r 23 .When eliminating intermediate X from N Y , we obtain network N YX on the bottom right.

Figure 9 .
Figure 9. Simplifying flux networks for an initial n-tuple of reactions r and a set of intermediate species I.

Figure 10 .
Figure 10.Elimination of intermediate species from flux networks in different orders is not confluent without factorization.

Corollary 2 .
The simplification relation F restricted to flux networks in the set {V | V n * F V} is confluent.

Figure 13 .
Figure 13.Systems of constrained equations with ODEs.

Figure 14 .
Figure 14.System of constrained equations of a constrained flux network V&C.

Figure 16 .
Figure 16.Simplification rules for n-ary constrained flux networks, with I the set of intermediate species and r the n-tuple of initial reactions.

Figure 17 .
Figure 17.Reaction networks for the Michaelis-Menten example.MMnet E and MMnet C are obtained from the initial network MMnet after removing E and C respectively.MMnet CE is obtained after removing both C and then E in this order.MMnet EC is obtained by inverting the order of elimination.

Figure 18 .
Figure 18.Example illustrating the need of Condition 2.

Figure 19 .
Figure 19.Example illustrating the need of Condition 3.

2 Figure 20 .
Figure 20.Example illustrating the need of Condition 4.

Figure 21 .
Figure 21.Example illustrating the need of Condition 4.
So now assume the v; e is the merging of a flux v p ; e p and v c ; e c ∈ W. Then we have Cons I (r v p ) = {Y} and Prod I (r v p ) = Cons I (r v c ) = {X}.Therefore, the linearity condition implies e p = Ye p and e c = Xe c , with for any Z ∈ I, Z / ∈ Spec(e p ), Spec(e c ).Then we have e = Ye p e c /cons, and the linearity condition is satisfied in W . Now let v; e ∈ W be a reaction such that Cons I (r v ) = ∅ and Prod I (r v ) = {Y}, and consider the third linearity condition.By linearity, v; e cannot be the simplification of a reaction of W with X as modifier.If we had v; e ∈ W, then the condition is satisfied in W by induction, and therefore in W too. Assume that v; e is the merging of a reaction v p ; e p and v c ; e c ∈ W. Then we have Cons I (r v p ) = ∅, Prod I (r v p ) = Cons I (r v c ) = {X}, and Prod I (r v c ) = {Y}.Therefore, the linearity conditions on W imply e c = Xe c , with for any Z ∈ I, Z / ∈ Spec(e p ), Spec(e c ).Then we have e = e p e c /cons, and the condition is satisfied in W .

Finally, consider the
stoichiometric condition.We only have to verify this property for new fluxes v; e that are the merging of a flux v p ; e p and v c ; e c .We have Prod I (r v p ) = Cons I (r v c ) = {X}.Moreover, by normalization, we haveProd I (r v p v c ) = Prod I (r v c )\Cons I (r v p ).Since |Prod I (r v c )| ≤ 1, we have |Prod I (r v p v c )| ≤ 1, and similarly for Cons I (r v p v c ).

Figure 22 .
Figure 22.Network W and its simplifications.(top left) Network W. (top right) Network W XYZ after eliminating X, Y and Z (in this order).(bottom left) Network W XZY after eliminating X, Z and Y. (bottom right) Network W XYZd after eliminating X, Y, Z, and the dependent reaction.The new parameter is K = k 2 k 3 + k 3 k 4 + k 4 k 5 .

Figure 25 .
Figure25.Simplified networks from W. Both networks have the same structure, and the kinetic expressions are defined in the table.The network W 1 is obtained by removing in order S4 n , S24 n , S24 c , S4 c and the dependent fluxes.The network W 2 is obtained by removing in order S4 n , S24 n , S4 c , S24 c and the dependent fluxes.

Figure 26 .
Figure26.In red the elementary mode v red , in blue v blue , in green v green , and in magenta v magenta .

Figure 27 .
Figure 27.Simplification rules for systems of equations.

Proof.
The rule (C-SIM) for the networks is directly imitated by the rule (E-SIM) for the systems.Since we have rv = ∑ 1≤i≤k a i r v i , the equation is similar to: Ȧ = ∑ 1≤i≤k a i r v i (A)e + ∑ 1≤i≤k r v i (A)e i + ∑ v ;e ∈V r v (A)e = ∑ 1≤i≤k r v i (A)(e i + a i e) + ∑ v ;e ∈V r v (A)e .

kC
W (for some k > 0) and W C W. If W C-MOD W or W C-SIM W, only the kinetics are modified between W and W. Therefore, the Lemma is still true in W. It remains to investigate the cases W C-DEP W and W C-INTER W. (C-DEP) Assuming that W C-DEP W, we prove that each point of the Lemma is satisfied by W.
and v 1 depends on v 2 , with coefficient a, but v 2 does not depend on v 1 (or conversely), we have v 1 = ∑ Let W be a network such that W C-DEP W 1 and W C-MOD W 2 .Then ∃W i such that W i * This case is quite trivial.Let v; e be the dependent flux, X the modifier, and v ; e another reaction such that v depends on v with factor a. If we remove X first and then v, the flux v ; e is simplified into v ; e [X := X(0)] + ae[X := X(0)].Otherwise, it is simplified into v ; (e + ae)[X := X(0)].The two expressions are trivially similar.Let W be a network such that W C-DEP W 1 and W C-INTER W 2 .Then ∃W i such that W i * Let W be a network such that W C-MOD W i for i ∈ {1, 2}.Then ∃W i such that W i * j 2 .If we remove v 1 , v 2 ; e 2 becomes v 2 ; e 2 + ae 1 , and can be removed.The fluxes obtained are of the form v j ; e j + a j 1 e 1 + a j 2 (e 2 + ae 1 ).If we remove v 2 first, then we can remark that v 1 is still dependent, with v 1 = ∑ Proof.Proof.The full proof is given in Appendix C. The main idea is to use Proposition 3 to prove that if X is in the dependent flux v d , it is also in one of the fluxes v i that v d depends on.Therefore, if we eliminate X and combine v d with another flux v , we also merge v i with v .Then v d X v is still dependent on v i X v and other fluxes, and can be removed.Lemma 10.