Skip to main content

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

Abstract

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 \(C = \{ v\in \mathbb {R}^{\text {Rxn}} \mid Sv = 0, \; l \leq v \leq u\}\). Here \(S\in \mathbb {R}^{\text {Met} \times \text {Rxn} }\) denotes the stoichiometric matrix, given a set of (internal) metabolites Met and a set of reactions Rxn. The vectors vC are called (feasible) flux vectors and can be interpreted as steady-state flux distributions of the metabolic network. The vectors \(l,u \in \mathbb {R}^{\text {Rxn}}_{\pm \infty }\) define lower and upper bounds on the fluxes, where \(\mathbb {R}_{\pm \infty } := \mathbb {R}\cup \{\pm \infty \}\). By IrrevRxn we denote the set of irreversible reactions, which can carry flux in only one direction, i.e., v i ≥0, for all iIrrev. For simplicity, we assume l i ≥0, for all iIrrev. Reactions in Rev=RxnIrrev 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:

  1. 1.

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

  2. 2.

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

  3. 3.

    Set \(\mathcal {F}\) of protected functionalities (or phenotypes) for the reduced network. We assume that any protected functionality \(f \in \mathcal {F}\) can be described by a corresponding system of linear inequalities: D f vd f .

  4. 4.

    Minimum degrees of freedom: dofdof 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
figure1

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 i corresponds to reaction r i

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
figure2

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

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 lvu. Each functionality \(f\in \mathcal {F}\) is described by a system of linear inequalities: D f vd f . 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 i {0,1} to indicate whether or not reaction i carries flux in the subnetwork. Thus we need the relationship a i =0 if and only if v i =0. For an irreversible reaction iIrrev, this can be achieved using constraints of the form

$$\begin{array}{*{20}l} & \delta \, a_{i} \leq v_{i} \leq Ma_{i}. \end{array} $$
(1)

For reversible reactions, we use another binary variable \(\bar {a}_{i}\) and the constraints

$$\begin{array}{*{20}l} &\delta \, a_{i} - M \, \bar{a}_{i} \leq v_{i} \leq M \, a_{i} - \delta \, \bar{a}_{i}, \quad a_{i} + \bar{a}_{i} \leq 1. \end{array} $$
(2)

This type of constraints is called a big M constraint, where M0 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 i =1 for all iP Irr=P Rxns∩Irrev. Enforcing flux through a protected reversible reaction can be realized in a similar way with the constraints \( a_{i} + \bar {a}_{i} = 1,\) for all iP Rev=P Rxns∩Rev.

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

$$ \sum_{i \in \text{Rxn}_{m}} a_{i} + \sum_{i \in \text{Rev}_{m}} \bar{a}_{i} \geq 1, \quad \forall m \in {p^{Met}_{0}}, $$
(3)

where \(p^{Met}_{0} = \left \{ m \in P^{\text {Met}} \mid \text {Rxn}_{m} \cap P^{\text {Rxns}} = \emptyset \right \}\).

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 i , which indicate whether a reaction carries flux. The resulting MILP is the following:

$$\begin{aligned} & \mathbf{(MinNW-0)}\ \min \sum_{i \in\text{Rxn}} a_{i} + \sum_{k \in\text{Rev}} \bar{a}_{k} \\ & S v = 0, \ l \leq v \leq u & & \\ & D_{f} v \leq d_{f} & \forall f\in \mathcal{F} & \\ & \delta \, a_{i} \leq v_{i} \leq Ma_{i} & \forall i \in\text{Irrev} & \\ &\delta \, a_{i} - M \, \bar{a}_{i} \leq v_{i} \leq M \, a_{i} - \delta \, \bar{a}_{i} & \forall i\in\text{Rev}& \\ & a_{i} + \bar{a}_{i} \leq 1 & \forall i\in\text{Rev} & \\ & a_{i} = 1, \enspace a_{k} + \bar{a}_{k} = 1 & \forall i \in P^{\text{Irr}}, \forall k \in P^{\text{Rev}} & \\ & \sum_{i \in \text{Rxn}_{m}} a_{i} + \sum_{i \in \text{Rev}_{m}} \bar{a}_{i} \geq 1, & \forall m \in {p^{Met}_{0}} & \\ & v_{i} \in \mathbb{R}, \enspace a_{i} \in \{0,1\} & \forall i \in \text{Rxn} & \\ &\bar{a}_{k} \in \{0,1\} & \forall k \in \text{Rev} & \end{aligned} $$

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 \(\phantom {\dot {i}\!}r_{\mathrm {o}_{2}}\) that transports O 2 into the network. We would need a vector v with \(\phantom {\dot {i}\!}v_{\mathrm {o}_{2}} \geq \delta \) and \(\phantom {\dot {i}\!}v_{\mathrm {o}_{2}} = 0\) 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 \(j \in \mathcal {F}\) we search for a flux vector v j satisfying D j v jd j 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 \({a_{i}^{j}}\) for \({v_{i}^{j}}\). With this, the Eqs. (1) and (2) become

$$ \begin{aligned} \delta \, {a^{j}_{i}} \leq {v^{j}_{i}} \leq M{a^{j}_{i}} ~~~~~~~~~~~~~~~~~~~~~~~~~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall i\in\text{Irrev}, \end{aligned} $$
(4)
$$ \begin{aligned} \delta \, {a^{j}_{i}} - M \, \bar{a}^{j}_{i} \leq {v_{i}^{j}} \leq M \, {a^{j}_{i}} - \delta \, \bar{a}^{j}_{i}, ~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall i\in\text{Rev}, \end{aligned} $$
(5)
$$ \begin{aligned} {a^{j}_{i}} + \bar{a}^{j}_{i} \leq 1 ~~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, ~~~\forall i\in\text{Rev}. \end{aligned} $$
(6)

Using the new variables \({a_{i}^{0}}\), we reformulate the constraints regarding the protected reactions: \( {a^{0}_{i}} = 1\), for all iP Irr and \( {a^{0}_{i}} + \bar {a}^{0}_{i} = 1\), for all iP Rev. Finally, the constraints regarding the protected metabolites become

$$ \sum_{i \in \text{Rxn}_{m}} {a^{0}_{i}} + \sum_{i \in \text{Rev}_{m}} \bar{a}^{0}_{i} \geq 1, \quad \forall m \in {p^{Met}_{0}}. $$
(7)

To obtain a minimum subnetwork, we have to minimize the total number of active reactions. Thus, we need binary variables a i with the property

$$\begin{aligned} a_{i} &= 0~ \text{if and only if}~{a_{i}^{j}}\\ &=0~ \text{for all~} j\in \{ 0, \ldots, |\mathcal{F}|\}, \text{or equivalently} \\ a_{i} &= 1~ \text{if and only if~} {a_{i}^{j}} =1~ \text{for some~} j\in \{ 0, \ldots, |\mathcal{F}|\}. \end{aligned} $$

For irreversible reactions, this can be encoded by the constraints

$$\begin{array}{*{20}l} a_{i} \enspace \leq \enspace \sum_{j=0}^{|\mathcal{F}|} {a_{i}^{j}} \enspace \leq \enspace (1+|\mathcal{F}|) \cdot a_{i}, \quad \forall i\in\text{Irrev}, \end{array} $$
(8)

and for reversible reactions we get

$$\begin{array}{*{20}l}{} a_{i} \enspace \leq \enspace \sum_{j=0}^{|\mathcal{F}|} \left({a_{i}^{j}} + \bar{a}_{i}^{j} \right) \enspace \leq \enspace (2 + 2 |\mathcal{F}|) \cdot a_{i}, \quad \forall i\in\text{Rev}. \end{array} $$
(9)

The resulting MILP is the following:

$$\begin{aligned} &\mathbf{(minNW)}\ \min \sum_{i \in\text{Rxn}} a_{i} & \\ & S v^{j} = 0, \enspace l \leq v^{j} \leq u \qquad\qquad\qquad\qquad\qquad~~~~~~ \forall j\in \{ 0, \ldots, |\mathcal{F}|\}& \\ & D_{j} v^{j} \leq d_{j} \qquad\qquad\qquad\qquad\qquad\qquad \qquad~~~~\forall j\in \{ 1, \ldots, |\mathcal{F}|\} & \\ & \delta \, {a^{j}_{i}} \leq {v^{j}_{i}} \leq M{a^{j}_{i}} \qquad\qquad \qquad\qquad~~~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall i\in\text{Irrev} & \\ &\delta \, {a^{j}_{i}} - M \, \bar{a}^{j}_{i} \leq {v_{i}^{j}} \leq M \, {a^{j}_{i}} - \delta \, \bar{a}^{j}_{i} \qquad \qquad~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall i\in\text{Rev}& \\ & {a^{j}_{i}} + \bar{a}^{j}_{i} \leq 1 \qquad \qquad\qquad \qquad\qquad~~~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall i\in\text{Rev} & \\ & {a^{0}_{i}} = 1, \enspace {a^{0}_{k}} + \bar{a}^{0}_{k} = 1 \qquad \qquad\qquad\qquad~~~~~~~~ \forall i \in P^{\text{Irr}}, \forall k \in P^{\text{Rev}} & \\ & \sum_{i \in \text{Rxn}_{m}} {a^{0}_{i}} + \sum_{i \in \text{Rev}_{m}} \bar{a}^{0}_{i} \geq 1, \qquad\qquad\qquad\qquad\qquad ~~~~~\forall m \in {p^{Met}_{0}} & \\ & a_{i} \enspace \leq \enspace \sum_{j=0}^{|\mathcal{F}|} {a_{i}^{j}} \enspace \leq \enspace (1+|\mathcal{F}|) \cdot a_{i}, \qquad \qquad\qquad\qquad~~~\forall i\in\text{Irrev} & \\ & a_{i} \enspace \leq \enspace \sum_{j=0}^{|\mathcal{F}|} \left({a_{i}^{j}} + \bar{a}_{i}^{j} \right) \enspace \leq \enspace (2 + 2 |\mathcal{F}|) \cdot a_{i}, \qquad\qquad\qquad \forall i\in\text{Rev} & \\ & {v_{i}^{j}} \in \mathbb{R}, \enspace {a^{j}_{i}}, a_{i} \in \{0,1\}, \qquad\qquad\qquad~~~~ \forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall i \in \text{Rxn} & \\ &\bar{a}^{j}_{k} \in \{0,1\} \qquad\qquad\qquad\qquad\qquad~~~~~ \forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall k \in \text{Rev} & \end{aligned} $$

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:

$${} \sum\limits_{i \in\text{Rxn}} \left(1-{Z^{k}_{i}}\right) a_{i} + \sum\limits_{i \in\text{Rxn}} {Z^{k}_{i}} (1-a_{i}) \geq 1, \quad k = 1,2,\dots $$
(10)

where \({Z^{k}_{i}} = 1\) if reaction i carries flux in the subnetwork which was computed in the k-th step, otherwise \({Z^{k}_{i}}=0\). 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 [1922]. We first recall some basic definitions from flux coupling analysis (FCA).

A reaction rRxn is called blocked if v r =0 for all \(v \in C_{0} = \left \{ v\in \mathbb {R}^{\text {Rxn}} \mid Sv = 0, v_{i} \geq 0, \forall i \in \text {Irrev} \right \}\). 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,sRxn, we say r is partially coupled to s, and write rs, if v r =0v s =0, for all vC 0. The relation rs 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]={sRxnrs }. It follows \(\text {Rxn} = \bigcup _{[r]\in \overline {\text {Rxn}}} \; [r]\), where \(\overline {\text {Rxn}}\) 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 rs. 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 rs, 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:

$$\begin{aligned} &\mathbf{{(minNW)_{rep}}} \min \sum_{[r] \in \overline{\text{Rxn}}} |\left[r\right]|~ a_{[r]} & \\ & S v^{j} = 0, l \leq v^{j} \leq u \qquad\qquad\qquad\qquad\qquad~~~~~~~~ \forall j\in \{ 0, \ldots, |\mathcal{F}|\} & \\ & D_{j} v^{j} \leq d_{j} \qquad\qquad\qquad\qquad\qquad\qquad\qquad~~~~~ \forall j\in \{ 1, \ldots, |\mathcal{F}|\} & \\ & \delta \, a^{j}_{[r]} \leq {v^{j}_{s}} \leq Ma^{j}_{[r]} \qquad \qquad\quad~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, [r]\in \overline{\text{Irrev}}, s\in [r] & \\ & \delta \, a^{j}_{[r]} - M \, \bar{a}^{j}_{[r]} \leq {v_{s}^{j}} \leq M \, a^{j}_{[r]} \,-\, \delta \, \bar{a}^{j}_{[r]} ~~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, [r]\in \overline{\text{Rev}}, s\in [r] & \\ & a^{j}_{[r]} + \bar{a}^{j}_{[r]}\leq 1 \qquad\qquad\qquad\qquad~~~~~~\forall j\in \{ 0, \ldots, |\mathcal{F}|\}, {[r]}\in\overline{\text{Rev}} & \\ & a^{0}_{[r]} = 1, a^{0}_{[r']} + \bar{a}^{0}_{[r']} = 1 \qquad\qquad\qquad~~~~ \forall {[r]} \in \overline{P^{\text{Irrev}}}, {[r']} \in \overline{P^{\text{Rev}}} & \\ & \sum_{[r] \in \overline{\text{Rxn}_{m}}} a^{0}_{[r]} + \sum_{[r] \in \overline{\text{Rev}_{m}}} \bar{a}^{0}_{[r]} \geq 1, \qquad\qquad\qquad\qquad~~~~~\forall m \in \overline{p^{Met}_{0}} & \\ & a_{[r]}\leq \sum_{j=0}^{|\mathcal{F}|} a_{[r]}^{j} \leq (|\mathcal{F}|+1) \cdot a_{[r]} \qquad\qquad\qquad\qquad~~~\forall {[r]}\in\overline{\text{Irrev}} & \\ &a_{[r]} \leq \sum_{j=0}^{|\mathcal{F}|} \left(a_{[r]}^{j} + \bar{a}_{[r]}^{j} \right) \leq (2+2 |\mathcal{F}|)\cdot a_{[r]} \qquad\qquad~~~~~~\forall {[r]}\in\overline{\text{Rev}} & \\ &a^{j}_{[r]}, a_{[r]} \in \{0,1\}, {v^{j}_{s}} \in \mathbb{R} \qquad~~~~ \forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall [r] \in \overline{\text{Rxn}}, s \in \text{Rxn} & \\ &\bar{a}^{j}_{[r']} \in \{0,1\} \qquad\qquad\qquad\qquad~~~~~~~ \forall j\in \{ 0, \ldots, |\mathcal{F}|\}, \forall [r'] \in \overline{\text{Rev}} & \end{aligned} $$

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. \(\overline {\text {Irrev}}\) denotes the representatives of the irreversible reactions, and \(\overline {\text {Rev}} \) those of the reversible reactions. Similarly, \(\overline {P^{\text {Irr}}} \) resp. \(\overline {P^{\text {Rev}}}\) is the set of representatives of protected irreversible resp. protected reversible reactions. With \(\overline {p^{Met}_{0}}\) we denote the representatives which include a protected metabolite.

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

$$ \begin{aligned} \sum_{[ r] \in\overline{\text{Rxn}}}\left(1\,-\,Z^{k}_{[r]}\right) a_{[r]} + \sum_{[r] \in\overline{\text{Rxn}}} Z^{k}_{[r]} \left(1\,-\,a_{[r]}\right) \geq 1, \,\,\,\, k = 1,2,\dots \end{aligned} $$
(11)

Using representatives we need only \(|\overline {\text {Rxn}}|\) 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])

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 δ avM a is replaced by a=0v=0,a=1vδ, 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

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

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: Helicobacter pylori 26695.

Table 4 Computational results using indicator variables

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:

  1. 1.

    There are 28 protected reactions.

  2. 2.

    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
figure3

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.

Abbreviations

EFM:

Elementary flux mode

FVA:

Flux variability analysis

FCA:

Flux coupling analysis

LP:

Linear program/programming

MILP:

Mixed integer linear program/programming

References

  1. 1

    O’Brien EJ, Monk JM, Palsson BO. Using genome-scale models to predict biological capabilities. Cell. 2015; 161(5):971–87.

    Article  PubMed  PubMed Central  Google Scholar 

  2. 2

    Lewis NE, Nagarajan H, Palsson BO. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nat Rev Microbiol. 2012; 10(4):291–305.

    CAS  PubMed  PubMed Central  Google Scholar 

  3. 3

    Schuster S, Hilgetag C. On elementary flux modes in biochemical reaction systems at steady state. J Biol Syst. 1994; 2(2):165–82.

    Article  Google Scholar 

  4. 4

    Klamt S, Gilles ED. Minimal cut sets in biochemical reaction networks. Bioinformatics. 2004; 20(2):226–34.

    CAS  Article  PubMed  Google Scholar 

  5. 5

    Erdrich P, Steuer R, Klamt S. An algorithm for the reduction of genome-scale metabolic network models to meaningful core models. BMC Syst Biol. 2015; 9(1):48.

    Article  PubMed  PubMed Central  Google Scholar 

  6. 6

    Pfeiffer T, Sanchez-Valdenebro I, Nu J, Montero F, Schuster S. METATOOL: for studying metabolic networks. Bioinformatics. 1999; 15(3):251–7.

    CAS  Article  PubMed  Google Scholar 

  7. 7

    Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng. 2003; 5(4):264–76.

    CAS  Article  PubMed  Google Scholar 

  8. 8

    Vlassis N, Pacheco MP, Sauter T. Fast reconstruction of compact context-specific metabolic network models. PLoS Comput Biol. 2014; 10(1):e1003424.

    Article  PubMed  PubMed Central  Google Scholar 

  9. 9

    Burgard AP, Vaidyaraman S, Maranas CD. Minimal reaction sets for Escherichia coli metabolism under different growth requirements and uptake environments. Biotechnol Prog. 2001; 17(5):791–7.

    CAS  Article  PubMed  Google Scholar 

  10. 10

    Jonnalagadda S, Srinivasan R. An efficient graph theory based method to identify every minimal reaction set in a metabolic network. BMC Syst Biol. 2014; 8(28):1.

    Google Scholar 

  11. 11

    Edwards JS, Palsson BO. Robustness analysis of the Escherichia coli metabolic network. Biotechnol Prog. 2000; 16(6):927–39.

    CAS  Article  PubMed  Google Scholar 

  12. 12

    Segre D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA. 2002; 99(23):15112–7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13

    Wilhelm T, Behre J, Schuster S. Analysis of structural robustness of metabolic networks. Syst Biol. 2004; 1(1):114–20.

    CAS  Article  Google Scholar 

  14. 14

    Behre J, Wilhelm T, von Kamp A, Ruppin E, Schuster S. Structural robustness of metabolic networks with respect to multiple knockouts. J Theor Biol. 2008; 252(3):433–41.

    CAS  Article  PubMed  Google Scholar 

  15. 15

    Acuna V, Chierichetti F, Lacroix V, Marchetti-Spaccamela A, Sagot MF, Stougie L. Modes and cuts in metabolic networks: complexity and algorithms. Biosystems. 2009; 95(1):51–60.

    Article  PubMed  Google Scholar 

  16. 16

    Burgard AP, Pharkya P, Maranas CD. Optknock: A bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng. 2003; 84(6):647–57.

    CAS  Article  PubMed  Google Scholar 

  17. 17

    Haus UU, Klamt S, Stephen T. Computing knock-out strategies in metabolic networks. J Comput Biol. 2008; 15(3):259–68.

    CAS  Article  PubMed  Google Scholar 

  18. 18

    Tamura T, Takemoto K, Akutsu T. Finding minimum reaction cuts of metabolic networks under a Boolean model using integer programming and feedback vertex sets. Comput Knowl Disco Bioinformatics Res. 2012; 1:240–258.

    Article  Google Scholar 

  19. 19

    Burgard AP, Nikolaev EV, Schilling CH, Maranas CD. Flux coupling analysis of genome-scale metabolic network reconstructions. Genome Res. 2004; 14(2):301–12.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  20. 20

    Larhlimi A, David L, Selbig J, Bockmayr A. F2C2: a fast tool for the computation of flux coupling in genome-scale metabolic networks. BMC Bioinforma. 2012; 13(1):57.

    Article  Google Scholar 

  21. 21

    Goldstein Y, Bockmayr A. Double and multiple knockout simulations for genome-scale metabolic network reconstructions. Algorithm Mol Biol. 2015; 10:1.

    Article  Google Scholar 

  22. 22

    Röhl A, Goldstein Y, Bockmayr A. EFM-Recorder - faster elementary mode enumeration via reaction coupling order. In: Strasbourg Spring School on Advances in Systems and Synthetic Biology: 2015. p. 91–100.

  23. 23

    CPLEX. http://www-01.ibm.com/software/commerce/optimization/cplex-optimizer/.

  24. 24

    BiGG Models. http://bigg.ucsd.edu/.

  25. 25

    Marashi SA, Bockmayr A. Flux coupling analysis of metabolic networks is sensitive to missing reactions. BioSystems. 2011; 103:57–66.

    Article  PubMed  Google Scholar 

Download references

Acknowledgments

We would like to thank Alexandra Reimers and Yaron Goldstein for constructive discussions.

Funding

Not applicable.

Availability of data and materials

The dataset supporting the conclusions of this article is available in the minimal network repository, Project name: minimal network, Project home page: https://sourceforge.net/projects/minimalnetwork/, Archived version: 2016-07-27, Operating system(s): Platform independent, Programming language: MATLAB, Other requirements: no, License: GNU GPL, Any restrictions to use by non-academics: no.

Authors’ contributions

AR developed the theoretical part, implemented the algorithms and drafted the manuscript. AB supervised the study and participated in writing the document. Both authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Author information

Affiliations

Authors

Corresponding author

Correspondence to Annika Röhl.

Additional file

Additional file 1

Comparison with NetworkReducer, FASTCORE and new experiments. Detailed description of the requirements on the networks and the bounds on the flux rates. (PDF 120 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Röhl, A., Bockmayr, A. A mixed-integer linear programming approach to the reduction of genome-scale metabolic networks. BMC Bioinformatics 18, 2 (2017). https://doi.org/10.1186/s12859-016-1412-z

Download citation

Keywords

  • Constraint-based modeling
  • Model reduction
  • Stoichiometric models
  • Mixed-integer linear programming
  • Metabolic networks