Exploiting stoichiometric redundancies for computational efficiency and network reduction

Abstract

Abstract

Analysis of metabolic networks typically begins with construction of the stoichiometry matrix, which characterizes the network topology. This matrix provides, via the balance equation, a description of the potential steady-state flow distribution. This paper begins with the observation that the balance equation depends only on the structure of linear redundancies in the network, and so can be stated in a succinct manner, leading to computational efficiencies in steady-state analysis. This alternative description of steady-state behaviour is then used to provide a novel method for network reduction, which complements existing algorithms for describing intracellular networks in terms of input-output macro-reactions (to facilitate bioprocess optimization and control). Finally, it is demonstrated that this novel reduction method can be used to address elementary mode analysis of large networks: the modes supported by a reduced network can capture the input-output modes of a metabolic module with significantly reduced computational effort.

1 Introduction

The stoichiometry matrix provides a conciserepresentation of the structure of a chemical reaction network. This matrix is the primary system description when analyzing metabolic networks, complemented by constraints on reaction rates or full reaction kinetics when these are available.

The analysis in this paper begins with the stoichiometry matrix, N. This matrix is invariably not full rank; redundancies in the matrix reflect dependancies within the network. Linearly dependent rows of N correspond to chemical species whose dynamic behaviour is dependent on other species in the network. The reduction in system complexity through elimination of these redundant species is standard in stoichiometric analysis [3, 11], and the resulting computational advantages are well recognized.

This paper shows that a description of the redundancies among the columns of N can be harnessed to improve analysis of steady-state behaviour in a number of ways. We begin with the observation that the steady-state balance condition depends only on the form of the redundancies within N. Consequently, a great deal of the content of the stoichiometry matrix plays no role in any investigation of steady-state behaviour. Dropping this unnecessary content results in streamlined steady-state analysis. The techniques of metabolic control analysis and metabolic balancing are taken asexamples. The relationship between the stoichiometry matrix and steady-state system behaviour is not one-to-one. This observation allows us to define classes of networks that are equivalent from the steady-state point of view. We then relax this equivalence notion to produce network reductions that respect the stoichiometric structure of the original network. This novel reduction procedure complements similar methods that describe metabolic networks in terms of macro-reactions [2, 10]. These methods facilitate the application of bioprocess optimization and control techniques to models of intracellular metabolism. In contrast with previously published approaches, our reduction strategy allows for targeted elimination of specific reactions, which can be used to generate simplified descriptions of the flux profiles within an input-output network module. In particular, this method can be used to significantly reduce the computational effort required to determine the set of elementary flux modes through network modules, complementing existing approaches to this problem [5, 8].

This paper is organized as follows. Section 2 describes how dependancies among the rows and columns of a stoichiometry matrix can be identified, and then used to formulate a factored stoichiometric matrix. Section 3 briefly outlines how this matrix factorization leads to computational efficiencies in two commonly used analysis techniques. In Section 4, this factored matrix formulation is used to characterize classes of networks that share identical steady-state behaviour. In Section 5, these classes of networks are generalized to include network reductions whose steady state behaviour is consistent with the original. It is demonstrated that flux mode analysis of a reduced network provides a meaningful description of flux through the original network. An E. coli core metabolism network is analyzed as an example. Concluding thoughts are contained in Section 6. Finally, an appendix details how the novel network reduction method described in this paper is, in certain cases, equivalent to the network reduction approach described in [2]. A preliminary description of some of these results appeared in [4].

2 Rank deficiencies

Consider a network of n chemical species involved in m reactions in a fixed volume. The concentrations of the species are the elements of the n-vector s. The reaction rates are collected into the m-vector rt, which depends on the species concentrations and on a vector of parameters pr. The network topology is described by the n × m stoichiometry matrix N, whose i, j-th element indicates the net number of molecules of species i produced in reaction j (negative values indicate consumption). The system dynamics are described at each time t by

$\frac{d}{dt}\text{s}(t)={\text{N}}_{\text{V}}(\text{s}(t),)p).$(1)

Analysis of this system can be simplified by exploiting the dependencies among the rows and columns of N. Metabolic networks typically have highly redundant stoichiometry matrices. (An example is the metabolic map from Escherichia coli from [12] which has a 770×931 stoichiometry matrix of rank 733.)

Linear dependencies among the rows of N correspond to structural conservations, which most often occur as conserved moieties. Linear dependencies among the columns of N correspond to steady-state flux distributions. We next formalize these stoichometric features.

2.1 Deficiencies in row rank

This subsection reviews the standard procedure given in [11] for treating the structural conservations characterized by linearly dependent rows of N. Let r denote the rank of N and re-index the species so that the first r rows of N are linearly independent. Identify the first r species as the independent species, denoted _si, and the remaining n - r as the dependent species, denoted _sd, so that $\text{s}=[{\text{s}}_i^T,{\text{s}}_d^T{]}^T$. Let NR denote the submatrix of N consisting of the first r rows. We can then write N = LiNR, where the matrix Li, referred to as the row link matrix, has the form

$L=\left[\begin{array}{c}\hfill I_r\hfill \\ \hfill L_0\hfill \end{array}\right].$(2)

This row link matrix can be constructed from a matrix G whose rows form a basis for the left nullspace of N. Such a matrix can be chosen of the form

$G=[-L_0{I}_{n-r}],$

which is referred to as the conservation matrix. The rows of this matrix correspond directly to combinations of species concentrations that are conserved. Computational approaches to finding this conservation matrix are reviewed in [13] (see also [18]).

Following [11], note that

$\frac{d}{dt}\text{s}(t)=\frac{d}{dt}\left[\begin{array}{c}\hfill {\text{s}}_i(t)\hfill \\ \hfill {\text{s}}_d(t)\hfill \end{array}\right]=L{N}_R\text{v}(\text{s}(t)p).$

Equation (2) then gives

$\frac{d}{dt}{\text{s}}_d(t)=\frac{d}{dt}{L}_0{\text{s}}_i(t)\text{for}\text{all}t\ge 0.$

Integrating gives an explicit dependence among the species concentrations: ${\text{s}}_d\left(t\right)={\text{L}}_0{\text{s}}_i\left(t\right)+\tilde{\text{T}}$ for all time, where $\tilde{\text{T}}={\text{s}}_d\left(0\right)-{\text{L}}_0{\text{s}}_i\left(0\right)$. Finally, concatenating $\tilde{\text{T}}$ with the r-vector of zeros 0 _r , define $\text{T}=[0_r^T,{\tilde{\text{T}}}^T{]}^T$ so that

$\text{s}(t)=L{\text{s}}_i(t)+T.$(3)

As a consequence of this decomposition, dynamic analysis can be restricted to a reduced version of equation (1), namely

$\frac{d}{dt}{\text{s}}_i(t)={\text{N}}_{\text{R}}\text{v}(L{\text{s}}_i(t)+T\text{,}p)\text{.}$(4)

The partitioning of species into independent and dependent classes is not unique; specific methods for partitioning are discussed in [3].

2.2 Deficiencies in column rank

Redundancy in the columns of N can be used to arrive at a factored description of system behaviour, as follows.

Mirroring the analysis in the previous subsection, label the reactions so that the first m - r columns of N are linearly dependent on the remaining r, and partition the reaction rates accordingly into m - r independent rates, denoted _rti, and r dependent rates, denoted _rtd, so $\text{v}=[{\text{v}}_i^T,{\text{v}}_d^T{]}^T$.

In analogy to the construction of the row link matrix, let NC denote the submatrix of N consisting of the last r columns, from which N can be recovered as N = NCPP, where the column link matrix P is of the form

$P=[\begin{array}{cc}\hfill P_0\hfill & \hfill I_r\hfill \end{array}].$

The column link matrix, which is dual to the row link matrix, can be determined through construction of a matrix K of the form

$\text{K}=\left[\begin{array}{c}\hfill I_{m-r}\hfill \\ \hfill -P_0\hfill \end{array}\right],$

whose columns span the (right) nullspace of N. The linearly independent columns of K correspond to steady state flows within the network (irrespective of any reversibility constraints).

To arrive at a factored description of system dynamics, write

$\begin{array}{l}\frac{d}{dt}s(t)=\text{Nv}(\text{s}(t),p)\\ ={\text{N}}^{\text{C}}P\text{v}(\text{s}(t),\text{p})\\ ={\text{N}}^{\text{C}}[P_{\text{0}}{I}_r]\text{v}(\text{s}(t),p)\end{array}$(5)

This description is consistent with the standard method by which redundancy in columns is exploited. As described in, e.g. [3], the system flux vector J (i.e. the steady state reaction rate vector) can be partitioned into dependent and independent components $\text{J}=[{\text{J}}_i^T,{\text{J}}_d^T{]}^T$ corresponding to the partitioning of the reaction rates described above. The dependence is given explicitly by

$J=K{J}_i=\left[\begin{array}{l}{I}_{m-r}\\ -{\text{P}}_0\end{array}\right]{\text{J}}_i.$(6)

(In [3] the submatrix -_PP0 is referred to as _K0.) Equation (6) can be derived from Equation (5) as follows. At steady state, the fact that NC has full column rank gives _{PP0 rti} + _rtd = 0. Consequently _rtd = - _{PP0 rti}, and so rt = K _rti . Equation (6) provides a dependence among the steady-state reaction rates that is dual to the dependence among the species concentrations in Equation (3).

The partitioning of independent and dependent reaction rates is not unique. A procedure for choosing independent reaction rates as the entry and exit points from the network is outlined in [20].

2.3 Complete stoichiometric reduction

The previous two subsections describe complementary system decompositions based on linear dependencies. These can be combined to arrive at a fully factored description of system dynamics:

$\begin{array}{l}\frac{d}{dt}\text{s}(t)=L{N}_R^{C}P\text{v}(\text{s},\text{p})\\ =\left[\begin{array}{l}{\text{I}}_r\\ {\text{L}}_0\end{array}\right]{\text{N}}_{\text{R}}^{\text{C}}[{\text{P}}_0{\text{I}}_r]\text{v}(\text{s,}\text{p})\end{array}$(7)

where NRC is the invertible upper right r × r submatrix of N. In the following discussion, the submatrix NRC will be referred to as a stoichiometric core of the system. (This factorization of N is complementary to the singular value decomposition, discussed in this context in [1].)

Restricting attention to the dynamics of the independent species, it follows from equation (7) that

$\frac{d}{dt}{\text{s}}_i(t)={\text{N}}_{\text{R}}^{\text{C}}\text{Pv}(\text{L}{\text{s}}_{\text{i}}\text{+T,}\text{p})$(8)

Because the stoichiometric core NRC is invertible, a concise description of steady state behaviour is:

$0=\text{Pv}(\text{L}{\text{s}}_{\text{i}}(t)+\text{T,}\text{p})$(9)

Consequently, when analyzing steady state behaviour, matrix storage cost can be reduced by retaining only L and P, rather than the full stoichiometry matrix N. (In the case of the E. coli metabolic map from [12] this results in a reduction of over 75% – more than half a million matrix entries.) Of course, the storage cost for N, _Li0 and _PP0 can typically be reduced by exploiting their sparsity, and so the storage savings afforded by this factorization may not be significant.

The characterization of steady state in Equation (9) indicates that the full rank stoichiometric core NRC plays no role in steady state analysis. This observation was used by Wagner to formulate the “nullspace” approach to identification of elementary flux modes [19]. The following section briefly demonstrates computational efficiencies afforded by this observation in two commonly used analytic techniques.

3 Computational efficiencies

3.1 Metabolic control analysis

Metabolic Control Analysis (MCA) provides a theory of local parametric sensitivity analysis that takes advantage of stoichiometric structure. The direct approach [11] begins with the steady state condition from equation (4). Differentiating with respect to the parameter pr yields (unscaled) concentration sensitivity coefficients (called response coefficients) of the form

$\frac{d\text{s}}{d\text{p}}=-\text{L}{\left({\text{N}}_{\text{R}}\frac{\partial v}{\partial s}\text{L}\right)}^{-1}{\text{N}}_{\text{R}}\frac{\partial \text{v}}{\partial \text{p}}$(10)

A streamlined version of this formula can be derived by differentiating equation (9), resulting in

$0=\left[\text{P}\frac{\partial \text{v}}{\partial \text{s}}\text{L}\frac{d{\text{s}}_i}{d\text{p}}+\text{P}\frac{\partial \text{v}}{\partial \text{p}}\right].$

When the matrix $\text{P}\frac{\partial \text{v}}{\partial \text{s}}\text{L}$ is presumed invertible (by the standard assumption of asymptotic stability), one arrives at the (unscaled) response coefficients in the form

$\frac{d\text{s}}{d\text{p}}=-\text{L}{\left(\text{P}\frac{\partial \text{v}}{\partial \text{s}}\text{L}\right)}^{-1}\text{P}\frac{\partial v}{\partial \text{p}}.$(11)

The (unscaled) concentration control coefficients are then

${\text{C}}^{\text{s}}=-\text{L}{\left(\text{P}\frac{\partial \text{v}}{\partial \text{s}}\text{L}\right)}^{-1}\text{P}.$

Flux response coefficients and flux control coefficients can be derived accordingly.

Formula (11) is similar to formula (10), but involves the column link matrix P in place of the partially reduced stoichiometry matrix N_R. The consequence of this streamlined description is a reduced computational effort and a clear role for stoichiometric redundancy in the form of the sensitivity coefficients.

3.2 Metabolic flux analysis: Metabolic balancing

In metabolic balancing [7, 16], a metabolic network is assumed to be at steady state and a subset of the fluxes are measured. One then attempts to infer the remaining fluxes. Following [7], partition the reaction vector rt into unknown rates, denoted ^rtn, and known rates, denoted ^rtb. The stoichiometry matrix is likewise partitioned into N ⁿ and N ^b and the steady state balance equation (from equation (1)) can be written as

${\text{N}}^n{\text{v}}^n=-{\text{N}}^{\text{b}}{\text{v}}^{\text{b}}.$(12)

One then attempts to solve for ^rtn. In the simplest case N ⁿ is invertible and ^rtn is found exactly. More typically, one can make use of the pseudo-inverse of N ⁿ . The paper [7] contains a discussion of the cases in which the given measurements are insufficient or are inconsistent (in which cases the system (12) is under- or over-determined, respectively).

If the known and unknown rates are each partitioned into independent and dependent components (with the columns of N re-ordered accordingly) one arrives at a decomposition of rt into four components:

$v=\left(\begin{array}{c}\hfill {\text{v}}_i^n\hfill \\ \hfill {\text{v}}_i^b\hfill \\ \hfill {\text{v}}_d^n\hfill \\ \hfill {\text{v}}_d^b\hfill \end{array}\right).$

The steady state equation (9) is simply

$\text{Pv}=[{\text{P}}_0{\text{I}}_r]\text{v}=0.$

Partitioning ${\text{P}}_0=\left[{\text{P}}_0^n{\text{P}}_0^b\right]$ allows us to write thisbalance equation as

$\left[{\text{P}}_0^n{\text{P}}_0^b\begin{array}{cc}\hfill \text{I}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \text{I}\hfill \end{array}\right]\left(\begin{array}{c}\hfill {\text{v}}_i^n\hfill \\ \hfill {\text{v}}_i^b\hfill \\ \hfill {\text{v}}_d^n\hfill \\ \hfill {\text{v}}_d^b\hfill \end{array}\right)=0.$

where the blocks in the identity matrix have dimensions corresponding to ${\text{v}}_d^n$ and ${\text{v}}_d^b$. Isolating the unknown rates leads to a restatement of equation (12):

$\left[{\text{P}}_0^n\begin{array}{c}\hfill \text{I}\hfill \\ \hfill 0\hfill \end{array}\right]\left(\begin{array}{c}\hfill {\text{v}}_i^n\hfill \\ \hfill {\text{v}}_d^n\hfill \end{array}\right)=-\left[{\text{P}}_0^b\begin{array}{c}\hfill 0\hfill \\ \hfill \text{I}\hfill \end{array}\right]\left(\begin{array}{c}\hfill {\text{v}}_i^b\hfill \\ \hfill {\text{v}}_i^b\hfill \end{array}\right).$

One is now faced with inverting $\left[\begin{array}{cc}\hfill {\text{P}}_0^n& \hfill \end{array}\begin{array}{c}\hfill \text{I}\\ \hfill 0\end{array}\right]$ as opposed to the original task of inverting N ⁿ . Although these matrices are of the same size, the former will typically be easier to deal with. In particular, the inverse of $\left[\begin{array}{cc}\hfill {\text{P}}_0^n& \hfill \end{array}\begin{array}{c}\hfill \text{I}\\ \hfill 0\end{array}\right]$ (if it exists) can be expressed directly in terms of the inverse of a submatrix of ${\text{P}}_0^n$, since in general

${\left[\begin{array}{cc}\hfill \text{A}\hfill & \hfill \text{I}\hfill \\ \hfill \text{B}\hfill & \hfill 0\hfill \end{array}\right]}^{-1}=\left[\begin{array}{cc}\hfill 0\hfill & \hfill {\text{B}}^{-1}\hfill \\ \hfill \text{I}\hfill & \hfill -\text{A}{\text{B}}^{-1}\hfill \end{array}\right]$

if B is invertible. Furthermore, the same submatrix can be used to simplify computation of the pseudo-inverse of $\left[\begin{array}{cc}\hfill {\text{P}}_0^n& \hfill \end{array}\begin{array}{c}\hfill \text{I}\\ \hfill 0\end{array}\right]$ through the use of the Schur complement [15]. The pseudo-inverse of

$\left[\begin{array}{cc}\hfill \text{A}\hfill & \hfill \text{I}\hfill \\ \hfill \text{B}\hfill & \hfill 0\hfill \end{array}\right]$

is

${\left({\left[\begin{array}{cc}\hfill \text{A}\hfill & \hfill \text{I}\hfill \\ \hfill \text{B}\hfill & \hfill 0\hfill \end{array}\right]}^T\left[\begin{array}{cc}\hfill \text{A}\hfill & \hfill \text{I}\hfill \\ \hfill \text{B}\hfill & \hfill 0\hfill \end{array}\right]\right)}^{-1}\left[\begin{array}{cc}\hfill \text{A}\hfill & \hfill \text{I}\hfill \\ \hfill \text{B}\hfill & \hfill 0\hfill \end{array}\right]$

where

$\begin{array}{ccc}\hfill & \hfill \hfill & {\left({\left[\begin{array}{cc}\hfill \text{A}& \hfill \text{I}\\ \hfill \text{B}& \hfill 0\end{array}\right]}^T\left[\begin{array}{cc}\hfill \text{A}& \hfill \text{I}\\ \hfill \text{B}& \hfill 0\end{array}\right]\right)}^{-1}\hfill \\ \hfill & \hfill \hfill & =\left[\begin{array}{cc}\hfill ({\text{B}}^T\text{B}{)}^{-1}\hfill & \hfill -({\text{B}}^T\text{B}{)}^{-1}{\text{A}}^T\hfill \\ \hfill -\text{A}({\text{B}}^T\text{B}{)}^{-1}\hfill & \hfill \text{I}+\text{A}({\text{B}}^T\text{B}{)}^{-1}{\text{A}}^T\hfill \end{array}\right].\hfill \end{array}$

The remainder of this paper addresses a novel analytic approach facilitated by the factored description of system behaviour in Equation (7).

4 Steady-state-equivalent networks

Returning to the steady state balance equation (9), recall that this condition depends only on the left and right nullspaces of the stoichiometry matrix N. It follows that any network whose stoichiometry matrix shares these nullspaces will be described by identical steady state conditions. This observations allows the set of all n × m stoichiometry matrices to be partitioned into equivalence classes that share the same steady-state behaviour: two stoichiometry matrices are members of the same equivalence class when they share the same left and right nullspaces.

As a simple example, consider the three-species network shown in Fig. 1A.

Here

$N=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill \end{array}\right]$

which has rank two. Following the procedure in the previous section, we take the stoichiometric core NRC as the upper right 2 × 2 submatrix, resulting in link matrices

$\text{L}=\left[\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 1\hfill \end{array}\right]\text{and}\text{P}=\left[\begin{array}{ccc}\hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right]$

Choosing any invertible 2 × 2 matrix as an alternative stoichiometric core leads to a steady-state equivalent network. For example

${\text{N}}_{\text{R}}^{\text{C}}=\left[\begin{array}{cc}\hfill -1\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill -1\hfill \end{array}\right]\text{gives}\text{N}=\text{L}{\text{N}}_{\text{R}}^{\text{C}}\text{P}=\left[\begin{array}{ccc}\hfill 2\hfill & \hfill -1\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill \end{array}\right]$

(Figure 1B) while

${\text{N}}_{\text{R}}^{\text{C}}=\left[\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 2\hfill & \hfill -1\hfill \end{array}\right]\text{gives}\text{N}=\text{L}{\text{N}}_{\text{R}}^{\text{C}}\text{P}=\left[\begin{array}{ccc}\hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill \end{array}\right]$

(Figure 1C).

The members of this equivalence class are networks composed of three species (S ₁, S ₂, S ₃) involved in three reactions (rates v ₁, v ₂, v ₃) for which (i) a conservation [S ₁] - [S ₂] + [S ₃] = T holds for all time, and (ii) steady-state occurs when v ₁ = v ₂ = v ₃. All such networks can be generated by the choice of different stoichiometric cores, and the results of a steady state analysis performed for any one representative of this equivalence class will apply equally to every member of the class.

As another example, consider the network shown in Fig. 2A, (from [6], with reactions re-indexed).

The stoichiometry matrix is

$\text{N}=\left[\begin{array}{cccccccccc}\hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{array}\right]$(13)

which has full row rank. The reaction indexing was chosen so that the last six columns form a full-rank submatrix, which can be taken as a stoichiometric core NRC. The kernel of N is spanned by the columns of

$\text{K}=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 2\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]$

from which the column link matrix can be constructed.

In [6] this network is used to illustrate the method of metabolic balancing. Suppose that the following flux measurements are available: v ₇ = 2, v ₉ = 2 and v ₁₀ = 1 (units arbitrary). It can then be verified that v ₁ = 1, v ₄ = 1, v ₆ = 1, v ₈ = 1, while v ₂, v ₃ and v ₅ are constrained only by the two conditions v ₂ = v ₃ = 1 - v ₅.

Consider now another member of the equivalence class represented by this network, generated by the alternative stoichiometric core

${\text{N}}_{\text{R}}^{\text{C}}=\left[\begin{array}{cccccc}\hfill -1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \end{array}\right]$

The steady-state behaviour of the corresponding alternative network (generated by multiplication by the column link matrix), shown in Fig. 2B, is identical to the original. For instance, as the reader can verify, the metabolic balancing result cited above (for the original) carry over precisely to this alternative network.

4.1 Relaxation of steady-state constraints

The structural conservations and steady-state flux distributions in a system are described by the row and column link matrices L and P, respectively. Each choice of an invertible stoichiometric core leads to a stoichiometry matrix composed of linearly independent combinations of the rows of L and the columns of P, and so all of the corresponding stoichiometry matrices share the same left- and right-nullspaces.

Alternatively, a stoichiometric core could be chosen that is not full rank (but is still of the appropriate dimension: r × r). In this case, the corresponding stoichiometry matrix has the same number of rows and columns as the original, but the dimension of the nullspaces will have increased (because the combinations of rows of L and the columns of P are not independent). Such a stoichiometry matrix characterizes a network in which the original constraints have been relaxed.

As an example, consider again the network shown in Fig. 2A. Recall that the original stoichiometric core is the upper-right 6 × 6 submatrix of the stoichiometry matrix (13). Consider an alternative core:

${\text{N}}_{\text{R}}^{\text{C}}=\left[\begin{array}{cccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]$

which has rank 4. The corresponding stoichiometry matrix characterizes the network shown in Fig. 3.

In the original network, the steady-state constraints on the reaction rates are:

$\begin{array}{l}{\upsilon }_6={\upsilon }_1{\upsilon }_5={\upsilon }_1-{\upsilon }_3{\upsilon }_7=2{\upsilon }_1+{\upsilon }_2-{\upsilon }_3\\ {\upsilon }_{10}={\upsilon }_1{\upsilon }_9={\upsilon }_1+{\upsilon }_1{\upsilon }_8=-{\upsilon }_2+{\upsilon }_3+{\upsilon }_4.\end{array}$

Figure 3 makes it clear that most of these constraints hold in the alternative network. The exceptions are that v ₅ is absent, and the two constraints on v ₉ and v ₁₀ have been merged into the single requirement that v ₉ = v ₁₀ + v ₄.

Arbitrary relaxation of steady-state constraints cannot be expected to yield meaningful networks. In contrast, we next present a strategy for targeted network reduction via specific choice of a non-invertible stoichiometric core.

5 Steady-state-consistent network reduction

The network in Fig. 3 is simpler than the original (Fig. 2A): there are two fewer species and one less reaction. By tailoring the choice of a non-invertible stoichiometric core, we can direct this simplification process to arrive at a meaningful reduction of the original network.

Previously published algorithms for network reduction have been motivated by the need to replace uncharacterized intracellular reactions with a network of macro-reactions that connect extracellular substrates to extracellular products [2, 10]. In this section, we describe targeted network reductions that can be reached by appropriate choice of low rank stoichiometric cores. As shown in the Appendix, this procedure recapitulates the technique of Haag et al. [2] in certain cases.

Low rank stoichiometric cores can be generated by zeroing selected columns of the original stoichiometric core, thus eliminating the corresponding reaction from the network. This procedure leads to a reduced network that is consistent with the original behaviour at steady state.

An algorithm for network reduction is as follows:

Determine the rank r of the given n × m stoichiometry matrix N.

Identify a set of r linearly independent columns of N; re-index the reactions so that these are the rightmost columns of N. From this set, select q < r columns to be eliminated in the network reduction.

Generate the row link matrix L and the column link matrix P.

Take the upper right r × r submatrix of N as the stoichiometric core NRC.

Generate a reduced rank core by replacing the q columns of NRC that have been targeted for elimination with columns of zeros.

Generate a reduced stoichiometry matrix as the product LiNRCPP.

The choice of reactions to be eliminated (step 2) will depend on the aim of the reduction. In some cases, there may be specific reactions whose elimination is preferred. In others, this choice may be dictated by desired characteristics of the reduced network. An example of this latter case is addressed in the next section.

As an example, consider once more the network in Fig. 2A. Recall that the original stoichiometric core is the upper right 6 × 6 submatrix of the stoichiometry matrix (13). Taking an alternative stoichiometric core resulting from zeroing the first column of NRC (which corresponds to reaction 5) leads to the alternative stoichiometry indicated in Fig. 4A.

Alternatively, elimination of reaction 6, by zeroing the second column of NRC, leads to the topology shown in Fig. 4B.

5.1 Exchange-equivalent steady-state-consistent network reductions

When the reduction procedure described in the previous subsection is applied to a network, each of the original reactions is treated in one of three ways:

Reactions targeted for elimination are eliminated.

Reactions corresponding to the retained columns of the stoichiometric core are unchanged in the reduced network.

All other reactions are subject to lumping.

Analysis of metabolic networks is often focused on input-output behaviour: the processing of substrates to end products. To arrive at a network reduction that is consistent with the input-output behaviour of the original network, the exchange reactions can be chosen as part of the stoichiometric core to be retained unchanged in the reduction. The resulting reduced network can be interpreted as a minimal input-output description of the original network.

As an example, starting with the network in Fig. 2, consider an alternative stoichiometric core generated by zeroing the first and second columns, retaining only the exchange fluxes (i.e. the last four columns). The resulting minimal input-output description is shown in Fig. 5.

All internal species have been eliminated. The fluxes have been reduced to a minimal number while maintaining stoichiometric consistency among the exchange species.

5.2 Flux mode analysis

One reason for carrying out network reduction is to reduce the number of elementary flux modes (EFMs) through the network. EFMs characterize the potential flux distributions within a network, and are useful in a range of analysis techniques [14]. The number of EFMs supported by a network grows rapidly with network size, limiting the use of EFM-based analysis for large networks. Recent efforts to alleviate this problem have focused on techniques for identifying EFMs through subnetworks [5, 8]. As discussed in [8], the choices made in isolating a subnetwork can impact the resulting flux distribution, which may not properly represent the behaviour of the subnetwork within the original network. The reduction technique introduced in this paper provides a means to reduce the number of EFMs in a subnetwork while retaining all exchange fluxes, so that the interface with surrounding network components is unchanged.

A significant drawback of the reduction procedure described here is that it does not directly account for constraints on reaction rates. EFMs are constructed with tools from convex analysis, which allow irreversibility constraints to be incorporated into the description of feasible flux distributions. When the reduction technique described above is applied to a network, irreversibility constraints on the eliminated reactions appear as inequality constraints on linear combinations of the reaction rates in the reduced network. For instance, in the reduction in Fig. 5, irreversibility constraints on the eliminated reactions 5 and 6 would result in the following constraints on the reactions in the reduced network:

$\begin{array}{l}{\upsilon }_5\ge 0\Rightarrow {\upsilon }_1-{\upsilon }_3\ge 0\\ {\upsilon }_6\ge 0\Rightarrow {\upsilon }_1\ge 0\end{array}$

The constraint v ₆ ≥ 0 corresponds directly to an irreversibility constraint on v ₁, which is retained in the reduced network. However, the constraint on v ₅ leads to a more complex inequality constraint. Consequently, the set of feasible flux distributions in the reduced network is not a cone, and so cannot be analyzed directly using the standard techniques for analysis of flux distributions.

This issue can be side-stepped by targeting only reversible reactions for elimination. In that case,standard elementary mode analysis can be applied to the reduced network (in which any irreversibility constraints will be unaltered).

As an example, consider again the procedure applied to the network in Fig. 2A to arrive at the reduced network in Fig. 5. The software package efmtool [17] was used to analyze three cases of irreversibility constraints on the original network:

If all reactions are reversible, then the original network supports 21 EFMs. The reduced network supports 16.

If the exchange reactions (7–10) are irreversible (oriented as in Fig. 2) and all internal reactions are reversible, then the original network supports 15 EFMs. The reduced network supports 11.

If reactions 5 and 6 are reversible and all other reactions are irreversible (oriented as in Fig. 2), then the original network supports 7 EFMs. The reduced network (fully irreversible in this case) supports only 5.

The EFMs for the third case (all reactions irreversible except for reactions 5 and 6) are illustrated in Fig. 6. The first five elementary modes are clearly comparable between the original and reduced networks. The sixth has no counterpart in the reduced network; it involves only internal reactions, and has been eliminated. The seventh elementary mode of the original network corresponds to a non-elementary flux mode in the reduced network. The input-output flux distributions can be described as follows:

$\begin{array}{l}\text{EFM}1:{\upsilon }_7\Rightarrow {\upsilon }_9\text{EFM2}:{\upsilon }_8\Rightarrow {\upsilon }_9\\ \text{EFM}4:2{\upsilon }_7\Rightarrow {\upsilon }_9\&{\upsilon }_{10}\text{EFM2}:2{\upsilon }_8\Rightarrow {\upsilon }_9\&{\upsilon }_{10}.\end{array}$

These input-output behaviours are effectively captured in the network reduction.

These examples demonstrate that when targeting internal reversible reactions for elimination, there will generally be considerable freedom in terms of the choice of reactions to eliminate. The number of reactions eliminated in the reduction is constrained by the rank of the stoichiometry matrix. In some cases (e.g. the first two cases above), only a subset of the internal reversible reactions can be eliminated; the choice of which subset could be dictated by further criteria. Alternatively, it may be that after targeting all internal reversible reactions for elimination, there will be an opportunity to eliminate further reactions (which would cause changes in input-output structure or rate constraints).

5.3 Example: core metabolism of E. coli

To better illustrate the potential of this reduction method for facilitating elementary mode analysis, we consider a more complex metabolic network: the core E. coli metabolism network described in [9]. This network is composed of 63 species involved in 77 reactions. The stoichiometry matrix has rank 57, indicating that there are 6 structural conservations within the system. Of the 77 reactions, 14 are reversible exchange reactions, and one (biomass production) can be taken as an irreversible exchange reaction. Of the remaining 62 internal reactions, 27 are irreversible and 35 are reversible. The network admits 2 295 967 elementary modes (as determined by efmtool [17]).

To generate a reduced version of the network, the species and reactions were indexed so that (i) the top 57 rows of the stoichiometry matrix have full rank, (ii) the 15 exchange reactions correspond to the rightmost columns of the stoichiometry matrix, (iii) the columns corresponding to the 35 reversible internal reactions appear among the rightmost 57 columns of the stoichiometry matrix, and (iv) the rightmost 57 columns of the stoichiometry matrix have full rank. The upper-right 57 × 57 submatrix of the stoichiometry matrix then forms a full rank stoichiometric core, from which the original stoichiometry matrix N can be recovered from the row and column link matrices L and P.

To generate a reduced network, an alternative stoichiometric core NRC was constructed as follows. The columns in the original core that correspond to reversible internal reactions (35 in all) were replaced by columns of zeros. The remaining columns of the original core (corresponding to either exchange reactions or irreversible internal reactions) were unchanged. A reduced stoichiometry matrix was generated as N _red = LiNRCPP. The reduced network is composed of 52 species involved in 42 reactions, of which 15 are exchange reactions and 27 are irreversible internal reactions. This network admits 18 480 elementary modes (fewer than 1% of the original). This set of flux modes provides a full characterization of the system’s steady-state input-output behaviour. As in the analysis in the previous section, the modes that are not retained in the reduction are either (i) internal cycles or (ii) input-output modes that are not elementary in the reduced network.

The seven irreversible internal reactions that were included in the stoichiometric core of the original network were chosen arbitrarily (within the constraint that the core be invertible). The potential advantages or disadvantages of such a choice are not addressed in this paper.

6 Conclusion

The results in this paper extend the toolkit for taking advantage of stoichiometric structure in the analysis of network behaviour. Complete factorization of the stoichiometry matrix (equation (7)) can be easily automated [13, 18] and is likely a worthwhile pre-processing step if steady-state network analysis demands a significant computational effort.

The model reduction approach that follows from this factorization may be useful in dealing with large-scale (e.g. genome-wide) networks. In particular, when targeting reversible internal reactions for elimination, this approach can facilitate elementary mode analysis of large input-output modules. Alternatively, the reduction procedure could be applied to achieve other goals, such as elimination of a pre-determined network submodule. Further efforts will be required to explore the potential and limitations of this approach and to develop a systematic methodology for its application.

Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council of Canada.

7 Appendix: Comparison of reduction methods

In [2], a method is described to eliminate redundancies in a network by writing all reaction rates in terms of a subset of measured reaction rates (taking advantage of metabolic balancing). That procedure can be described as follows.

Given a metabolic network, partition the species into extracellular metabolites (concentrations c) and intracellular metabolites (concentrations x). Writing the reaction rate vector as v, the dynamics of the extracellular concentrations can be described by

$\frac{d}{dt}\text{c}(t)={\text{K}}_1\text{v}(t)X_{\upsilon }(t)+\text{F}(\text{c}(t),t),$(14)

where K ₁ is the stoichiometry matrix for the extracellular species, X _v  (t) is the viable cell concentration at time t, and the function F describes mass transfer in and out of the reaction vessel. Denoting the stoichiometry matrix for the intracellular species as K ₂, a steady state assumption for these intracellular species leads to the balance equation

${\text{K}}_2\text{v}=0$

The assumption is then made that constraints are imposed on the steady-state reaction flux of the form

${\text{K}}_3\text{v}=\text{f}(\text{c})$(15)

where the function f is known, such that the matrix

$\left[\begin{array}{c}\hfill {\text{K}}_2\hfill \\ \hfill {\text{K}}_3\hfill \end{array}\right]$

is invertible. (This condition requires that (i) there are no dependent internal species (i.e. K ₂ has full row rank) and (ii) the number of rows of matrix K ₃ is equal to the difference between the total number of fluxes and the number of internal species.) The steady-state reaction fluxes are then given by

$\text{v}={\left[\begin{array}{c}\hfill {\text{K}}_2\hfill \\ \hfill {\text{K}}_3\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \text{I}\hfill \end{array}\right]\text{f}(\text{c})$

and the intracellular network can be eliminated by writing the dynamics for the extracellular species as

$\frac{d}{dt}\text{c}(t)=\text{Kf}(\text{c}(t))X_u(t)+\text{F}(\text{c}(t),t),$

where the reduced stoichiometry matrix K is given by:

$\text{K}={\text{K}}_1{\left[\begin{array}{c}\hfill {\text{K}}_2\hfill \\ \hfill {\text{K}}_3\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \text{I}\hfill \end{array}\right].$

In the simplest case, the constraints (15) are given as measurements of a set of reaction rates, so that f is constant and K ₃ has the form

${\text{K}}_3=[\begin{array}{cc}\hfill \text{I}\hfill & \hfill 0\hfill \end{array}]$

(where the reactions have been indexed so the measured rates appear as the first components in the vector v). In this case, the reduced stoichiometry matrix can be derived explicitly, as follows. Partitioning matrices K ₁ and K ₂ as

${\text{K}}_1=\left[\begin{array}{cc}\hfill {\text{K}}_{11}\hfill & \hfill {\text{K}}_{12}\hfill \end{array}\right]\text{and}{\text{K}}_2=[{\text{K}}_{21}{\text{K}}_{22}]$(16)

(where the dimensions of the submatrices align with the partitioning of K ₃), we can write the reduced stoichiometry matrix K as

$\begin{array}{l}{\text{K}}_1{\left[\begin{array}{c}\hfill {\text{K}}_2\hfill \\ \hfill {\text{K}}_3\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \text{I}\hfill \end{array}\right]\\ =[{\text{K}}_{11}{\text{K}}_{12}]{\left[\begin{array}{cc}\hfill {\text{K}}_{21}\hfill & \hfill {\text{K}}_{22}\hfill \\ \hfill \text{I}\hfill & \hfill 0\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \text{I}\hfill \end{array}\right]\\ =[{\text{K}}_{11}{\text{K}}_{12}]\left[\begin{array}{cc}\hfill 0\hfill & \hfill \text{I}\hfill \\ \hfill {\text{K}}_{\text{22}}^{\text{-1}}\hfill & \hfill {\text{K}}_{22}^{\text{-1}}{\text{K}}_{21}\hfill \end{array}\right]\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \text{I}\hfill \end{array}\right]\\ =[{\text{K}}_{11}{\text{K}}_{12}]\left[\begin{array}{c}\hfill \text{I}\hfill \\ \hfill \text{-}{\text{K}}_{\text{22}}^{\text{-1}}{\text{K}}_{21}\hfill \end{array}\right]\\ ={\text{K}}_{11}-{\text{K}}_{12}{\text{K}}_{22}^{-1}{\text{K}}_{21}.\end{array}$

This formula for the reduced stoichiometry matrix can be recovered using the reduction method described in this paper, as follows.

Replacing the transport processes F (from equation (14)) with a collection of inward-oriented single-species exchange reactions, we can formulate a stoichiometry matrix for the full network (both extracellular and intracellular species) as

$\text{N}=\left[\begin{array}{cc}\hfill {\text{K}}_1\hfill & \hfill \text{I}\hfill \\ \hfill {\text{K}}_2\hfill & \hfill 0\hfill \end{array}\right].$

Partitioning K ₁ and K ₂ as in equation (16) above, we have

$\text{N}=\left[\begin{array}{ccc}\hfill {\text{K}}_{11}\hfill & \hfill {\text{K}}_{12}\hfill & \hfill \text{I}\hfill \\ \hfill {\text{K}}_{21}\hfill & \hfill {\text{K}}_{22}\hfill & \hfill 0\hfill \end{array}\right].$

To apply a column factorization, we first solve for a nullspace matrix, W, which satisfies

$\text{NW}=\left[\begin{array}{ccc}\hfill {\text{K}}_{11}\hfill & \hfill {\text{K}}_{12}\hfill & \hfill \text{I}\hfill \\ \hfill {\text{K}}_{21}\hfill & \hfill {\text{K}}_{22}\hfill & \hfill 0\hfill \end{array}\right]\text{w}=0.$

Partitioning W, and assuming the desired form, we have

$\text{NW}=\left[\begin{array}{ccc}\hfill {\text{K}}_{11}\hfill & \hfill {\text{K}}_{12}\hfill & \hfill \text{I}\hfill \\ \hfill {\text{K}}_{21}\hfill & \hfill {\text{K}}_{22}\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill I\hfill \\ \hfill {\text{W}}_1\hfill \\ \hfill {\text{W}}_2\hfill \end{array}\right]=0.$

This gives

$\begin{array}{l}{\text{K}}_{11}+{\text{K}}_{12}{\text{W}}_1+{\text{W}}_2=0\\ {\text{K}}_{21}+{\text{K}}_{22}{\text{W}}_1=0.\end{array}$

Solving, we have

${\text{W}}_1=-{\text{K}}_{22}^{-1}{\text{K}}_{21}{\text{W}}_2=-{\text{K}}_{11}+{\text{K}}_{12}{\text{K}}_{22}^{-1}{\text{K}}_{21}.$

The nullspace matrix is then

$\text{W}=\left[\begin{array}{c}\hfill \text{I}\hfill \\ \hfill -{\text{K}}_{22}^{-1}{\text{K}}_{21}\hfill \\ \hfill -{\text{K}}_{11}+{\text{K}}_{12}{\text{K}}_{22}^{-1}{\text{K}}_{21}\hfill \end{array}\right],$

from which we can construct the column link matrix as

$P=\left[\begin{array}{ccc}\hfill {\text{K}}_{22}^{-1}{\text{K}}_{21}\hfill & \hfill \text{I}\hfill & \hfill 0\hfill \\ \hfill {\text{K}}_{11}-{\text{K}}_{12}{\text{K}}_{22}^{-1}{\text{K}}_{21}\hfill & \hfill 0\hfill & \hfill \text{I}\hfill \end{array}\right].$

Taking the stoichiometric core as

$N_R^C=\left[\begin{array}{cc}\hfill {\text{K}}_{12}\hfill & \hfill \text{I}\hfill \\ \hfill {\text{K}}_{22}\hfill & \hfill 0\hfill \end{array}\right],$

the stoichiometry matrix satisfies N = NRCPP.

Eliminating the internal reactions from the stoichiometric core results in a reduced core:

${\text{N}}_{\text{R}}^{\text{C}}=\left[\begin{array}{cc}\hfill 0\hfill & \hfill \text{I}\hfill \\ \hfill 0\hfill & \hfill \text{0}\hfill \end{array}\right].$

The reduced stoichiometry is then given by the product

${\text{N}}_{\text{R}}^{\text{C}}\text{P}=\left[\begin{array}{ccc}\hfill {\text{K}}_{\text{11}}-{\text{K}}_{\text{12}}{\text{K}}_{\text{22}}^{\text{-1}}{\text{K}}_{21}\hfill & \hfill 0\hfill & \hfill \text{I}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right].$

Removing the zero columns, we arrive at a stoichiometry that is identical to the reduction from [2] (along with the exchange reactions).

Comparing the two methods, we find they are complementary. The method of Hagg et al. has the advantage that complex constraints can be applied. The method described in this paper has the advantage that the user can choose how many reactions to eliminate, and specific reactions can be targeted for elimination.