A mixed-integer linear programming approach to the reduction of genome-scale metabolic networks.

BACKGROUND: Constraint-based analysis has become a widely used method to study metabolic networks. While some of the associated algorithms can be applied to genome-scale network reconstructions with several thousands of reactions, others are limited to small or medium-sized models. In 2015, Erdrich et al. introduced a method called NetworkReducer, which reduces large metabolic networks to smaller subnetworks, while preserving a set of biological requirements that can be specified by the user. Already in 2001, Burgard et al. developed a mixed-integer linear programming (MILP) approach for computing minimal reaction sets under a given growth requirement.
RESULTS: Here we present an MILP approach for computing minimum subnetworks with the given properties. The minimality (with respect to the number of active reactions) is not guaranteed by NetworkReducer, while the method by Burgard et al. does not allow specifying the different biological requirements. Our procedure is about 5-10 times faster than NetworkReducer and can enumerate all minimum subnetworks in case there exist several ones. This allows identifying common reactions that are present in all subnetworks, and reactions appearing in alternative pathways.
CONCLUSIONS: Applying complex analysis methods to genome-scale metabolic networks is often not possible in practice. Thus it may become necessary to reduce the size of the network while keeping important functionalities. We propose a MILP solution to this problem. Compared to previous work, our approach is more efficient and allows computing not only one, but even all minimum subnetworks satisfying the required properties.

Background

In computational systems biology, genome-scale metabolic network reconstructions are used to build in silico models of cellular metabolism [1]. To analyze these models, a large variety of constraint-based methods has been developed over the years [2].

Typically, the metabolic network is assumed to be in steady-state, i.e., the production and consumption of the internal metabolites has to be balanced. This leads to a flux space of the form . Here denotes the stoichiometric matrix, given a set of (internal) metabolites Met and a set of reactions Rxn. The vectors v∈C are called (feasible) flux vectors and can be interpreted as steady-state flux distributions of the metabolic network. The vectors define lower and upper bounds on the fluxes, where . By Irrev⊆Rxn we denote the set of irreversible reactions, which can carry flux in only one direction, i.e., v ≥0, for all i∈Irrev. For simplicity, we assume l ≥0, for all i∈Irrev. Reactions in Rev=Rxn∖Irrev are called reversible.

Some constraint-based analysis methods can be applied to genome-scale network reconstructions with several thousands of reactions. Others are limited to small or medium-sized models, like the computation of elementary flux modes [3] or minimal cut sets [4]. In such situations, a natural question is whether it is possible to reduce the given large network to a meaningful smaller one of practical size.

In 2015, Erdrich et al. [5] introduced a method called NetworkReducer, which reduces large metabolic networks to smaller subnetworks, while preserving relevant biological properties of interest. The algorithm in [5] is divided into two parts: network pruning and network compressing. In the compressing step, reactions belonging to the same enzyme subset [6] are lumped together. In the pruning step removable and non-removable reactions are identified such that the reduced network consisting of the non-removable reactions fulfills four requirements, which can be specified by the user:

Set of protected metabolites P Met: all metabolites in P Met must be retained in the reduced network.

Set of protected reactions P Rxns: all reactions in P Rxns must be retained in the reduced network.

Set of protected functionalities (or phenotypes) for the reduced network. We assume that any protected functionality can be described by a corresponding system of linear inequalities: D v≤d .

Minimum degrees of freedom: dof≥dof min. Here, the degrees of freedom dof correspond to the dimension of the null space of the stoichiometric matrix S, i.e., dof=|Rxn|−rank(S).

The overall goal of NetworkReducer is to find a smaller subnetwork such that all requirements (1) – (4) can be satisfied by a suitable flux vector. An example is given in Fig. 1.

Fig. 1

Solid arcs correspond to active reactions, dotted arcs to inactive reactions. In a, the flux vector satisfies the functionality of carrying flux through the biomass reaction while having oxygen uptake. In b, the functionality is carrying flux through the biomass reaction while there is no oxygen uptake. Combining the two flux vectors leads to the network in c, which contains 7 active reactions. A minimum subnetwork enabling both functionalities with only 6 reactions is given in (d). The corresponding binary variables for 1d would have the following values: a 1=1,a 2=1,a 3=1,a 4=1,a 5=1,a 6=0,a 7=0,a 8=1, where a corresponds to reaction r

The method of Erdrich et al. [5] searches for a suitable subnetwork by iterating over the reactions. In every iteration, the flux rate through one particular reaction is set to zero and a linear program (LP) is solved to check if the remaining reactions still form a feasible subnetwork. Feasibility means that there exists non-zero flux vectors satisfying the steady-state condition and the other requirements. To identify the reaction to be eliminated a flux variability analysis (FVA) [7] is done and a reaction with smallest overall flux range is selected. Thus in every iteration, an LP is solved and an FVA is performed. Each FVA involves solving up to 2n LPs, where n is the number of reactions.

An important aspect of the method in [5] is that it does not necessarily compute a minimum subnetwork (with respect to the number of active reactions), see Fig. 2 for an example. The method that we develop here will always find a feasible subnetwork with a minimum number of active reactions. A subnetwork satisfying the requirements (1) – (3) can be obtained by solving only one mixed-integer linear program (MILP). If this subnetwork does not fulfill the dof-requirement (4), we exclude this subnetwork and compute a new subnetwork by solving the MILP again. This method turns out to be much faster than the algorithm introduced in [5]. More importantly, we are guaranteed to obtain a minimum subnetwork regarding the number of active reactions, which is not the case for NetworkReducer. However, due to the minimality condition, our method cannot preserve flux variability in the same way as NetworkReducer does.

Fig. 2

If in the first step of the pruning procedure the flux through reaction 1 is set to zero, reaction 1 is removable and reactions 2 and 3 are non-removable. If in the first step reaction 2 or 3 is set to zero, both of them would be removable and reaction 1 would be non-removable. The resulting subnetwork is then smaller than the first one

A second related work is the FASTCORE algorithm of N. Vlassis et al. [8]. This method is also based on solving several LPs but without performing an FVA in between. Thus it is a very fast approach. However, the resulting subnetworks are not minimum and only protected reactions can be specified, but no protected metabolites, functionalities, or degrees of freedom.

An early approach for network reduction was introduced by Burgard et al. [9] already in 2001 and later improved in 2014 by Jonnalagadda and Srinivasan [10]. This method also allows computing minimum subnetworks using an MILP approach. However, only one functionality can be formulated and not several ones like in NetworkReducer. Other related work can be found in [11-18].

Altogether, our method can be seen as a network reduction algorithm that merges features from NetworkReducer and the method in [9], such that we can specify biological requirements like in [5] and compute all minimum subnetworks using an MILP, similar to [9].

The organization of this paper is as follows. In the Methods section we develop the underlying MILP methods. We start with the basic algorithm and then describe several improvements. In the Results and discussion section we compare our MILP approach with the existing methods NetworkReducer and FASTCORE. Furthermore, we apply it to a collection of genome-scale network reconstructions and discuss the results. The last section is Conclusion.

A software tool implementing the algorithms described in this paper is available at https://sourceforge.net/projects/minimalnetwork/.

Methods

Basic MILP to compute a minimum subnetwork

We always assume that our network is in steady-state, i.e., Sv=0, with bounds on the reaction rates l≤v≤u. Each functionality is described by a system of linear inequalities: D v≤d . For example, we may require that the biomass reaction has to carry at least 99% of its maximal rate: v Bio≥0.99· max(v Bio).

We will use binary variables a ∈{0,1} to indicate whether or not reaction i carries flux in the subnetwork. Thus we need the relationship a =0 if and only if v =0. For an irreversible reaction i∈Irrev, this can be achieved using constraints of the form

For reversible reactions, we use another binary variable and the constraints

This type of constraints is called a big M constraint, where M≫0 is a sufficiently large constant, e.g. some upper bound on the flux rates. With δ>0 we denote a threshold indicating above which flux rate a reaction is considered to be active. Practically δ will be chosen between 1e−06 and 1e−04.

To force protected irreversible reactions to carry flux, we use the constraints a =1 for all i∈P Irr=P Rxns∩Irrev. Enforcing flux through a protected reversible reaction can be realized in a similar way with the constraints for all i∈P Rev=P Rxns∩Rev.

For any protected metabolite m∈P Met, let Rxn be the set of reactions involving m. By Rev we denote the set of reversible reactions in Rxn. If Rxn contains at least one protected reaction r, metabolite m will be protected by reaction r. However, if Rxn∩P Rxn=∅, further constraints are needed to protect m:

where .

In [5], an additional requirement is to specify a minimum number of active reactions. Here we do not include this restriction for the following reasons. First, we will search for the minimum number of active reactions such that all the other requirements are fulfilled. Second, in [5] the minimum number of active reactions is always set to 1. Since there exist reactions which are forced to carry flux, this constraint is redundant.

To find a subnetwork which contains the minimum number of active reactions, we minimize over the sum of the binary variables a , which indicate whether a reaction carries flux. The resulting MILP is the following:

Conflicting functionalities

In the case study considered in [5], the resulting subnetwork should keep two desired functionalities: under both aerobic and anaerobic conditions at least 99.9% of the maximal growth rate should be maintained. These two requirements cannot be realized with the same flux vector v because they imply two opposite states of the reaction that transports O 2 into the network. We would need a vector v with and at the same time, which is not possible.

MinNW-0 computes one feasible flux vector v of the network. But, to get a subnetwork which fulfills the two functionalities we need one flux vector which fulfills the aerobic condition and another one for the anaerobic condition, see Fig. 1. To realize this with a single MILP we have to modify MinNW-0. First, we search for a flux vector v 0 which contains the protected metabolites and protected reactions. Additionally, for each functionality we search for a flux vector v satisfying D v ≤d and corresponding binary variables. For example, in Fig. 1, we would have a 1=1 in case 1a) and a 1=0 in case 1b). Due to (1) and (2), this would imply a 1=1 and a 1=0 at the same time, which is not possible. Thus we have to use different binary variables for . With this, the Eqs. (1) and (2) become

Using the new variables , we reformulate the constraints regarding the protected reactions: , for all i∈P Irr and , for all i∈P Rev. Finally, the constraints regarding the protected metabolites become

To obtain a minimum subnetwork, we have to minimize the total number of active reactions. Thus, we need binary variables a with the property For irreversible reactions, this can be encoded by the constraints

and for reversible reactions we get

The resulting MILP is the following: minNW computes a subnetwork with a minimum number of active reactions while satisfying all the requirements.

Example for minNW

The network in Fig. 1 a fulfills the functionality regarding the aerobic condition, while the network in Fig. 1 b fulfills the anaerobic condition. The combination of the minimum subnetworks corresponding to each functionality does not lead to a minimum subnetwork for both, see Fig. 1. The minimum subnetwork for this example is given in Fig. 1 d.

Computing all minimum subnetworks

There are scenarios where we have to compute more than one subnetwork. For instance, consider the case where the minimum dof (requirement (4)) is larger than 1. If the subnetwork computed with minNW does not have the required dof, we have to compute a different subnetwork. Furthermore, the computed minimum subnetwork need not be unique. Thus there may exist different subnetworks which all fulfill the requirements and have the same number of active reactions. So we may be interested in finding all minimum subnetworks. To compute different subnetworks we can use the MILP minNW in an iterative way. Whenever a minimum subnetwork is found, we formulate a constraint which excludes this subnetwork as a feasible solution and solve the (extended) MILP again. For that purpose we formulate the following constraints:

where if reaction i carries flux in the subnetwork which was computed in the k-th step, otherwise . Thus (10) guarantees that at least one inactive reaction will become active, or at least one active reaction will become inactive in the new solution.

Solving minNW iteratively and adding the constraints (10) in each step, we are now able to enumerate all minimum subnetworks.

Reducing the number of binary variables

To further improve efficiency, we will make use of flux coupling information [19-22]. We first recall some basic definitions from flux coupling analysis (FCA).

A reaction r∈Rxn is called blocked if v =0 for all . In a pre-processing step, blocked reactions will be removed from the network, which is also done in [5]. Thus we assume from now on that the network contains only unblocked reactions.

Given two unblocked reactions r,s∈Rxn, we say r is partially coupled to s, and write r⇔s, if v =0⇔v =0, for all v∈C 0. The relation r⇔s is reflexive, transitive and symmetric and therefore defines an equivalence relation on Rxn. This means that the set of reactions Rxn can be partitioned into equivalence classes [r]={s∈Rxn∣r⇔s }. It follows , where denotes the set of all equivalence classes. An equivalence class can be represented by any of its elements. We say that r is a representative of [r] or that [r] is the coupling class of r. Note that [r]=[s] iff r⇔s. Biologically, coupling classes can be interpreted as subsets of reactions that are always active together at steady-state, similarly to the notion of enzyme subsets in [6].

The main advantage of introducing coupling classes is that, if one reaction in a class is not carrying flux, no other reaction in the class does, and vice versa. Therefore, in every approach where binary variables are used to indicate if a reaction appears or not, it suffices to consider one reaction from every coupling class instead of considering all of them. Depending on the number of reactions and associated coupling classes, this may significantly reduce the number of required variables.

Based on the equivalence relation r⇔s, we now use binary variables corresponding to the coupling classes [r] instead of having binary variables for each individual reaction. Thus we can rewrite the algorithm minNW in the following way:

Here, |[r]| denotes the cardinality of the coupling class [r]. Thus, we compute the smallest subnetwork with respect to the number of active reactions and not with respect to to the number of active representatives. denotes the representatives of the irreversible reactions, and those of the reversible reactions. Similarly, resp. is the set of representatives of protected irreversible resp. protected reversible reactions. With we denote the representatives which include a protected metabolite.

To exclude previously enumerated subnetworks the constraints (10) can be adapted in the following way:

Using representatives we need only instead of |Rxn| binary variables. For many genome-wide networks, this reduces the number of 0-1 variables by about 1/2, see the examples in Table 1.

Table 1

Number of representatives for different genome-wide metabolic networks (computed with F2FC [20])

Model	Reactions	Unblocked	Coupling classes
Mus musculus	3726	2436	1489
E. coli iJO1366	2583	2369	1399
S. Typhimurium LT2	2545	1620	1047
S. boydii CDC 3083-94	2592	1546	1016
K. pneumoniae MGH 78578	2262	1223	804
Y. pestis CO92	1961	1065	639
S. cerevisiae S288c	1577	885	558
G. metallireducens GS-15	1285	845	330
M. tuberculosis iNJ661	1025	800	412
B. subtilis 168	1250	658	342
P. putida KT2440	1056	652	282
C. ljungdahlii DSM 13528	785	526	215
H. pylori iIT341	554	501	209
M. barkeri str. Fusaro	690	484	174
S. aureus iSB619	743	465	224
T. maritima MSB8	652	385	148

Results and discussion

We implemented our MILPs in MATLAB and used CPLEX [23] as a solver like in [5]. For NetworkReducer resp. FASTCORE we used the implementation provided by the authors of [5] resp. [8]. All computations were done on a desktop machine with two processors Intel(R) Core(TM) i5-2400S, CPU 2.50GHZ, each 1 thread. For algorithm minNW rep, we computed the coupling classes for partially coupled reactions using the software F2C2 [20].

Indicator variables

We implemented two versions of our algorithms. In one version we used the big M constraints from the original MILP formulation in the Methods section. We observed that the solutions are highly dependent on the given tolerances in the MILP solver. To increase numerical stability, we implemented a second version using indicator variables and some other features of CPLEX [23]. The use of indicator variables is straightforward. For example, the big M constraint δ a≤v≤M a is replaced by a=0⇒v=0,a=1⇒v≥δ, where a∈{0,1} is the indicator variable. MILP solvers using indicator variables handle them in two different ways. They may reformulate the given indicator constraints into big M constraints or branch on the indicator variables. CPLEX chooses one of these two methods depending on M. If M is small, it will formulate big M constraints, otherwise it will use branching. For the results and the running time we only applied the version where indicator variables were used, due to numerical instability of the big M formulation. While indicator variables drastically increase the running time, we still outperform the algorithm in [5].

Comparison with NetworkReducer

In a first experiment, we ran our implementations on the two metabolic network reconstructions and functionalities considered in [5]. Table 2 shows the running time for calculating a subnetwork with the desired properties.

Table 2

Time (in seconds) needed to compute a subnetwork with given requirements resp. constraints

Algorithm	Synechocystis sp. PCC 6803		E. coli iAF1260
	Time	reactions	Time	reactions
NetworkReducer	324	462	21987	455
minNW	31	453	4074	416

NetworkReducer: The algorithm introduced in [5]. minNW: The MILP introduced here, using indicator variables

For Synechocystis sp. PCC 6803, the subnetwork computed by NetworkReducer [5] contains 462 reactions and thus 9 reactions more than the minimum subnetwork with 453 reactions obtained by our method. The two subnetworks have 413 reactions in common. 49 reactions in the larger subnetwork cannot be found in the minimum subnetwork, while 40 reactions in the minimum subnetwork do not appear in the larger one.

Regarding E. coli iAF1260 we get similar results. The subnetwork computed by NetworkReducer contains 39 reactions more than the minimum subnetwork obtained by our method. Both networks have 424 reactions in common. There are 51 reactions that can only be found in the subnetwork computed with NetworkReducer, while there are 12 reactions which appear only in the minimum subnetwork.

Comparison with FASTCORE

FASTCORE [8] is a heuristic algorithm which is much faster than our method. However, the computed subnetworks are not minimum as can be seen from Table 3. The subnetwork computed with our method is not contained in the subnetwork computed with FASTCORE.

Table 3

rxns: number of unblocked reactions in the original network

		rxns	time	rxns	time
Models	rxns	FASTCORE	FASTCORE	minNW	minNW
M. tuberculosis iNJ661	1025	134	0.12	62	8727
H. pylori	501	319	0.6	306	123
26695

rxns FASTCORE: number of the reactions in the subnetwork computed with FASTCORE. time FASTCORE: running time in seconds of FASTCORE. rxns minNW: number of the reactions in the subnetwork computed with minNW. time minNW: running time in seconds for the algorithm minNW using indicator variables

For H. pylori 26695 there are 22 reactions that appear only in the FASTCORE subnetwork and 9 reactions which can be found only in the minimum subnetwork. Similarly, for M. tuberculosis iNJ661, there are 78 reactions that appear only in the FASTCORE subnetwork and 6 reactions which can be found only in our subnetwork. The names of the reactions for both examples are given in the (Additional file 1).

Network reduction for genome-scale metabolic networks

As a proof of concept we applied our methods to compute minimum subnetworks for 16 metabolic network reconstructions taken from BiGG Models [24] under different scenarios. For each type of organism in BiGG we considered one model (except for human recon because there is no biomass reaction). An overview of the results is given in Table 4. In some cases we had only one minimum subnetwork, while for some models and scenarios we found different ones. For example, in the case of H. pylori 26695, we get 16 distinct minimum subnetworks, which will be discussed in the section Case study: .

Table 4

Computational results using indicator variables

			ess	rxns	mets	reps	time	time
Models	rxns	mets	rxns	in SNW	in SNW	in SNW	minNW	minNW rep	SNWs
Mus	2436	1665	247	351	351	241	49085	2949	1
musculus
E. coli	2369	1159	363	562	601	262	704	587	1
iJO1366
S. Typhimurium	1620	1098	305	458	455	277	1565	1507	1
LT2
S. boydii CDC	1546	1019	441	445	450	209	15898	463	1
3083-94
K. pneumoniae	1223	830	203	338	340	188	299	194	1
MGH 78578
Y. pestis	1065	761	279	339	339	171	7544	5970	1
CO92
S. cerevisiae	885	639	262	290	289	195	1225	720	2
S288c
G. metallireducens	845	710	544	557	567	153	257	49	1
GS-15
M. tuberculosis	800	580	314	427	425	168	4065	811	1
iNJ661
B. subtilis	658	500	270	296	300	134	16854	10027	1
168
P. putida	652	539	300	344	348	116	3784	827	7
KT2440
C. ljungdahlii	526	448	369	383	389	118	26	7.6	44
DSM 13528
H. pylori	501	381	265	321	323	89	9.8	9.8	16
26695
M. barkeri str.	484	417	289	364	369	90	25.38	24.3	20
Fusaro
S. aureus	465	387	71	122	127	75	28	27	1
iSB619
T. maritima	385	331	267	282	280	87	14	5.59	28
MSB8

rxns: number of unblocked reactions in the original network. mets: number of metabolites in the original network after removing dead-end-metabolites. ess rxns: number of essential reactions in the original network. rxns in SNW: number of reactions in the subnetwork. mets in SNW: number of metabolites in the subnetwork. reps in SNW: number of representatives remaining in the subnetwork. time minNW: running time in seconds for the algorithm minNW. time minNW rep: running time in seconds for the algorithm minNW rep. SNWs: number of minimum subnetworks which exist and fulfill all the requirements. For detecting the running time, only one subnetwork was computed

Following [4], we call a reaction essential if after removing this reaction it is no longer possible to achieve at least p % of the maximal biomass production rate. Like in [4], we choose p=20. A minimum subnetwork where it is possible to achieve a maximal biomass rate constitutes a subnetwork where all essential reactions must be active and so all essential reactions have to be included in the subnetwork. We will present the number of essential reactions for the different models to give an idea how many reactions are additionally needed to have a functional minimum subnetwork including all essential reactions.

The scenarios for the different networks and some conclusions are given next, full details can be found in the (Additional file 1). The bounds on the flux rates are those from BiGG Models.

For the networks Mus musculus, E. coli iJO1366, S. Typhimurium LT2, S. boydii CDC 3083-94, and K. pneumoniae MGH 78578 the requirements are that at least 99.9% of the maximal biomass rates for the aerobic and anaerobic case can be realized by the subnetwork. For Y. pestis CO92 the requirements are that at least 99.9% of the maximal growth rate with glycine uptake and the maximal growth rate without glycine uptake can be realized by the subnetwork. For S. cerevisiae S288c the maximal biomass rate with and without ethanol exchange has to be realized by the reduced subnetwork. For G. metallireducens GS-15, C. ljungdahlii DSM 13528, and T. maritima MSB8 the maximal biomass rate with H2O uptake and without H2O exchange has to be realized by the reduced subnetwork. For M. tuberculosis iNJ661 one requirement is that at least 99.9% of the maximal growth rate can be achieved. Additionally we defined 36 protected reactions. For B. subtilis 168 the requirements are that at least 99.9% of the maximal growth rate with hydrogen uptake and the maximal growth rate without hydrogen uptake can be realized by the subnetwork. For P. putida KT2440 one requirement is that at least 99.9% of the maximal growth rate can be achieved. Additionally we defined protected reactions to keep the TCA cycle. For H. pylori 26695 one requirement is that at least 99.9% of the maximal growth rate can be achieved. Additionally we defined 28 protected reactions. A detailed discussion of this test case will be given in the next subsection. For M. barkeri str. Fusaro the requirements are that at least 99.9% of the maximal growth rate with ammonia uptake and the maximal growth rate without ammonia uptake can be realized by the subnetwork. For S. aureus N315 at least 99.9% of the maximal biomass rate with glucose uptake and without glucose uptake has to be realized by the subnetwork.

Case study: Helicobacter pylori 26695

In this section we discuss the results for computing several minimum subnetworks for the metabolic network H. pylori 26695 using indicator variables. The requirements are the following:

There are 28 protected reactions.

The maximal biomass yield is 20.2606, and the subnetworks should be able to produce at least 99.9% of this yield.

In total we computed 16 subnetworks each containing 321 reactions, which is the minimum number needed to fulfill the requirements. The time needed to compute all these minimum subnetworks was 127 seconds with minNW and 33 seconds with minNW rep. Altogether the 16 minimum subnetworks use 329 different reactions, which can be found in the (Additional file 1). 311 reactions are present in every subnetwork, among them all the 265 essential reactions of H. pylori. Only 18 reactions are not present in every subnetwork: CCP, G3PD1, D-Amino acid dehydrogenase, FUMt3, Glycerol-3-phosphate dehydrogenase (NADP), SUCFUMt, L-alanine dehydrogenase, Anthranilate synthase, Formate-tetrahydrofolate ligase, D-Alanine exchange, D-alanine transport via proton symport, L alanine reversible transport via proton symport, L-Alanine exchange, ANS2, GAR transformylase-T, NO3t2, NO3t3, Catalase. Figure 3 shows the distributions of these reactions in the 16 subnetworks.

Fig. 3

The two illustrations show the distribution of the reactions which are not present in all subnetworks. In Fig. 3 a each reaction (x-axis) has a bar. The bar indicates in how many subnetworks the reaction can be found. For example, reaction CCP can be found in every subnetwork except 1 (there are in total 16 subnetworks) and reaction CAT can be found in only one subnetwork. Fig. 3 b illustrates where the reactions are found. Again the x-axis corresponds to the reactions. Thus a dot at (1,CCP) means that CCP appears in subnetwork 1. CCP can be found in every subnetwork except in the second one, whereas CAT can be found only in the second one

Additional insight can be obtained by analyzing co-occurrence patterns of the 18 non-essential reactions. Some of these reactions are mutually exclusive regarding the minimum subnetworks. For example, all subnetworks that contain reaction CCP do not contain CAT and vice versa. The same holds for the pair FTHFLi and GART, and the pair ANS and ANS2. Regarding the functionalities of these reaction pairs, one can easily check that the two reactions in each pair do basically the same. Therefore, it is sufficient if only one of them is present. In opposite to this, we can see that DALAt2r and EX_ala__D(e) form a cycle since they always appear in the same subnetworks. The same holds for ALAt2r and EX_ala__L(e). Both cycles also seem to be mutually exclusive, thus only one of them is present in the subnetworks. Similar observations can be made for the cycle formed by NO3t2 and NO3t3, which is mutually exclusive to the cycle formed by SUCFUMt and FUMt3.

One may ask whether the reactions that never appear together in the same subnetwork are also mutually exclusive regarding elementary flux modes (EFMs), i.e., whether or not there exists an EFM involving both reactions [25]. While this holds for the reaction pair FTHFLi and GART and the pair CCP and CAT, it is not true for the other reactions.

Conclusion

We developed an MILP approach to compute for a given large metabolic network one or more minimum subnetworks preserving biological requirements that can be specified by the user. Compared to previous work [5, 8, 9], our method guarantees minimality of the subnetwork regarding the number of active reactions while preserving all the given requirements. In case there exist several minimum solutions, we are able to enumerate all of them. This may give additional insight how the network is functioning and which reactions are really needed to satisfy the requirements. We applied our algorithms to several genome-scale metabolic networks and we always found all the minimum subnetworks in reasonable time.

Once these subnetworks have been computed, further analysis becomes possible by using methods that are not applicable to the original network. For example, one may compute elementary flux modes and minimal cut sets. In addition, one can take a closer look to the reactions involved in one or all minimum subnetworks in order to get a better understanding of their role in the network.