You are viewing a javascript disabled version of the site. Please enable Javascript for this site to function properly.
Go to headerGo to navigationGo to searchGo to contentsGo to footer
In content section. Select this link to jump to navigation

Exploiting stoichiometric redundancies for computational efficiency and network reduction

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.

1Introduction

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].

2Rank 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

(1)
ddts(t)=NV(s(t),)p).

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.1Deficiencies 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 s=[siT,sdT]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

(2)
L=IrL0.

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=[L0 Inr],

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

ddts(t)=ddtsi(t)sd(t) = LNRv(s(t)p).

Equation (2) then gives

ddtsd(t)=ddtL0si(t)      for all t0.

Integrating gives an explicit dependence among the species concentrations: sd(t)=L0si(t)+T˜ for all time, where T˜=sd(0)-L0si(0). Finally, concatenating T˜ with the r-vector of zeros 0 r , define T=[0rT,T˜T]T so that

(3)
s(t) = Lsi(t)+T.

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

(4)
ddtsi(t)=NRv(Lsi(t)+T,p).

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

2.2Deficiencies 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 v=[viT,vdT]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=[P0Ir].

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

K=ImrP0,

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

(5)
ddts(t) = Nv(s(t),p)          =NCPv(s(t),p)          = NC[P0 Ir]v(s(t),p)

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 J=[JiT,JdT]T corresponding to the partitioning of the reaction rates described above. The dependence is given explicitly by

(6)
J=KJi=ImrP0Ji.

(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.3Complete 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:

(7)
ddts(t)=LNRCPv(s, p)          =IrL0 NRC[P0 Ir]v(s, p)

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

(8)
ddtsi(t)=NRCPv(Lsi+T, p)

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

(9)
0=Pv(Lsi(t)+T, p)

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.

3Computational efficiencies

3.1Metabolic 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

(10)
dsdp=LNRvsL1 NRvp

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

0=Pvs L dsidp+Pvp.

When the matrix PvsL is presumed invertible (by the standard assumption of asymptotic stability), one arrives at the (unscaled) response coefficients in the form

(11)
dsdp=LPvsL1Pvp.

The (unscaled) concentration control coefficients are then

Cs=LPvsL1P.

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 NR. 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.2Metabolic 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 n and N b and the steady state balance equation (from equation (1)) can be written as

(12)
Nnvn=Nbvb.

One then attempts to solve for rtn . In the simplest case N n is invertible and rtn is found exactly. More typically, one can make use of the pseudo-inverse of N 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=vinvibvdnvdb.

The steady state equation (9) is simply

Pv=[P0 Ir] v=0.

Partitioning P0=[P0nP0b] allows us to write thisbalance equation as

P0nP0bI00I  vinvibvdnvdb=0.

where the blocks in the identity matrix have dimensions corresponding to vdn and vdb. Isolating the unknown rates leads to a restatement of equation (12):

P0nI0 vinvdn=P0b0I  vibvib.

One is now faced with inverting [P0nI0] as opposed to the original task of inverting N n . Although these matrices are of the same size, the former will typically be easier to deal with. In particular, the inverse of [P0nI0] (if it exists) can be expressed directly in terms of the inverse of a submatrix of P0n, since in general

AIB01=0B1IAB1

if B is invertible. Furthermore, the same submatrix can be used to simplify computation of the pseudo-inverse of [P0nI0] through the use of the Schur complement [15]. The pseudo-inverse of

AIB0

is

AIB0T AIB01 AIB0

where

([AIB0]T[AIB0])-1=[(BTB)-1-(BTB)-1AT-A(BTB)-1I+A(BTB)-1AT].

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

4Steady-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=110101011

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

L=100111   and  P=110101

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

NRC=1101  gives  N=LNRCP=211101110

(Figure 1B) while

NRC=1021  gives  N=LNRCP=110121011

(Figure 1C).

The members of this equivalence class are networks composed of three species (S 1, S 2, S 3) involved in three reactions (rates v 1, v 2, v 3) for which (i) a conservation [S 1] - [S 2] + [S 3] = T holds for all time, and (ii) steady-state occurs when v 1 = v 2 = v 3. 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

(13)
N=110010100001110001000010110000100001000000000100010001010010

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

K=1000010000100001101010002110011110011000

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 7 = 2, v 9 = 2 and v 10 = 1 (units arbitrary). It can then be verified that v 1 = 1, v 4 = 1, v 6 = 1, v 8 = 1, while v 2, v 3 and v 5 are constrained only by the two conditions v 2 = v 3 = 1 - v 5.

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

NRC=101000000100010000100000100001010010

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.1Relaxation 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:

NRC=001100011000010000000000000011000000

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:

υ6=υ1   υ5=υ1υ3   υ7=2υ1+υ2υ3υ10=υ1   υ9=υ1+υ1   υ8=υ2+υ3+υ4.

Figure 3 makes it clear that most of these constraints hold in the alternative network. The exceptions are that v 5 is absent, and the two constraints on v 9 and v 10 have been merged into the single requirement that v 9 = v 10 + v 4.

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.

5Steady-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:

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

  • 2. 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.

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

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

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

  • 6. 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.1Exchange-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:

  • 1. Reactions targeted for elimination are eliminated.

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

  • 3. 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.2Flux 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:

υ50υ1υ30υ60υ10

The constraint v 6 ≥ 0 corresponds directly to an irreversibility constraint on v 1, which is retained in the reduced network. However, the constraint on v 5 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:

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

  • 2 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.

  • 3 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:

EFM1:  υ7υ9            EFM2:    υ8υ9EFM4:  2υ7υ9&υ10  EFM2:    2υ8υ9&υ10.

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.3Example: 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.

6Conclusion

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.

7Appendix: 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

(14)
ddtc(t)=K1v(t)Xυ(t)+F(c(t), t),

where K 1 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 2, a steady state assumption for these intracellular species leads to the balance equation

K2v=0

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

(15)
K3v=f(c)

where the function f is known, such that the matrix

K2K3

is invertible. (This condition requires that (i) there are no dependent internal species (i.e. K 2 has full row rank) and (ii) the number of rows of matrix K 3 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

v=K2K31  0I f(c)

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

ddtc(t)=Kf(c(t))Xu(t)+F(c(t), t),

where the reduced stoichiometry matrix K is given by:

K=K1K2K31 0I.

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

K3=[I0]

(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 1 and K 2 as

(16)
K1=K11K12   and  K2=[K21  K22]

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

K1K2K310I=[K11 K12]K21K22I01 0I=[K11 K12]0IK22-1K22-1K21 0I=[K11  K12] I-K22-1 K21=K11K12K221K21.

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

N=K1IK20.

Partitioning K 1 and K 2 as in equation (16) above, we have

N=K11K12IK21K220.

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

NW=K11K12IK21K220  w=0.

Partitioning W, and assuming the desired form, we have

NW=K11K12IK21K220  IW1W2=0.

This gives

K11+K12W1+W2=0    K21+K22W1=0.

Solving, we have

W1=K221K21    W2=K11+K12K221K21.

The nullspace matrix is then

W=IK221K21K11+K12K221K21,

from which we can construct the column link matrix as

P=K221K21I0K11K12K221K210I  .

Taking the stoichiometric core as

NRC=K12IK220,

the stoichiometry matrix satisfies N = NRCPP.

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

NRC=0I00.

The reduced stoichiometry is then given by the product

NRCP=K11K12K22-1K210I000.

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.

References

1 

Famili I, Palsson BØ2003Systemic metabolic reactions are obtained by singular value decomposition of genome-scale stoichiometric matricesJournal of Theoretical Biology2248796

2 

Haag JE, Vande Wouver A, Bogaerts P2005Systematic procedure for the reduction of complex biological reaction pathways and the generation of macroscopic equivalentsChemical Engineering Science60459465

3 

Heinrich R, Schuster S1996Chapman & HallThe Regulation of Cellular SystemsNew York

4 

Ingalls BP2009Stoichiometric redundancy: Computational efficiencies and steady-state-consistent network reductions, ColoradoProceedings of the Conference on Foundations of Systems Biology in Engineering (FOSBE) Denver

5 

Kaleta C, de Figueiredo LF, Schuster S2009Can the whole be less than the sum of its parts? Pathway analysis in genome-scale metabolic networks using elementary flux patternsGenome Research1918721883

6 

Klamt S, Stelling J Systems Modeling in Cellular Biology 2006 Stoichiometric and constraint-based modeling Szallasi Z, Stelling J, Periwal V MIT Press Cambridge Massachusetts

7 

Klamt S, Schuster S, Gilles ED2002Calculability analysis in underdetermined metabolic networks illustrated by a model of the central metabolism in purple nonsulfer bacteriaBiotechnology and Bioengineering77734751

8 

Marashi S-A, David L, Bockmayr A2012Analysis of metabolic subnetworks by flux cone projectionAlgorithms for Molecular Biology717

9 

Palsson BØ2006Systems Biology: Properties of reconstructed networksCambridge University PressCambridge Massachusetts

10 

Provost A, Bastin G2004Dynamic metabolic modelling under the balanced growth conditionJournal of Process Control14717728

11 

Reder C1988Metabolic control theory: A structural approachJournal of Theoretical Biology135175201

12 

Reed JL, Vo TD, Schilling CH, Palsson BØ2003An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR)Genome Biology4R54.1R54.12

13 

Sauro HM, Ingalls BP2004Conservation analysis in biochemical networks: Computational issues for software writersBiophysical Chemistry109115

14 

Schuster S, Dandekar T, Fell DA1999Detection of elementary flux modes in biochemical networks: A promising tool for pathway analysis and metabolic engineeringTrends in Biotechnology175360

15 

Serre D2002Matrices: Theory and ApplicationsNew YorkSpringer-Verlag

16 

Stephanopoulos GN, Aristidou AA, Nielsen J1998Metabolic Engineering: Principles and MethodologiesAcademic PressSan Diego

17 

Terzer M, Stelling J2008Large scale computation of elementary flux modes with bit pattern treesBioinformatics10.1093/bioinformatics/btn401

18 

Vallabhajosyula RR, Chickarmane V, Sauro HM2006Conservation Analysis of large biochemical networksBioinformatics22346353

19 

Wagner C2004Nullspace approach to determine the elementary modes of chemical reaction systemsJournal of Physical Chemistry B10824252431

20 

Westerhoff HV, Hofmeyr J-H.S, Kholodenko BN1994Getting to the inside of cells using metabolic control analysisBiophysical Chemistry50273283

Figures and Tables

Fig.1

Steady-state equivalent networks. All three networks have identical conservations and steady-state flux distributions. They are members of a single equivalence class as described in the text. (Multiple-headed arrows are used to indicate production of multiple copies of molecules from a single reaction event.

Steady-state equivalent networks. All three networks have identical conservations and steady-state flux distributions. They are members of a single equivalence class as described in the text. (Multiple-headed arrows are used to indicate production of multiple copies of molecules from a single reaction event.
Fig.2

Metabolic balancing example. A. Original network.B. Steady-state-equivalent network. In both cases, with flux measurements of v 7 = 2, v 9 = 2 and v 10 = 1 (units arbitrary), it follows that v 1 = 1, v 4 = 1, v 6 = 1, v 8 = 1, while v 2, v 3 and v 5 cannot be determined (beyond the constraints that v 2 = v 3 = 1 - v 5)

Metabolic balancing example. A. Original network.B. Steady-state-equivalent network. In both cases, with flux measurements of v
7 = 2, v
9 = 2 and v
10 = 1 (units arbitrary), it follows that v
1 = 1, v
4 = 1, v
6 = 1, v
8 = 1, while v
2, v
3 and v
5 cannot be determined (beyond the constraints that v
2 = v
3 = 1 - v
5)
Fig.3

Reduced version of the network in Fig. 2A. The double-headed arrow indicates that the reaction with rate v 1 produces two molecules of species A.

Reduced version of the network in Fig. 2A. The double-headed arrow indicates that the reaction with rate v
1 produces two molecules of species A.
Fig.4

Reductions of the network in Fig. 2A. A. Reaction 5eliminated. B. Reaction 6 eliminated.

Reductions of the network in Fig. 2A. A. Reaction 5eliminated. B. Reaction 6 eliminated.
Fig.5

Minimal exchange-equivalent version of the network in Fig. 2A.

Minimal exchange-equivalent version of the network in Fig. 2A.
Fig.6

Elementary modes of the network in Fig. 2A (original) and the reduced version in Fig. 5. Reaction rates are either zero (thin arrows), one (thick arrows), two (double-thick arrows) or negative one (reaction 5 in elementary modes 3 and 5).

Elementary modes of the network in Fig. 2A (original) and the reduced version in Fig. 5. Reaction rates are either zero (thin arrows), one (thick arrows), two (double-thick arrows) or negative one (reaction 5 in elementary modes 3 and 5).