Stability of Metabolic Networks via Linear-In-Flux-Expressions

The methodology named LIFE (Linear-in-Flux-Expressions) was developed with the purpose of simulating and analyzing large metabolic systems. With LIFE, the number of model parameters is reduced by accounting for correlations among the parameters of the system. Perturbation analysis on LIFE systems results in less overall variability of the system, leading to results that more closely resemble empirical data. These systems can be associated to graphs, and characteristics of the graph give insight into the dynamics of the system. This work addresses two main problems: 1. for fixed metabolite levels, find all fluxes for which the metabolite levels are an equilibrium, and 2. for fixed fluxes, find all metabolite levels which are equilibria for the system. We characterize the set of solutions for both problems, and show general results relating stability of systems to the structure of the associated graph. We show that there is a structure of the graph necessary for stable dynamics. Along with these general results, we show how stability analysis from the fields of network flows, compartmental systems, control theory and Markov chains apply to LIFE systems.


1.
Introduction. Quantitative Systems Pharmacology (QSP) aims to gain more information about a potential drug treatment on a human patient before the more expensive stages of development begin [22]. QSP models allow us to perform in silico experiments on a simulated metabolic system that predicts the response of perturbing a flux. A drug may be metabolized differently by various patients, and modelers working in pharmacology must anticipate these differences. Building a profile of how the drug affects different classes of simulated patients will help the developers of new drugs understand the viability of a treatment and acquire insight into the mechanisms by which the drug acts.
A recent advancement in QSP modeling called Linear-in-Flux-Expressions (LIFE) is a method of analyzing systems of Ordinary Differential Equations (ODEs) [19,18]. Originally, LIFE was designed to analyze metabolic systems, which are composed of Fluxes and Metabolites. Fluxes in the metabolic system are the rates of chemical reactions in the human body, and they determine the dynamics on the metabolites, which are the various chemical compounds involved in metabolism. Modeling these systems depends on choosing fluxes, which are difficult to measure directly, so that the system effectively simulates human metabolism.
To implement LIFE on a metabolic network, the network must be written as a directed graph [18]. The edges of the graph represent the reaction rates (fluxes), and the vertices represent quantities of chemical compounds (metabolites). From the graph we construct the stoichiometric matrix of the system. This stoichiometric matrix is not the classical one mentioned by [21,24]. The LIFE method is also different from QSP models whose dynamics traditionally depend on a matrix containing information about the flux of the system. In these classical QSP models the dynamics of the metabolites are linear with respect to the metabolites. By contrast, systems using the LIFE method are linear in fluxes and have a stoichiometric matrix that is dependent on the metabolites.
Initially, the LIFE method was developed using the human cholesterol metabolism network [19]. LIFE enables us to simply describe the correlations among the fluxes of the model at steady state. There are generally many correlations among fluxes, and maintaining these correlations leads to a more consistent response to perturbing the fluxes in the system. This was advantageous to QSP modelers, who previously analyzed flux perturbations with little to no consideration to relationships among fluxes [1]. Now, we expand our study of these systems, showing that with few assumptions, systems that are linear in the flux is stable.
The LIFE method evolved from methods in systems biology [21]. Systems biology, in conjunction with network flows [12,9], Markov chains [6], laplacian dynamics [13], control theory [3], and compartmental systems [4,14] allow us to better understand biological networks on which pharmacology models are based. The field of compartmental systems focuses on models based on directed graphs. Vertices of the graph represent quantities whose dynamics are determined by the edges of the graph, which represent fluxes among compartments. Markov chains study dynamics on directed graphs as well, but by contrast, this field focuses on stochastic processes. Control theory studies the way an external agent can alter the natural evolution of a system, given a set of admissible controls. In pharmacology, metabolism follows its natural evolution, and drugs serve as our controls. These fields have much to contribute to systems pharmacology, and we summarize useful results. We identify assumptions which are usually satisfied by real metabolic networks, guaranteeing stability of the metabolic system at a unique equilibrium.
The paper is organized as follows. In section 2, we describe the model system for the LIFE approach in the form of a system of ODEs associated to a metabolic networks, then show existence of positive solutions and provide results of equilibria under general assumptions. Also special classes of LIFE systems are introduced. Section 3 investigates the flow vectors for which a given metabolite vector x is an equilibrium of the network, including the extreme pathways approach. On the other side, Section 4 studies the opposite problem: find the metabolite vectors which are equilibria of the network for a fixed flow vector. This is done first investigating the relationships between linear LIFE systems and Markov chains, Laplacian dynamics and linear compartmental systems. Then we deal with special classes of nonlinear LIFE systems. Finally, a comparison between zero-deficiency theory is discussed. The paper ends with conclusions in Section 5 and an Appendix containing examples. 2. System model.

Notation and preliminaries.
We indicate by R + = [0, +∞) the set of positive real numbers, by R n the Euclidean real space of dimension n and by M n×m the set of n × m matrices with real entries. Given a matrix S, we indicate by S T its transpose. Given d 1 , . . . , d n ∈ R, diag(d 1 , . . . , d n ) is the diagonal matrix with entries d i on the diagonal. We denote by 1 a column vector with all entries equal to 1, of size clear from the context.
We introduce some terminology commonly used in graph theory. A directed graph is a couple G = (V, E), with V = {v 1 , . . . , v n } the set of vertices and E ⊂ V × V the set of edges. For a graph with n vertices and m edges, ordering the edges lexicographically, the incidence matrix is a matrix, Γ ∈ M n×m such that Γ ij = 1 if the jth edge is (v k , v i ) for some vertex v k , Γ ij = −1 if the jth edge is (v i , v k ) for some vertex v k , and Γ ij = 0 otherwise. A path is a sequence of distinct vertices v i1 · · · v i k , with (v ij , v ij+1 ) ∈ E for j = 1, . . . , k − 1. A graph is strongly connected if there exists a path between every pair of vertices. A strongly connected component of a directed graph is a maximal strongly connected subgraph.
A terminal component of a directed graph G = (V, E) is a strongly connected component G ′ = (V ′ , E ′ ), with V ′ ⊂ V , E ′ ⊂ E, such that there exists no edge e = (v ′ , v), with v ′ ∈ V ′ and v ∈ V \ V ′ . An undirected path is a sequence of distinct vertices v i1 , · · · , v i k , with either (v ij , v ij+1 ) ∈ E or (v ij+1 , v ij ) ∈ E for j = 1, . . . , k − 1. A directed graph is weakly connected if there exists an undirected path between every pair of vertices. A weakly connected component of a directed graph is a maximal weakly connected subgraph. A directed graph G = (V, E) is weakly reversible if every weakly connected component is also strongly connected.
2.2. LIFE model. We indicate by x ∈ R n the metabolite variables and by f ∈ R m the flux variables. A general system of ODEs which governs the quantities of x and f is written as where F : R n ×R m → R n and G : R n ×R m → R m . In [11,15], the authors show that the dynamics described by (1) evolve over a much smaller time-scale than (2). This is referred to as "time-scale separation". Based on time-scale separation arguments of metabolic systems, we approximate the dynamics of the fluxes with G ≈ 0, and our work focuses on the dynamics of the metabolites (1), with the fluxes playing the role of constant parameters. The dynamics (1) is very general and we restrict to special system, which are linear in the fluxes, thus can be written as: where S : R n → M n×m is called the stoichiometric matrix. One constructs the stoichiometric matrix from the metabolites and the reactions that comprise a biochemical system. Each reaction corresponds to a flux f that connects two distinct metabolites or represents an intake or an excretion from the network. Each row of S corresponds to a metabolite and each column of S corresponds to a flux. We now illustrate how to construct a directed graph from the metabolic network for the system (3). We represent metabolites with vertices V = {v 1 , . . . , v n }. We construct a set of edges E ⊂ V × V to represent reactions; each edge is associated to a flux from one metabolite to another, notice that we do not have loops. To represent intakes and excretions, we introduce two virtual vertices, v 0 and v n+1 , not associated with any metabolite but rather representing the external environment. We denote by I, X the set of vertices attached to v 0 , v n+1 , the vertices in I and X are called intake vertices and excretion vertices, respectively. We also introduce edges (v 0 , w) with w ∈ I ⊂ V representing intakes, and (w, v n+1 ), w ∈ X ⊂ V representing excretions. We use the extended graphG = (Ṽ ,Ẽ) defined byṼ = V ∪ {v 0 , v n+1 } = {v 0 , v 1 , . . . , v n , v n+1 } andẼ collecting edges in E together with intake and excretion edges. The rows of the matrix S can be indexed by vertices in V and the columns by edges inẼ, thus we write S ve for the entry corresponding to vertex v and edge e. Moreover we denote by x v the metabolite corresponding to vertex v and by f e the flux corresponding to edge e. Example 2.1. To illustrate the concepts of graph with virtual vertices and stoichiometric matrix related to a metabolic network, we provide a toy example with linear dynamics. Consider the system given by the following stoichiometric matrix and fluxes vector: Then the corresponding graph is represented in Figure 2.
Remark 1. It is worth a reminder that our stoichiometric matrix is different from the traditional one defined by [21,24], in which entries are stoichiometric coefficients, i.e. do not depend on metabolites.
To correctly represent the reactions corresponding to fluxes (which take always strictly positive values), we assume: (A) For x ∈ (R + ) n , it holds where H e : R n → R is a positive continuous function.
namely all columns of S have zero sum, except those corresponding to intakes and excretions, which have positive and negative sum, respectively. Under Assumption (A), the dynamics (3) can be interpreted as mass conservation law. Indeed, rewriting (3) entrywise and using (A), we havė which is the mass balance for metabolite x v : its variation is given by the sum of the incoming flows, minus the sum of the outgoing flows. This is the analogous Kirchhoff's current law for electrical networks, with the difference that currents are allowed to take negative values as well, while here metabolite variables are nonnegative. The total mass in the system is m = v∈V x v . From (5) we havė Clearly, in the case without intakes nor excretions,ṁ = 0, i.e. the total mass of a closed system is constant in time.
Another remark which is useful is that, under Assumption (A), S(x) can be rewritten as S(x) = ΓD(x), where D(x) is a diagonal matrix of size m × m, with diagonal entries given by H e (x)'s, and Γ is obtained from the incidence matrix ofG by removing the first and last rows (corrseponding to v 0 and v n+1 ). In the particular case without intakes nor excretions, Γ is the incidence matrix of G.
In the remainder of this section we study the dynamics (3) under the very general Assumption (A), while later in the paper we add other assumptions, restricting our attention to systems for which stronger statements can be obtained. A first important general property is that positivity of solution is guaranteed:

Proposition 1. Consider a system (3) satisfying (A) and the Cauchy problem with
Proof. Existence follows from Lipschitz condition, while positivity of solution follows from the invariance of the set {x : x v ≥ 0}.
In the next Proposition we show that existence of nontrivial equilibria implies some structure on the network: every vertex v for which there is a directed path from some w ∈ I to v, must also have a directed path from v to some y ∈ X. This result refines the space of networks with which we are concerned. More precisely: Proposition 2. Consider a system (3) satisfying (A). Assume there exists an equilibriumx ∈ (R + ) n for a flux vector f such that f e > 0 for every e ∈Ẽ. Then for every vertex v ∈ V for which there exists a path from I to v, there exists a path from v to X.
Proof. Assume there exists an equilibriumx as in the statement and, by contradiction, a vertex v for which there exists a path from w ∈ I to v, but there exists no y ∈ X to which v is connected. Since there is no path from v to X, either v belongs to a terminal component with no excretion, or there is a path from v to a terminal component with no excretion. Denote by G T = (V T , E T ) such a terminal component. Since there are a path from w ∈ I to v and a (possibly trivial) path from v to V T , then there is also a path from v 0 to V T . Denote by v 0 , v 1 = w, . . . , v ℓ−1 , v ℓ ∈ V T one such a path, such that v ℓ−1 / ∈ V T (possibly the path is a single edge, in case w ∈ V T ).
Finally we have a terminal component G T = (V T , E T ) with no excretion, and an The system associated to the graph in Figure 2 provides an explicit example of the contradiction argument of Proposition 2. The vertices v 3 and v 4 violate the condition of Proposition 2, since there is a path from the intake vertex v 1 to v 3 and v 4 , and there is no path from v 3 and v 4 to the excretion vertex v 2 ; hence, the system does not admit any equilibrium. The proof argument, specialized to this example, is to notice that v 3 and v 4 form a terminal component with no excretion, and to look at the path v 0 , v 1 , v 2 , v 4 . Assuming that there is an equilibriumx, one shows first thatx v1 > 0 due to the intake, and then thatx v1 > 0 implies alsox v2 > 0. Notice thatẋ v3 +ẋ v4 ≥ S v2,(v2,v4) f v2,v4 (in this particular example, equality is actually true), and the latter is > 0 sincex v2 > 0, thus contradicting the fact thatx is an equilibrium: the mass x v3 + x v4 grows unbounded.
The same system of Figure 2 shows the necessity of the assumption that H e (x) = 0 when x v = 0, e = (v, w) for Proposition 2 to hold. Indeed assume that H (v2,v5) (x) = 1 for every x such that x v2 = 0 and that the other functions H e , e = (v, w), satisfy H e (x) = x v . Let f be a flux vector with strictly positive components such Another example not satisfying the assumptions of Proposition 2 is shown in Figure 3. This system admits an equilibrium under Assumption (A) even though the vertex v 3 does not have a path to v 5 . This is because there does not exist a path  Figure 4. Under Assumption (A), the dynamics is written as: . . .
. We want to show that for any dynamical system in R n on a compact set there exists an equivalent dynamics defined on the cycle graph of Figure 4. In other words, Asumption (A) is so general that we give Proposition 3 to show arbitrary dynamics can be defined, and we focus on more specialized cases (Assumption (B), and (C)). More precisely for every general dynamics . . .
we look for a choice of the functions H e and the fluxes f e realizing such equivalence. Note that from (4), we have x vn = C − x 1 − · · · − x n−1 which implies thaṫ Define the set T = ∪ n i=1 {x : x j ≥ 0 for all j, x i = 0 and x 1 + · · · + x n ≤ 1} then we have the following: Proposition 3. Consider system (7) and assume F 1 = · · · = F n = 0 on the set T . Then there exist functions H e and fluxes f e such that the dynamics (6) is equivalent to (7) on the bounded set delimited by T .
Proof. Assign the functions H (vi,vj ) (x) according to the following rule: where [F ] + = max{F, 0} and [F ] − = max{−F, 0}. It is easy to verify thatẋ vi = F i (x) for i < n andẋ vn = − n i=1 F i = F n (x) using (8). Moreover, because of the assumption on F i 's, the functions H e satisfy (A), thus we are done.

Special LIFE systems.
Here we introduce a special class of systems of type (3) with simplified dynamics. We consider Assumption (A), and we impose further restriction on the functions H e (x): for an edge e = (v, w), we assume H e to depend on x v only, and moreover we impose the scalar function H e (x v ) to be strictly increasing. More precisely, Assumption (B) is the following. (B) It holds and each H e is a positive, differentiable, and strictly increasing function H e : R → R, with H e (0) = 0. A typical example of a system verifying (B) is given by metabolic networks with Hill functions representing reactions, i.e.
A further simplification occurs in the case where H e (x v ) is the same function H v (x v ) for all edges e having v as an initial vertex, i.e. such that e = (v, w) for some w. This gives Assumption (C), as follows.
and each H v is a positive, differentiable, and strictly increasing function H v : R → R, with H v (0) = 0. Under Assumption (C), the systemẋ = S(x)f can be equivalently re-written aṡ where J(f ) ∈ M n×n is defined by Finally, the simplest class of LIFE models we consider are linear systems, namely systems satisfying Example 2.1 (continued). This example is a linear LIFE system. We can equivalently re-write its dynamicsẋ = S(x)f asẋ = J(f )x + φ, with 3. Equilibria for fixed metabolites. In this section we consider equilibrium solutions of the system (3) satisfying Assumption (A). In general one is interested in conditions guaranteeing existence of an equilibrium and also in conditions necessary for uniqueness and stability of such an equilibrium. There are two problems: for fixed metabolite concentrations x find all flux vectors f for which x is an equilibrium, i.e.ẋ = S(x) · f = 0, and, vice versa, for a fixed flux vector f , find all x that are equilibria. In this section we focus on the first, while the latter is investigated in Section 4. The set of flux vectors for which x is an equilibrium formed by all vectors f that solve the equation S(x) · f = 0, i.e. the nullspace N (S(x)) of S(x). First, we discuss the dimension of N (S(x)). Then, since fluxes must be positive (to have a correct biological meaning), we focus on describing the cone 's as entries and Γ is obtained from the incidence matrix ofG by removing the first and last rows; in case there are no intakes nor excretions, then Γ is the incidence matrix of G. Assuming that x has strictly positive entries, (A) implies that all diagonal elements of D(x) are strictly positive, so that D(x) is invertible. Hence, for Γ · (D(x) · f ) = 0 to have non trivial solutions, the nullspace of Γ must have dimension greater than zero.
We extend one of the results of [18] to get:

Proposition 4. Consider the system (3) with no intakes nor excretions satisfying Assumption (A) and let
G be the the associated graph. Let x ∈ R n have strictly positive entries and S(x) ∈ M n×m be the stoichiometric matrix, then where ℓ is the number of weakly connected components of G.
Proof. We have S(x) = ΓD(x). As shown in Proposition 4.3 of [2] the rank of an incidence matrix Γ is rank(Γ) = n − ℓ. Since D(x) is full rank we have rank(S(x)) = n − ℓ.
Next we consider systems that contain both intakes and excretions. It was shown in [18] that if intakes and excretions are added to a graph satisfying the conditions of Proposition 4 then rank(S(x)) = n. We extend this result to get: Proposition 5. Consider the system (3) satisfying Assumption (A) and let G be the the associated graph. Let x ∈ R n have strictly positive entries and S(x) ∈ M n×m be the stoichiometric matrix, then rank(S(x)) = n − k, where k is the number of weakly connected components containing neither intake nor excretion vertices.
Proof. It was shown in [18] that if G is weakly connected and contains an intake vertex then rank(S(x)) = n. The same argument from [18] can be also used for an excretion vertex and so if G is weakly connected and contains an excretion vertex we also have rank(S(x)) = n. In the case where G is not weakly connected S(x) can be rewritten as a block diagonal matrix where each diagonal element S i (x) is the stoichiometric matrix of a weakly connected component. If S i (x) contains an intake or excretion vertex then rank(S i (x)) = n i , where n i is the number of metabolites in S i (x). If S i (x) contains neither intake nor excretion vertices then rank(S i (x)) = n i − 1. This implies that S(x) has n i linearly independent rows for each S i (x) with intake or excretion vertices and n i − 1 linearly independent rows for each S i (x) with neither intake nor excretion vertices. Let k represents the number of S i (x) with neither intake nor excretion vertices, the total number of linearly independent rows in S is now expressed ( Notice that the rank of S(x) depends on the number of weakly connected components of the graph, which is the same irrespective of the orientation of edges. However, when we focus on the existence of non-trivial positive flows admitting an equilibrium with strictly positive entries, the orientation of edges does matter, as we can see e.g. from Proposition 2, where a necessary condition is given for existence of equilibria, in terms of existence of suitable paths.
3.1. Network flows. The problem of finding positive flows admitting an equilibrium with strictly positive entries has been extensively studied in the operations research literature under the name of 'network flow' problems [12]. In network flow problems one considers a directed graph where edges represent flows between the vertices. The maximum flow which an edge can support is called the capacity of the edge. In addition to capacity each edge may also have a cost associated to it. Network flow problems assume that the system is at equilibrium with respect to the vertices i.e. that the flow entering and leaving a vertex must be the same. A flow that satisfies this assumption is called a feasible flow. Finding feasible flows is the same as finding flows in the cone N (S(x)) ∩ (R + ) m with the further constraint that each flux must be less than its capacity.
A flow is a mapping from f : [17]. One of the most common network flow problems is to find the maximum flow of the network, i.e. the largest amount of total flow from a source to a sink. In [9] the authors consider a network which contains exactly one source (v 0 ) and one sink (v n+1 ). They then prove the following result, known as max-flow min-cut theorem, which characterizes the maximum flow as the minimum cut cacapity, where a cut set is a set of edges whose removal disconnects the source from the sink, and the capacity of a cut is the sum of the capacities of its edges.
If for all v ∈ I there is a path from v to X (and hence a path from v to v n+1 ), the minimum cut is the removal of all edges (v 0 , v). The maximum flow ϕ * then has ϕ * v0,v =f v0,v , thus also ensuring the existence of an equilibrium flow f * satisfying the same.
If for some v i ∈ I there is no path from v to X, then all feasible flows ϕ satisfy ϕ v0,v = 0, and hence all equilibrium flows f satisfy f v0,v = 0 which contradicts the assumption.
3.2. Extreme pathway algorithm for calculating the positive basis. There are many standard methods for computing a basis of the nullspace of a matrix [20], and hence to describe N (S(x)). Since we are interested in positive fluxes, we are rather interested in the cone N (S(x)) ∩ (R + ) m . One method to describe this cone is to look for a positive basis, which is a minimal set of vectors generating the whole cone via linear combinations with positive scalars, called positive combinations. The use of positive combinations ensures that all generated vectors belong to the cone. Minimality is equivalent to ask for the vectors to be positively linearly independent, i.e. no vector can be expressed as a positive combination of the others. Since we allow only positive combinations (and not every linear combination), typically a positive basis has more vectors than a traditional basis. In [23] the authors prove that a flux cone has a positive basis, and that the vectors of the positive basis are unique up to multiplication by a positive scalar.
Finding a positive basis for a cone is not as easy as finding a traditional basis. In [21], Palsson describes a method to find a positive basis by using extreme pathways of a metabolic network. The extreme pathways are a unique set of positively independent flux vectors that represent the edges of the cone N (S(x)) ∩ (R + ) m [21,23].
Here we summarize a method to build extreme pathways described by [23]. We start considering the equation V 0 · S(x) T = C 0 , and the solution given by V 0 = I m , the m × m identity matrix, S(x) the stoichiometric matrix and Algorithm 1. The algorithm consists of the following steps: Step 1. Select from C i the first non-zero column (lexicographically), whose corresponding metabolite is neither an intake vertex nor an excretion vertex, say j i . If no such column exists jump to Step 5.
Step 2. Add to C i all possible new rows, obtained using positive linear combinations of two other rows, which have a zero on the column j i . Then define C i+1 by removing the old rows used to generate the new ones. Notice that C i+1 has all zeros in column j i .
Step 3. Define V i+1 by adding to V i the rows generated by the same combination as those of Step 2 (to keep the validity of the equation V · S(x) T = C) and removing the old rows.
Step 4. Remove positively linearly dependent rows from V i+1 . Remove the corresponding rows from C i+1 .
Step 5. Repeat the steps 1-4 now considering first excretion vertices and then intake vertices until all columns are 0.
If the algorithm stops at step m, then C m = 0 (the matrix with all zero entries). The rows of V m are the vectors of a positive basis. Notice also that the number of rows n i may increase or decrease during the various steps.
Originally in [21] the extreme pathways method was proposed for stoichiometric matrices not dependent on x. We apply it here to the most general case of systems (3) with Assumption (A), by fixing a desired metabolite variables vector x ∈ (R + ) n and using the extreme pathways method to characterize all fluxes vectors with positive entries such that the corresponding LIFE system (3) admits x as an equilibrium.
To illustrate the process we report an example using the Reverse Cholesterol Transport Network (RCT) from [19]. The RCT network is shown in Figure 5.
metabolic network Figure 5. Reverse Cholesterol Transport Network from [19]. This network contains 6 vertices which represent metabolites, 10 edges which represent fluxes and 2 virtual vertices v 0 , v n+1 . There are three intake vertices v 1 , v 2 , v 3 and 1 excretion vertex v 6 .
Example 3.1. Finding Extreme Pathways of RCT network We use Palsson's algorithm for finding extreme pathways on the stoichiometric matrix from the RCT network. The stoichiometric matrix for RCT is given by: For space, the complete calculations are continued in the Appendix.
3.3. Extreme pathways provide a positive basis. Let us first recall Farkas' Lemma [10]. Here we use the notation x ≥ 0 to indicate that every entry in the vector x is positive. 1. There exists x ∈ R m such that Ax = b and x ≥ 0.
2. There exists y ∈ R n such that A T y ≥ 0 and b T y < 0.
Proof. There are three statements which must be shown to prove this Proposition.

Statement 1:
The extreme pathway vectors (v i ) are positively independent (systematically independent in terminology used in [21]). Statement 2: The span of the extreme pathways vectors is contained in Finally it remains to prove statement 3, in order to show (N (S(x)) ∩ (R + ) m ) ⊂ λi≥0 λ i v i we proceed as follows. First we show that after one iteration of the algorithm the rows of V 1 form a positive basis for vectors w such that the j 1 -th element of S(x)w equals zero. After n iterations, we show that V n is a positive basis for vectors w ∈ N (S(x)).
For simplicity of notation, we assume j 1 = 1. Define W 1 to be the set of vectors w such that the first entry of S(x)w equals zero. It is clear from the construction of V 1 that its rows belong to W 1 . It remains to be shown that the rows of V 1 span W 1 . Since V 1 · S(x) T = C 1 , thus we can alternatively show that the columns of C T 1 span the subset of the range of S(x) with first entry equal to zero. In other words, we must show that if λi≥0 λ i c i,1 = 0, where c i,1 represents the first entry in the c i -th row of C 0 , that the sum λi≥0 λ i c i = 0 can also be represented using rows from C 1 . We assume an ordering such that c i,1 ≥ 0 for i ≤ ℓ and c i,1 < 0 for i > ℓ, each row in C 1 can be represented by c i + α i,j c j where i ≤ ℓ, j > ℓ and α i,j = −ci,1 cj,1 ≥ 0. We need to find µ i,j such that, We can split the sum on the right hand side by considering the positive and negative entries separately, For λ i , i ∈ {ℓ + 1, . . . , n} we can write: For λ i , i ∈ {1, . . . , ℓ} we can write: The equations (13), (14) can be written in the following way, where [c i,1 ] represents a 1×(n−ℓ) vector with all entries equal to c i,1 and diag(c i,1 ) ∈ M (n−ℓ)×(n−ℓ) a diagonal matrix with c i,1 on the diagonal. Let A, b be the matrix on the left-hand side and the vector on the right-hand of equation (15), respectively. We can then apply Farkas' Lemma. Consider a vector y ∈ R n such that A T · y ≥ 0, then for the first row we have c 1,1 y 1 − c 1,1 y ℓ+1 ≥ 0. Further for any i = 1, . . . , ℓ and j = ℓ + 1, . . . , n we have c i,1 y i − c i,1 y j ≥ 0 which implies y i ≥ y j . Settinḡ y = max(y ℓ+1 , . . . , y n ) we have y i ≥ȳ for i = 1, . . . , ℓ. Since all c j,1 < 0 for j = ℓ + 1, . . . , n and all λ j ≥ 0 we have that y j c j,1 λ j ≥ȳc j,1 λ j . Next we definē y = (ȳ, . . . ,ȳ) ∈ R n and note that b T y ≥ b Tȳ . Since λi≥0 λ i c i,1 = 0 we have b Tȳ = 0. Therefore condition 2 of Farkas' Lemma fails. This implies that condition 1 is true and thus there exists a set of positive µ's that solve the system (15). Since a set of µ's can be found we have that after the first iteration of the algorithm the rows of V 1 span W 1 . By similar argument the rows of the V i in subsequent iterations span the set of vectors w such that the first i entries of S(x)w equal zero. After n iterations, the rows of V n span N (S(x)) ∩ (R + ) m , thus we are done. 4. Equilibria and asymptotic behavior of metabolites for fixed fluxes. In this section we characterize equilibria and the stability of LIFE systems with a fixed flux vector. We start with linear systems, and then we consider special LIFE systems, with Assumptions (B) and (C). LIFE systems satisfying at least Assumption (A) are known as compartmental systems in the automatic control community [14,4]. We build upon the well-established results on linear compartmental systems to get a full understanding of linear LIFE systems, as well as for special LIFE systems satisfying Assumption (C). Notice that weakly connected components of G correspond to subsystems having no interaction with each other, so that we can study each weakly connected component separately from the others. Hence, throughout this section, we assumeG is weakly connected, without loss of generality.
4.1. Linear systems. In this section we study the properties of linear systems, i.e. special LIFE system (3) satisfying Assumption (C) with fluxes H v (x v ) = x v , or, equivalently, a system (11), with h(x) = x. We first focus on metabolic networks with no intakes nor excretions, recalling results from a rich literature from different communities. For the case with intakes and excretions, we make use of results for Compartmental Systems [4].

Linear system without intakes and excretions.
In the case with no intakes nor excretions (also known as free closed system), the linear LIFE dynamics becomė x = J(f )x. Notice that J(f ) is a Metzler matrix (i.e. has non-negative off-diagonal entries), and all its columns sum to zero (1 T J(f ) = 0 T ).
This system is known as Laplacian dynamics, because L = −J(f ) T is a weighted Laplacian of the graph G defined as follows. Given the graph G and given strictly positive weights f e associated with its edges (for us, the weights are the fluxes), the weighted adjacency matrix A is defined by The eigenvalues of the Laplacian, and in particular its eigenvalue 0 and the corresponding eigenspace which is actually the set of equilibria ofẋ = J(f )x, have been studied in graph theory (see e.g. [5]). The study of the spectrum of the Laplacian has received extensive attention also in the automatic control community, on one hand because of the recent interest in the so-called consensus dynamics (see [4,Chapter 7])ẋ = −Lx (which are different from the dynamicsẋ = J(f )x, since L = −J(f ) T , but the eigenvalues of −L and J(f ) are the same), and on the other hand because of the interest in the linear compartmental systemẋ = J(f )x (see [4,Chapter 9]).
Laplacian dynamics have been introduced in the mathematical biology literature by [13], together with a complete study of their equilibria and convergence properties. Notice that in [13] the matrix J(f ) itself is called a Laplacian, while in the graph theory and control theory communities the name Laplacian refers to The dynamicsẋ = J(f )x is also very related to continuous-time Markov chains (see e.g. the textbook [6]). Thus we elaborate on this, a homogeneous continuoustime Markov chain over a finite state V = {1, . . . , n} is a stochastic process X(t) taking values in V , that satisfies the Markov property: for all times t 0 ≤ t 1 ≤ · · · ≤ t h ≤ t h+1 , and which is homogeneous in time: P (t) is the matrix whose (i, j)-th entry represents the probability of transition from state i to state j in an interval of time of length t. P (t) is assumed to be rightdifferentiable, and with time derivative defined byṖ (t) = lim h→0 + (P (t + h) − P (t))/h. The Markov chain is fully described by its generator matrix (or transition rate matrix ) defined by Q =Ṗ (0). Indeed, P (t) is related to Q by the Kolmogorov forward equationṖ (t) = P (t)Q. The unique solution of the Cauchy problemṖ (t) = P (t)Q with P (0) = I is P (t) = e Qt . Denoting by π(t) a column vector whose ith entry is π i (t) = Pr(X(t) = i), we have that π T (t) = π T (0)P (t) = π T (0)e Qt and hence is a solution ofπ T = π T Q with given initial condition π(0).
The generator matrix Q is a Metzler matrix whose rows sum to 0 (Q1 = 0), and hence there is an immediate equivalence between the dynamicsπ T = π T Q and the linear LIFE systemẋ = J(f )x, simply by taking Q = J(f ) T . In Markov chains, one looks at π being a probability vector, namely having positive entries and i π i = 1. In LIFE systems, we are also interested in a vector x with positive entries, but the total mass could be arbitrary, so that one should set π i (0) = x vi (0) It is customary to 'normalize' the vectors π 1 , . . . , π k so that their entries sum to 1, so that they can be interpreted as probability vectors. Under this choice, the restriction of vector π i to the terminal component G i is the stationary distribution of the Markov chain restricted to such component, i.e. its unique equilibrium with mass 1 (uniqueness is obtained from the fact that the terminal component is a strongly connected component).
Notice that the dimension of the nullspace of an incidence matrix Γ (and hence of a stoichiometrix matrix S(x)) is related to the number of weakly connected components, while the dimension of the nullspace of a Laplacian matrix is related to strongly connected components, and not all components matter, but only the terminal ones.
In the case where G has a unique terminal component, the nullspace of J(f ) has dimension one: the equilibrium is unique, up to a multiplicative factor which is the initial total mass m 0 = i x vi (0). In the case where G is strongly connected, the whole graph is a single strongly connected component. In this case, not only the nullspace of J(f ) has dimension one, but also it is generated by a vector whose entries are all strictly positive.
The spectral properties of J(f ) given in Proposition 9 fully characterize the asymptotic behavior of the linear systemẋ = J(f )x, by standard theory of linear systems, giving the following: Proposition 10. Consider the linear LIFE system with no intakes nor excretions and denote by G 1 , . . . , G k the terminal components of the associated graph G. The following properties hold: • The total mass of the system m = v∈V x v = m 0 is constant in time.
• From any positive initial condition x(0), the system converges to an equilibrium, having strictly positive entries in correspondence of vertices of terminal components, and 0 elsewhere. • Moreover, if there is a unique terminal component, then the equilibrium is uniquely determined by the initial total mass m 0 .

Linear system with intakes and excretions.
We now focus on linear systems with intakes and/or excretions, using results from linear compartmental systems as summarized in [4,Chapter 9]. A first remark is that, having introduced a single v 0 from which all intake edges are originated, and a single v n+1 to which all excretion edges are headed, we have a single weakly connected component containing intakes and/or excretions. We restrict our attention to such a weakly connected component, while other possible components with no intakes nor excretions have a behavior described in Sect. 4.1.1. The dynamics are given byẋ = J(f )x + φ, where vector φ represents the intakes. Due to excretions J(f ) does not have all column-sums equal to 0. This means that −J(f ) T is not any more a Laplacian, but it is a grounded Laplacian. The term grounded Laplacian refers to a matrix obtained from a larger Laplacian matrix by deleting the row and column corresponding to a given vertex; the name 'grounded' has an interpretation for electrical networks, where this corresponds to connecting the given vertex to the ground. Consider the subgraph ofG where we have removed v 0 but not v n+1 , and consider the weighted Laplacian L ∈ M (n+1)×(n+1) of such a graph (with weights equal to the fluxes), then define L g by deleting the last row and last column (associated with v n+1 ). The resulting grounded Laplacian is L g is such that J(f ) = −L T g . The spectral properties of J(f ) are summarized as follows (see [4, Theorem 9.5 and Lemma 9.12]):

Proposition 11. Consider a linear system with intakes and/or excretions, then:
• All eigenvalues of J(f ) are either 0 or have strictly negative real part.
• The dimension of the nullspace of J(f ) is equal to the algebraic multiplicity of the 0 eigenvalue, and is equal to the number of terminal components not containing any excretion. In particular, the following are equivalent: (a) For every v ∈ V there is a path from v to X. To study the general case, where vertices might or might not have a path to some excretions, it is convenient to partition the system into two subsystems, as follows. Partition the vertex set as V = V 1 ∪ V 2 , with V 1 the set of v ∈ V such that there is a path from v to v n+1 . Without loss of generality, we can re-label vertices in V so as to have vertices v 1 , . . . , v r ∈ V 1 and v r+1 , . . . , v n ∈ V 2 . According to this decomposition, partition the vector x into two blocks x 1 corresponding to V 1 and x 2 corresponding to V 2 , and similarly partition φ as φ 1 , φ 2 . Notice that there is no edge from V 2 to V 1 , since an edge (w, v) with v ∈ V 1 implies that there is a path from v to v n+1 and also a path from w to v n+1 . Hence, we have: The first subsystem is called the reduced system, and its evolution is not affected by the second subsystem. The matrix J 1 is equal to the matrix J(f ) of the graph H obtained fromG as follows: for any edge has positive entries. For the second subsystem, it is easy to see that J 2 is the matrix J(f ) of the subgraph K of G corresponding to vertices in V 2 . Hence, L = −J T 2 is a Laplacian matrix, and by Proposition 9 its nullspace is generated by vectors π 1 , . . . , π k , where G 1 , . . . , G k are the terminal components in K (same as the terminal components with no excretions in G), and each vector π i has strictly positive entries corresponding to vertices in G i and is 0 elsewhere. Choosing each π i so that its entries sum to 1, the restriction of π i to G i is the stationary distribution of the corresponding Markov chain restricted to G i , namely the Markov chain whose generator matrix is the transpose of the submatrix of J(f ) corresponding to G i .
These remarks, together with standard tools of analysis of linear dynamical systems, lead to the following proposition: Example 4.1. In this example a simulation is used to verify the calculated equilibrium found using Proposition 12. Here we use the RCT network shown in 5 and calculate the equilibrium. Because every vertex of the RCT is connected to X we have thatẋ = J(f )x + φ with, The vector of fluxes and initial metabolite values were randomized to obtain Using these values the equilibrium was calculated to bex = (0.6354, 0.0824, 1.1040, 2.6211, 0.2015, 1.4122). The RCT network with the randomized initial values was simulated for 50 hours and the simulation results closely matchedx. The simulation results for the first 25 hours are shown in figure 6.

4.2.
Special LIFE systems. In this section, we focus on special LIFE systems. We first recall some interesting results from [16] valid under Assumption (B), and then we exploit them together with the spectral properties of J(f ) described in Sect. 4.1 in order to fully characterize equilibria and convergence of special LIFE systems under Assumption (C). We refer the reader to [14] for some results valid under great generality, with assumptions less restrictive than Assumption (B). As for linear systems, we first notice that weakly connected components correspond to subsystems which have no influence on each other, and hence can be studied separately; if there are any weakly connected components with no intake nor excretion, they are subsystems with constant total mass.

Proposition 13. ([16, Theorem 6]) Consider the special LIFE system under Assumption (B) with no intakes nor excretions. The following properties hold:
• The total mass of the system m = v∈V x v is constant in time.
• From any positive initial condition x(0), the system tends to the equilibrium set.

• Moreover, if there is a unique terminal component, then there exists a unique equilibrium with positive entries with the same mass as the initial mass, and the system converges to it.
For the general case with intakes and/or excretions, the following result holds on the asymptotic behavior of the dynamics. A caveat reported in [16] is that the second item in Proposition 14 does not rule out the possibility to have non-periodic oscillatory trajectories that approach the equilibrium set lying outside it, and approaching it rotating infinitely many times; this can happen in the case of a connected compact equilibrium set.
Remark 2. Recall that, by Proposition 2, if there exists a vertex v ∈ V such that there is a path from I to v and no path from v to X, then there exists no equilibrium. By the first item of Proposition 14, this further implies that all trajectories are unbounded.
The following results concern the existence and uniqueness of equilibria. • If there exists an equilibrium with positive entries, and if there exists a path from all v ∈ V to X, then the equilibrium is unique.
Notice that in the first item, the 'only if' part is true only when we require existence of equilibria with positive entries for completely arbitrary intakes: arbitrary intake set I and arbitrary positive values of the corresponding fluxes; this strong condition is not necessary to have an equilibrium with positive entries for a given set I and a given value of the intake fluxes. Also notice that this statement considers fixed fluxes f e , differently from Proposition 5, where x is fixed and fluxes f e are allowed to vary (except for intake fluxes); the crucial difference is that the product H e (x)f e can be made arbitrarily large in the context of Proposition 5, even in the case where H e (x) is bounded.
A remark about the second item is that, by Proposition 14, the existence and uniqueness of the equilibrium with positive entries further implies that all trajectories converge to such equilibrium.
To apply the first item of Proposition 14 or the second item of Proposition 15, one needs to already have the knowledge about existence of an equilibrium with positive entries. This happens for example in the case where one starts by fixing a desired equilibrium with positive entriesx, and then applies the extreme pathways technique in order to design suitable fluxes ensuring thatx is an equilibrium of the system. Under some assumptions on the graph, the above propositions then ensure uniqueness of the equilibrium, and its global asymptotic stability.

Special LIFE systems under Assumption (C).
In the remainder of this section, we consider Assumption (C). In this case, recall that the dynamics can be re-written asẋ = J(f )h(x) + φ. Different from linear systems, h(x) can contain non-linearities, but the matrix J(f ) has the same definition as for linear systems, so that its spectral properties are described by Propositions 9 and 11. Also notice that equilibria, i.e. solutions of J(f )h(x) = −φ, can be found by solving J(f )h = −φ, where h is an unknown vector in (R + ) n , and then solving h(x) = h. The latter is equivalent to solving H vi (x vi ) = h i , for i = 1, . . . , n.
We start by studying the case with no intakes nor excretions. In this case, by Proposition 9, J(f )h = 0 means h ∈ N (J(f )), where the nullspace N (J(f )) is generated by π 1 , . . . , π k , associated with terminal components G 1 , . . . , G k . Without loss of generality, we re-label vertices so that the first n 1 vertices encompass the first terminal component G 1 , the following n 2 vertices encompass the second terminal component G 2 , and so on, up to the last n k vertices encompass the last terminal component G k , and finally the remaining vertices are not in any terminal component (say there are n 0 of them). Denote the corresponding subblocks of vector h as h (1) , . . . , h (k) and h (0) for the non-terminal ones, and similarly define x (1) , . . . , x (k) and x (0) for vector x. By Proposition 9, the nullspace of J(f ) is N (J(f )) = {h ∈ (R + ) n : h (j) = α jπj , α j ∈ R and h (0) = 0 , for j ∈ {1, . . . , k}}, whereπ i is the restriction of π to the component G i , i.e. is a vector of size n i with strictly positive entries, representing the stationary distribution of the Markov chain whose generator is the transpose of the restriction of J(f ) to G i . Now we need to characterize the set of equilibria with positive entriesX := {x ∈ (R + ) n such that h(x) ∈ N (J(f ))}.
Recall that H v (x v ) are strictly increasing functions, being 0 when x v = 0. Denote by R v the range of H v (for x v ≥ 0), notice that either R v = [0, h max ), or R v = [0, +∞). Now denote by H j the set of vectors α jπj such that α j ≥ 0 and [α jπj ] v ∈ R v for all v in G j . Then denote by H −1 j h (j) the vector obtained from h (j) ∈ H j by applying entry-wise the inverse functions H −1 v . Finally we obtain X = {x ∈ (R + ) n such that x (j) = H −1 j h (j) , h (j) ∈ H j and h (0) = 0}. Having characterized the set of equilibria with positive entries, now recall that Proposition 13 applies, and trajectories remain bounded, with total mass constant in time, and approach the above-described equilibrium set. We focus on the case with intakes and/or excretions. The definition of the reduced system and the partitioning in two subsystems introduced for linear systems applies also to special LIFE systems under Assumption (C), with the only difference that now the dynamics are nonlinear: vectors h 1 (x 1 ) and h 2 (x 2 ) having replaced x 1 and x 2 in (16).
Hence, one can obtain the following analogous of Proposition 12.

Proposition 16. Consider a weakly connected special LIFE system satisfying Assumption (C). Denote by R v the range of the function H
, then from any initial condition, the reduced system converges to its unique equilibriumx 1 defined by [ . Otherwise, the system has no equilibrium, and lim t→∞ x(t) = +∞.
In the case where x 1 (t) converges tox 1 , if there are some vertices not connected to X, then: • If there exists a terminal component G T with no excretions such that there is a path from I to G T , then the mass in G T grows unbounded, and hence also lim t→∞ x 2 (t) = +∞ and lim t→∞ x(t) = +∞; • If for all terminal components with no excretions G 1 , . . . , G k there is no path from I to G i , then the mass of the system remains bounded, and moreover entries of x 2 (t) corresponding to non-terminal components converge to zero, while the restriction ofx 2 to a terminal component G i approaches the equilibrium set constructed as follows. The restriction of J 2 to vertices in G i has a nullspace generated by a single positive vectorπ i ; let H i denote the subset of such nullspace given by vectors h = απ i such that h v ∈ R v for all vertices v of G i ; the equilibrium set is given by x's such that there exist h ∈ H i verifying has positive entries; defineh to be the unique solution to J 1 h 1 = −φ 1 . Notice thath has all positive entries. If all entries ofh are within the range of the corresponding function H v (x v ), then we have a unique equilibrium with positive entriesx 1 obtained as in the statement of the proposition, and hence by Proposition 14 x 1 (t) converges tō x. Otherwise, there exists no equilibrium, and hence, by Proposition 14, the mass of system grows indefinitely.
If there is a terminal component with no excretions but connected to some intakes, then by Proposition 12 there is no equilibrium, and hence by Proposition 14 the mass of system grows indefinitely. If all terminal components with no excretions are not connected to intakes, then the equilibrium set is obtained from the properties of the nullspace of J 2 , given by Proposition 9. Moreover, Proposition 14 ensures that trajectories remain bounded and approach the equilibrium set.

4.3.
Zero-deficiency theory. In this section, we shortly recall the zero-deficiency theory, and compare it with our results on equilibria of LIFE systems. Zerodeficiency theory arises in the literature on chemical reaction networks, a seminal paper is [8]. Our short overview is based on [7].
A free closed chemical reaction network with m reactions between p complexes involving n species can be described byẋ = SΓR(x) where S ∈ M n×p , Γ ∈ M p×m is the incidence matrix of the network, and R(x) a column vector of size m. The deficiency of the chemical reaction network is the difference of the dimensions of the nullspaces of SΓ and Γ: δ = dim ker(SΓ) − dim ker(Γ), which is equivalent to δ = rank(Γ) − rank(SΓ) and to δ = p − ℓ − rank(SΓ), where ℓ is the number of weakly connected components. The first equivalence is due to rank-nullity Theorem, and the second to the fact that rank(Γ) = p − ℓ, since Γ is an incidence matrix (Proposition 4.3 of [2]). The reaction rate vector R(x) is governed by the 'massaction kinetics': R(x) = KΨ(x), where K ∈ M m×p and Ψ(x) is a vector of size p, defined as follows: K ej = k e if complex j is the reactant complex of reaction e, K ij = 0 otherwise; Ψ i (x) = n j=1 x sji j . Hence, the dynamics can be equivalently re-written asẋ = SΓKΨ(x). The following results hold:

Proposition 17 (Zero-deficiency Theorem). Consider a free closed chemical reaction network with mass-action kinetics. If its deficiency is zero, then: there exists an equilibrium with strictly positive entries if and only if the system is weakly reversible (i.e. each weakly connected component is also strongly connected).
Moreover, this strictly positive equilibrium is unique in each stoichiometric class (i.e. each weakly connected component has a space of equilibria which has dimension one, so that its equilibrium is unique up to a multiplicative constant representing the total mass in the component), and it is locally asymptotically stable.
Notice the close resemblance with the results for linear LIFE systems in the case with no intakes nor excretions, described in Propositions 9 and 10. It turns out that a particular class of closed free chemical reaction networks with mass-action kinetics and zero deficiency exactly coincides with linear LIFE systems with no intakes nor excretions. Indeed, let p = n, let S be the identity matrix of size n, let K be defined by K ej = f e if edge e is of the form e = (v j , w) and K ej = 0 otherwise; moreover choose exponents s ji to be zero and ones so that Ψ(x) = x. This gives exactly a linear LIFE system with no intakes nor excretions. It is immediate to see that the deficiency is zero, since S = I implies rank(SΓ) = rank(Γ).
In the more general case, free closed chemical reaction networks with mass-action kinetics and zero deficiency are a class of LIFE systems with different assumptions than those considered in this paper, e.g., they might not even satisfy Assumption (A), due to the presence of a rectangular matrix S left-multiplying Γ, and often do not satisfy Assumption (B), due to the dependence of entries of Ψ(x) on various entries of x.

5.
Conclusion. For general LIFE systems (Assumption (A)), we show the existence of positive solutions. Stability criteria have been shown for the structure of a graph associated to LIFE systems. We show that further conclusions cannot be drawn in the general case for LIFE systems because arbitrary dynamics may be defined.
For the problem of understanding equilibria for a fix set of metabolites, we show that the rank of the stoichiometric matrix (and the nullspace of this matrix) are determined by structural properties of the graph associated to the system, discuss the effect of the associated graph having intake vertices and excretion vertices on the rank of the stoichiometric matrix, and give necessary conditions on the structure of the graph associated to LIFE systems for the existence of equilibria. More biologically relevant equilibria of LIFE systems are those with all positive metabolites and fluxes and we prove that the extreme pathways method for finding a positive basis describes all such equilibria.
The field of network flows contributes to LIFE systems by way of the min cut max flow theorem. We show the capacity of edges in network flows relates to saturation of functions corresponding to edges of LIFE systems. The method of extreme pathways to calculate a positive basis is proven to include the entire intersection of the nullspace with the positive orthant. This basis is essential to describing equilibria of LIFE systems.
Equilibria and asymptotic behavior of metabolites for fixed fluxes are studied. We show that under stricter assumptions (Assumption (C)) we determine the eigenvalues of the jacobian of the LIFE system as well as structure of the graph which admit certain equilibria. We analyze the case with intakes and excretions versus the cases without intakes or excretions.
Equilibria for LIFE systems with terminal components exist, but these equilibria are dependent upon initial mass of the system. Furthermore, conditions are given for LIFE systems to tend to equilibrium from non-negative initial conditions. LIFE systems under Assumption (B) have bounded solutions if and only if there exits a non-negative equilibrium. Then for Assumption (B), we give conditions on the structure of the associated graph for which nontrivial equilibria exist.
We show that LIFE systems with no intakes or excretions have "zero deficiency." The rank of the stoichiometric matrix (defined in this work) gives information about the structure of the associated graph, specifically the connectivity of the graph and the existence of strongly connected components. Zero deficiency theory also provides information about the existence of equilibria with respect to this structure.
The structural conditions of graphs discussed in this work have implications on metabolic networks. We have included results about networks which are considered non-biological for the sake of completeness. These results allow a clear picture of the structure of metabolic systems which are capable of admitting a biologically relevant equilibrium. With this clear picture, we are able to determine metabolic networks for which it is most advantageous to analyze using Linear-In-Flux-Expressions.
Appendix. This Appendix contains the calculations for examples 3.1 and3.2. Example 3.1 illustrates finding the extreme pathways of RCT using the algorithm described in Section 3.2. Example 3.2 compares a standard basis and positive basis for RCT and shows that when only the positive orthant is considered they describe the same space.
Example 3.1 continued. As shown previously, the stoichiometric matrix for the RCT is given by: First let V 0 be the m× m identity matrix where m is the number of columns in S. Then define C 0 as follows, V 0 · S T = C 0 . To better illustrate the row combinations in the algorithm a column is appended to C 0 which labels the rows.
We have The next step is to identify any columns of C 0 that have no sources or sinks. For this example there are two columns, column 4 and column 5. Then we modify both V 0 and C 0 by first copying all the rows with a zero in the entry of the first identified column (column 4), then of the remaining rows add all possible positive combinations of two rows which produce a zero in this column. Note that V 0 is 10 × 10 and C 0 is 10 × 6.
For the stoichiometric matrix S(x) the basis for the null space was found. The span of the four basis vectors is shown below, To find the intersection of this span with the positive orthant, three conditions on the a i must hold.
We refer to the intersection of this span with the positive orthant as C. We show that under the conditions given, B = C.
First it is shown that B ⊂ C by showing that for an arbitrary set of b i ≥ 0, a i 's can be chosen to reach the same vector. Using the following substitutions for a i it is clear that if b i ≥ 0 the inequalities of (22) are satisfied This shows that any vector in B can be represented by vectors in C i.e. B ⊂ C.
Now it is shown that C ⊂ B. For arbitrary a's which satisfy (22), the following substitutions for b are used and the conditions of (20) are checked.
To insure that (20) is satisfied the following inequalities must hold, If choices for b 4 and b 6 can be found which satisfy these inequalities then we will have C ⊂ B.
We have the following two conditions:   (22) we have a 1 x 6 ≥ a 3 x 3 + a 4 x 2 , which means that (25) is also satisfied.  (22) we have a 1 x 6 ≥ a 3 x 3 + a 4 x 2 , which means that (25) is also satisfied.
Since the a ′ s were arbitrary and b ′ s are found which satisfy (20) this gives us that C ⊂ B as desired.