BACKGROUND
the depletion of fossil fuels and increasing concerns over environmental changes are key driving factors for the development of sustainable bioprocesses to produce chemicals and fuels from renewable resources  <cit> .
today, bioprocesses using microorganisms are being increasingly used for production of compounds with applications in food, agriculture, chemical and pharmaceutical industries  <cit> .
bioprocesses provide several advantages over traditional chemical processes including high specificity, low temperature, low pressure and reduced use of strong solvents; thus they are environmentally friendlier while reducing the dependency on fossil resources.
despite these advantages, the industry has not adapted bioprocesses extensively, because the viability of bioprocesses is often questionable due to low yield and productivity for desired compounds  <cit> .
in order to make bioprocesses economically viable, it is essential to engineer microbial strains that offer enhanced yield of the desired product  <cit> .
synthetic biology provides the tools and techniques to design and construct artificial cells with minimal functionality containing a minimal genome, but with all the essential genes for survival in a defined environment and possessing replication capabilities  <cit> .
such minimal cells provide a platform for efficient production of desired chemicals and decontamination of waste streams  <cit> .
strains with reduced genomes have been created by deleting large number of non-essential genes  <cit> .
these strains have shown to have equal or better growth performance compared to their parent strains  <cit> .
in biotechnology applications, improved performance has been reported by the strains with minimal metabolism created by blocking handful of reactions that drive the metabolic flux through the predefined minimal metabolic reactions  <cit> .
burgard et al.
<cit>  proposed a mathematical programming approach to find the minimal reaction sets under different uptake environments.
their study finds that minimal reaction sets are strongly dependent on medium constituents and cellular objectives.
this approach does not provide any indication on what reactions have to be blocked in order to construct the cell with minimal metabolism besides its computational complexity is high.
the approach used by trinh el at.
<cit>  identifies the reactions to be blocked to design the cell with minimal metabolism by considering the reduction in elementary flux modes  achieved by removing reactions from the metabolic model.
efms represent the various independent pathways available for the cell to achieve its cellular objectives  <cit> .
efms analysis has so far been employed only with small metabolic models representing the central metabolism but the computational complexity of efms analysis prevents its application to genome-scale models.
we have previously proposed a graph theory based approach for identifying minimal reaction set in metabolic networks  <cit> .
the approach exploits the network structure of metabolic networks and uses math programming efficiently, to identify the minimal reaction set.
significant reduction in the computational time has been achieved using the graph theory based approach compared to classical math programming.
the presence of redundant pathways in metabolic networks results in alternate optimal solutions and consequently create mismatch between model predictions and experimental observations  <cit> .
redundant pathways also lead to multiple minimal reaction sets with different biological significance.
several factors related to the physical and biological functioning of the designed cells including high substrate utilization, deregulated pathways, high tolerance to inhibitors, robust reproduction, predictable metabolic interactions, and physical robustness to sustain the stress and strain during fermentation have to be considered before creating minimal cells  <cit> .
though not all these factors may be equally important in designing minimal metabolic cells, some such as practically achievable metabolic fluxes, thermodynamically favourable pathways, and high substrate utilization should be incorporated.
also, the number of reactions to be knocked-out in order to create a minimal metabolism cell is an important factor.
each solution, minimal metabolism cell, found through computational analysis of metabolic network has different properties that may not be captured in the model used for computational analysis.
identifying all the minimal reaction sets would enable us to evaluate such non-quantifiable properties of different minimal reaction sets and select the one most suitable for experimental development.
in this paper, we propose a graph theory based recursive math programming approach to identify all the minimal reaction sets in the metabolic network.
methods
milp for finding minimal reaction set
the metabolic network of a given microorganism with n metabolites and m reactions is mathematically represented as  <cit> :
∑j=1msijvj= <dig> i= <dig> ,…,n
the stoichiometric matrix s captures interactions among reactions where sij is the stoichiometric coefficient of the ith metabolite in the jth reaction and vj is the flux  of reaction j. the zero in the right hand side is due to the steady-state assumption generally considered in metabolic network analysis.
the mathematical representation of metabolic networks enables analysis of the metabolism using optimization methods to identify internal flux distribution, metabolic capabilities, and strain improvement strategies through gene knock-out or insertion of non-native reactions  <cit> .
identification of minimal reaction set can be represented as an optimization problem given by  <cit> :
minimizez=∑j=1myjs.t∑j=1msijvj=0i= <dig> ,…,nvjmin·yj≤vj≤vjmax·yjj= <dig> ,…,myj∈ <dig> j= <dig> ,…,mvbiomass≥vbiomassmaxvj∈r
here, vjmin and vjmax represent the lower and upper bounds for the flux through reaction j. a binary variable yj is associated with each reaction with ‘1’ indicating the presence/activation of the reaction and ‘0’ its absence/deactivation.
cellular objectives are incorporated as constraints, for example, the objective in eq.
is to produce at least νmaxbiomass biomass.
although the mixed-integer linear programming  in eq.
has been reported to be successfully solved in some cases, the computational time increases exponentially with number of reactions  <cit> .
we previously reported an efficient approach that combines graph theory with math programming to solve this problem .
in our hybrid approach, a metabolic network is considered as an and-or graph where nodes represent metabolites and arcs represent reactions.
reactions that require multiple metabolites to proceed are considered to be related by a and-logic, while reactions that can produce or consume a metabolite using independent routes are considered to be conjoined by an or logic.
a depth is associated with each node and arc in the network starting with the extracellular metabolites and primary uptake reactions which are deemed to be of depth  <dig>  the depth of every other metabolite and reaction is assigned as an increment over its predecessor’s.
there are two phases in the hybrid approach .
based on the depth of reactions, phase  <dig> decomposes the metabolic network into sub-networks which are then analyzed in isolation using small milps to classify reactions as essential, extraneous or indeterminate.
essential reactions  are required for the cell to meet biological objectives and hence they are the part of every minimal set.
extraneous reactions  are not necessary for the cell.
indeterminate reactions  primarily consist of substitutable reactions i.e.
reactions which can be substituted with other reactions to achieve the cellular objectives.
these irs are holistically analysed in the subsequent phase  <dig>  using a milp with the same structure as that in  but smaller than the monolithic one.
through this, a subset of irs called additional reactions  necessary for the minimal metabolism cell are identified which together with srs identified in phase  <dig> form the minimal reaction set.
a substantial reduction in the computational time  required to identify one solution could be achieved through the hybrid approach compared to solving eq.
directly.
in this paper, we extend the above hybrid approach to identify all minimal reaction sets.
the theoretical basis of the proposed approach is discussed first.
reaction dependency and grouping
reactions in the metabolic network are dependent on each other since the network is an interconnected system designed to achieve the biological objectives of the cell.
two different kinds of dependencies can be identified – linear and flux dependency.
linear dependency arises between reactions due to the structure of the network where the product of a reaction feed into exactly one other reaction.
when a set of reactions are all linearly dependent, they can be considered to form a linear pathway.
instances of several reactions forming a linear pathway in metabolic networks are common.
for example, in the sample metabolic network shown in figure 2a, two external metabolites a_ext and b_ext enter into the cell and biomass is produced from them through the reaction network.
two linear pathways can be identified in this sample network {r <dig> r <dig>  r <dig> r5} and {r <dig> r6}.
in linear pathways, under steady-state assumption mentioned above, the flux through all the reactions has to be equal.
hence, deletion of any reaction in the linear pathway would result in the deletion of the whole pathway.
another kind of dependency, flux dependency, exists between reactions which are not structurally linear, but are required to co-exist to balance the fluxes.
if a reaction produces two products which are then consumed by two different reactions, then these three reactions are dependent on each other, since at steady state, the metabolites produced by the first reaction have to be consumed by the two down-stream reactions.
deletion of any one of these reactions would make all other reactions incapable of carrying flux at steady state.
in the sample network, although reactions {r <dig>  r <dig>  r11} are not linearly dependent, they have flux dependencies, since reactions {r <dig>  r11} consume the products of r <dig> and are dependent on each other through the flux balance requirement.
the linear and flux dependency among reactions in a metabolic network can be exploited to assemble reactions into groups for analysis, rather than analyzing them individually.
since deletion of any reaction would force all the dependent reactions to be excluded from the network, for network optimization using a milp, it is sufficient to associate a single binary variable with each group of reactions.
reduction of the number of binary variables reduces the search space and consequently reduces the computational cost of finding solutions.
identification of dependent reactions and simplification of metabolic networks using the reaction dependency has been reported in literature .
groups of dependent reactions are generally identified by comparing the rows in the null-space matrix of the stoichiometric matrix, s. the null-space represents all the possible steady-state flux distributions that satisfy eq.
<dig> and the dependent reactions are the rows in this matrix with same values after normalization with no contradictions in the directionality for irreversible reactions.
since this procedure depends on the directionality of reactions without considering the structural features of metabolic netwotk, it may not identify some dependent reaction groups due to imperfect assignment of reaction directionality in the metabolic networks.
also, identification of dependent reaction groups strictly based on flux distributions results in groups of structurally unrelated reactions which hinders interpretation.
we have developed a graph based algorithm that exploits the structure of metabolic network to identify groups of dependent reactions as described next.
given a metabolic network where the depth has been assigned to reactions as described in milp for finding minimal reaction set, the algorithm first creates a list of reactions, sorted in the ascending order of their depth.
reactions are grouped from this list using an iterative procedure.
in each iteration, a new group is created starting with the first reaction in the list, i.e.
reaction with the lowest depth.
dependent reactions are then added to this group step-by-step by searching the metabolic network in a breath-first manner.
in each step, all the reactions at the current depth + <dig> are collected and tested for linear or flux dependency.
reactions are added to the group if they dependent on other reactions that are already present in the group.
specifically, a single reaction that receives its reactants exclusively from another reaction in that group is deemed as linearly dependent.
similarly, multiple reactions are deemed to have flux dependency if all their reactants originate from one reaction already in the group.
the search continues until al the reactions are evaluated i.e., the highest depth in the network is reached or no dependent reaction can be found at a given depth.
once a group of dependent reactions is identified, all these reactions are removed from the reaction list.
in the subsequent iteration, the algorithm continues with the creation of a new group with the first reaction in the updated reaction list.
the algorithm terminates when the reaction list becomes empty.
we illustrate the algorithm using the sample network shown in figure 2a.
in the first iteration, the algorithm starts a new group with reaction r <dig> , identifies reaction r <dig> as linearly dependent at depth  <dig> in the first step, and then reaction r <dig> in the second step, and r <dig> in the third step.
in the fourth step, r <dig> is added to the group since one of its two reactants is exclusively from r <dig>  reactions {r <dig> r11} are then identified as having flux dependency since they exclusively receive their reactants from r <dig>  the search stops at this step since there are no further reactions in the list at higher depths.
thus the first group of dependent reactions is {r <dig> r <dig>  r <dig> r <dig> r <dig> r <dig> r11}.
the graph based algorithm thus groups reactions based on both linear and flux dependencies.
the number of reactions in a group is called its norm.
thus, the norm of this group is  <dig>  all these  <dig> reactions are removed from the reaction list.
the second iteration starts with reaction r <dig> , and identifies reaction r <dig> as its dependent.
the search for dependent reactions stops here since two reactions, {r <dig>  r8}, consume the product  of r <dig>  hence, the second group has  <dig> reactions {r <dig> r6}.
continuing in this fashion, two single reaction groups {r7} and {r8} are also identified.
thus, in total, there are four different groups in the sample network as shown in figure 2b.
once the groups of dependent reactions have been identified in the metabolic network, analysis can be carried out on these groups rather than on the individual reactions.
if a reaction from a group is essential for the cell, all the reactions in that group become essential since they are dependent on each other.
similarly, all the reactions in the group become extraneous or indeterminate if one of the reactions in the group is extraneous or indeterminate, respectively.
hence, the minimal reaction set identification problem is reduced to identification essential reaction groups , extraneous reaction groups , and indeterminate reaction groups .
recursive milp for finding all minimal reaction sets
the proposed recursive milp approach for identifying all the minimal reaction sets in metabolic network is shown in figure  <dig>  a given metabolic network is described using the groups of dependent reactions where a single binary variable is associated with each group.
then, phase  <dig> of the proposed approach classifies these groups into essential, extraneous and indeterminate groups using the algorithm described in jonnalagadda et al.
<cit> .
as described in milp for finding minimal reaction set, the essential reaction groups  are necessary for the cell to meet its cellular objectives and hence these groups have to be present in all minimal reaction sets.
extraneous reaction groups  are unnecessary for the cell and will be absent in every minimal reaction set.
indeterminate reaction groups  comprise substitutable reactions  which are the source of multiple optimal solutions hence, all minimal reaction sets can be identified by finding all the different additional reaction groups  from the irgs.
these multiple sets of args together with the srgs identified in phase  <dig> forms all the minimal reaction sets.
the algorithm for finding all args from irgs formulated as a recursive milp as shown in figure  <dig>  the first set of arg is found by solving the milp with the same constraints as given in eq.
<dig>  but considering only the irgs where binary variables have been associated for each group .
the objective function for the optimization is the minimization of ∑lirgwl·yl where the wl is the norm of the group and yl is the binary variable associated with that group.
the optimization procedure thus will identify the args such that the total number of reactions is minimal.
the args together with the srg from phase  <dig> forms the first minimal reaction set.
once an optimal solution is found, a constraint is added to the model to exclude that solution from the search space .
based on lee et al.
, the following constraint is added to eq.
2:
∑r∈nzwr·yr≤nz− <dig>
where nz is the groups in the optimal solution, and yr is the binary variable associated with the groups in nz.
eq.
means that at least one of the non-zero binary variable in the optimal solution is set to zero.
hence in the next recursion, nz is excluded from the search space and the optimizer is forced to find a new optimal solution.
this recursive procedure terminates when the optimizer returns a sub-optimal solution, i.e.
a solution with more reactions than that in the first solution.
in principle, all the minimal reaction sets can be identified by recursively solving the milp with a new constraint added to the model in each recursion.
however, the grouping of reactions offers insights which enable the computational cost to be reduced significantly by generating additional solutions without solving the milp, using group substitutability analysis .
group substitutability analysis for identifying solutions
in metabolic networks, some reactions share similar cellular functions such as producing, consuming a metabolite, or recycling co-factors.
for example, in the sample network shown in figure 2a reactions r <dig> and r <dig> consume metabolite d and produce metabolite g. the presence of reactions with similar functions enables the cell to survive under different conditions, stress, and malfunction of genes through substitution of reaction for another inactive reaction.
these reactions are considered substitutable since they result in alternate optima.
in the minimal reaction set identification problem, substitutable reactions lead to multiple minimal reaction sets.
the above recursive milp approach can be employed to identify all minimal reaction sets.
alternatively, many candidate solutions can be generated more efficiently by simply substituting a reaction in an optimal solution.
in this paper, we perform this substitution analysis on groups to efficiently identify alternate optimal solutions.
two types of group substitution are possible – single and multi-group substitution.
in single group substitution, a group is substituted with another group of the same norm.
so, the total number of reactions in the optimal solution remains unchanged.
for example, in figure 2b groups  <dig> and  <dig> are substitutable as both have the same metabolic function and have a norm of  <dig>  thus if group  <dig> is present in an optimal solution, another candidate solution can be generated by replacing it with group  <dig>  substitutability for single groups can be identified easily since all the groups would produce and consume the same metabolites, and can hence be identified by or gates in the graph representation of the metabolic network.
sets of groups could also be analysed for substitutability but this multi-group substitution is computationally complex and is beyond the scope of this paper.
group substitutability analysis is conducted in the proposed approach following the identification of a solution by the milp and candidate solutions generated.
not every candidate solution identified by the qualitative approach would meet the cellular objectives.
therefore, it is essential to verify the candidate solutions to ensure that the predefined biological objectives are satisfied.
this verification is conducted by solving a linear program  with the objective of maximizing the cellular objective , which is computationally efficient.
only candidate solutions that satisfy the objective are deemed as optimal solutions to the original milp and appended to the set of optimal solutions .
other candidate solutions are discarded.
a new constraint is also added to the model for each optimal solution thus identified to eliminate their re-identification in future recursions .
the algorithm then continues with solving the milp to find other optimal solutions.
RESULTS
we illustrate the proposed method by identifying all minimal reaction sets that support predefined growth for two systems – escherichia coli and saccharomyces cerevisiae.
case study 1: aerobic growth of escherichia coli on glucose
here, we identify all the minimal reaction sets from the e. coli metabolic network so as to meet cellular objective νmaxbiomass≥ <dig>  g/gdw∙h for a glucose uptake rate of 10 mmol /gdw∙h.
the network contains  <dig> metabolites and  <dig> reactions  <cit> .
these  <dig> reactions are first grouped based on dependency as described in reaction dependency and grouping.
there are in total  <dig> groups —  <dig> groups with norm  <dig>   <dig> groups with norm  <dig> and  <dig> groups with norm  <dig> i.e.
single reaction groups.
hence, the number of binary variables required for milp is reduced from  <dig> to  <dig>  the proposed recursive milp method is then employed to identify all the minimal reaction sets from this network.
phase  <dig> of the proposed approach classified  <dig> groups  as essential reaction groups of which  <dig> groups with norm  <dig> the remaining  <dig> groups of norm  <dig>  the method also denoted  <dig> groups  as extraneous.
there were  <dig> groups  identified as indeterminate containing  <dig> group with norm  <dig>   <dig> groups with norm  <dig> and the remaining  <dig> with norm  <dig>  hence, the number of binary variables defined in phase  <dig> is reduced from  <dig> to  <dig>  then the recursive milp is employed to identify irgs.
the first optimal solution contains  <dig> groups  —  <dig> groups of norm  <dig> and  <dig> group with norm  <dig>  the recursive milp found two more optimal solutions  that meet the predefined cellular objective.
in the fourth iteration, the optimizer found a sub-optimal solution with  <dig> ars and hence is terminated.
these three sets of additional reactions together with the  <dig> srs from phase  <dig> form the three different minimal reaction sets.
to cross-validate the results, we also implemented the classical monolithic milp approach with  <dig> binary variables.
the monolithic milp also identified the same three minimal reaction sets thus confirming the accuracy of the proposed approach.
the three minimal reaction sets identified in the e. coli metabolic network are shown in figure  <dig>  the reactions in the minimal reaction set are shown by thick solid line.
the three minimal reaction sets differ from each other by the presence of a single unique reaction while  <dig> of  <dig> reactions in the minimal reaction set are common to all three.
this indicates that  <dig> out of the  <dig> reactions identified in phase  <dig> by recursive milp are common to all three minimal sets.
however, these  <dig> reactions are deemed as indeterminate  in phase  <dig> since there exist alternative  pathways.
minimal reaction set  <dig> has a unique reaction phosphoenolpyruvate carboxykinase  while set  <dig> has pyruvate kinase  and set  <dig> has transhydrogenase .
the comparison of flux distributions from the different reaction sets reveals how the cell meets its biological objective while still staying minimal.
for example, minimal reaction set  <dig> contains ppck which converts oxaloacetate, produced from phosphoenol pyruvate through phosphoenolpyruvate corboxylase  reaction, back to phosphoenol pyruvate forming a cycle.
since such cycles may not generally be active at steady-state, considering thermodynamics  <cit> , this minimal reaction set may not be suitable for developing minimal metabolism cell.
similarly, minimal reaction set  <dig> has large flux through the transhydrogenase reaction that regenerates cofactors nadh, nadp from nad and nadph.
this set is also not desirable for developing minimal metabolism cell since such a high flux may not be practically possible in the organism.
in comparison, set  <dig> has a unique reaction pyk that converts phosphoenolpyruvate to pyruvate which is part of glycolysis pathway in aerobically growing e.coli and contains no coupled reactions ; hence, it is a suitable reaction set for developing the minimal metabolism cell.
we identify the number of reactions to be knocked-out from e.coli in order to develop minimal metabolism cell based on each minimal reaction set using the graph theory based approach  <cit> .
in brief, the procedure iteratively selects the reaction with lowest depth from the list of reactions not present in the minimal reaction set as the knock-out candidate.
in each iteration, all the reactions dependent on the selected reaction are excluded from the list.
the procedure continues until the list of reactions becomes empty.
thus all the reactions selected in this procedure have to be removed from the strain to achieve the minimal metabolism based on the selected minimal reaction set.
in this case study, all three minimal reaction sets require  <dig> reactions to be knocked-out from escherichia coli;  <dig> of these  <dig> reactions are the same in all cases.
hence, the minimal metabolism cell can be constructed by suitable blocking out the reactions corresponding to minimal reaction set  <dig>  in summary, finding multiple minimal sets enables us to develop the best minimal metabolism cell by selectively deleting the remaining two reactions.
case study 2: aerobic growth of saccharomyces cerevisiae on glucose
we now illustrate the computational efficiency of the proposed method by identifying all the minimal reaction sets for a genome-scale model of saccharomyces cerevisiae containing  <dig> metabolites and  <dig> reactions  <cit> .
the cellular objective is selected as νmaxbiomass of  <dig>  g /gdw∙h for glucose uptake rate of 1 mmol /gdw∙h.
the model is reduced to  <dig> reactions by removing  <dig> reactions that are not connected to the glycolysis pathway and  <dig> reactions which differ in a cofactor.
there are  <dig> groups of dependent reactions in the norm range  <cit> .
phase  <dig> of the proposed approach identified  <dig> groups  as essential,  <dig> groups  as extraneous, and  <dig> groups  as indeterminate.
the extraneous reactions are removed from further analysis.
unlike the e. coli model, the saccharomyces cerevisiae model has compartments.
out of the  <dig> indeterminate reactions,  <dig> reactions are involved in transporting metabolites among the compartments and inter-converting cofactor metabolites.
these are deemed to be essential reactions.
the remaining  <dig> groups containing  <dig> indeterminate reactions are further analyzed in phase  <dig> using recursive milp to find all additional reactions.
the results are given in table  <dig>  there are  <dig> reactions in the first solution that together with the  <dig> essential reactions from phase  <dig> form the first minimal reaction set with  <dig> reactions.
saccharomyces cerevisiae
case study
the proposed approach identified  <dig> minimal reaction sets within a reasonable amount of time.
based on the first solution,  <dig> other minimal reaction sets that meet the predefined cellular objective are identified through group substitutability analysis leading to a total of  <dig> minimal reaction sets in the first iteration.
these  <dig> optimal solutions were excluded from the search space through addition of new constraints.
in the seconds iteration,  <dig> more minimal reaction sets were identified —  <dig> from milp and  <dig> from substitution analysis.
the algorithm then continues with next iteration.
the results for each iteration are shown in table  <dig>  there are a total of  <dig> minimal reaction sets for this metabolic network.
the proposed recursive milp approach has to go through  <dig> iterations to identify all these optimal solutions.
further execution of milp resulted in a sub-optimal solution with  <dig> reactions, hence it terminated.
saccharomyces cerevisiae
case study for all the iterations in phase 2
to quantify the improvement achieved, we executed the milp with  <dig> binary variables for  <dig> indeterminate reactions.
the solver required  <dig> , <dig> cplex iterations and 1038 seconds to find the optimal solution.
the reduction of number of binary variables has resulted in a significant improvement with approximately 60% reduction in cplex iterations and 80% reduction in the time required for finding the optimal solution in phase  <dig>  we also compared the total time required for phase  <dig> &  <dig> to find the first solution by the proposed method with monolithic milp and graph theory based approach without grouping.
the results are given in table  <dig>  the time required by the proposed method is ~ 4% and 22% of the time required for the monolithic milp and graph theory based approach, accordingly.
for all  <dig> solutions, the proposed approach required 16311 seconds.
the large computational time required for monolithic milp restrained its use for finding all minimal reaction sets.
nonetheless, to validate the results monolithic milp was employed after excluding all the  <dig> minimal reactions from search space.
it found a sub-optimal solution with  <dig> reactions in the minimal set.
this guarantees that the proposed method identifies all minimal reaction sets.
discussion and
CONCLUSIONS
development of cells with minimal metabolic functionality is increasingly gaining importance.
the presence of redundant reactions in metabolic networks results in multiple minimal reactions sets that can meet the predefined cellular objectives.
in this paper, we proposed a graph theory augmented recursive milp approach to identify all the minimal reaction sets in a metabolic network.
the proposed method has been demonstrated by finding all the minimal reaction sets for escherichia coli and saccharomyces cerevisiae.
the proposed approach correctly identified all the minimal reaction sets in both the cases.
we also proposed the concept of grouping dependent reactions to reduce the number of binary variables for milp formulation.
in the present study, several groups of dependent reactions are identified in escherichia coli and saccharomyces cerevisiae and exploited to reduce the number of binary variables and consequently the solution time.
since the use of binary variables is very common in metabolic network analysis for identifying strain improvement strategies  <cit> , the reaction group concept will benefit the other applications as well.
here, we have developed a graph based algorithm that exploits the structure of the metabolic network to identify groups of dependent reactions.
we now compare the groups of dependent reactions identified by the proposed method with the previously reported approach based on steady-state flux distribution.
we used the metatool software  <cit>  to find the dependent reaction groups using steady-state flux distribution.
while the proposed graph based approach found  <dig> reaction groups  with norm more than  <dig> in the yeast model used in case study  <dig>  the flux based approach identified  <dig> dependent reaction groups .
the norm range of the reaction groups identified by the proposed approach is  <cit>  while that for the flux based approach is  <cit> .
the lower range of the proposed approach could be due to the strict use of the structure of the metabolic network in contrast to the flux based approach which can place structurally unconnected reactions in a group.
nonetheless, the use of structure of metabolic network identified some dependent reaction groups that were not identified by the flux based approach due to the use of reaction directionality as described in reaction dependency and grouping.
while there is a large overlap between the dependent reaction groups found by the two approaches, several groups are uniquely identified by each approach.
hence, combining the two approaches could further simplify the metabolic network.
generally, different minimal reaction sets have different structural and functional properties in terms of robustness, predictable metabolic interactions, practically achievable metabolic fluxes, and thermodynamically non-favourable pathways and cycles.
we identified one such cycle in one of the solutions identified in e.coli model in case study  <dig> and very high flux through co-factor recycling reaction in another solution.
thus, of the three, only one solution has practically achievable fluxes without any coupled reactions  and is suitable for developing the minimal metabolism cell.
this clearly shows the importance of identifying all minimal reaction sets.
in case study  <dig>  the proposed method identified  <dig> solutions.
while some of them use glycolysis pathway to produce dhap others use pentose phosphate pathway.
other differences observed include use of different reactions for the consumption of dhap and formate, recycling of nad, nadh, nadp, and nadph, and production of ammonium.
a preliminary analysis shows that some minimal reaction sets require complex set of reactions with impractical high fluxes compared to others.
further analysis may reveal the details of these pathways that qualify them to developing into minimal metabolism cells.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
both sj and rs contributed to the concept and methodology development.
sj implemented the methodology and conducted the analysis of results and biological interpretation.
rs supervised the study and assisted in implementation.
sj drafted the manuscript.
both authors read and approved the final manuscript.

 the depletion of fossil fuels and increasing concerns over environmental changes are key driving factors for the development of sustainable bioprocesses to produce chemicals and fuels from renewable resources. today, bioprocesses using microorganisms are being increasingly used for production of compounds with applications in food, agriculture, chemical and pharmaceutical industries. bioprocesses provide several advantages over traditional chemical processes including high specificity, low temperature, low pressure and reduced use of strong solvents; thus they are environmentally friendlier while reducing the dependency on fossil resources. despite these advantages, the industry has not adapted bioprocesses extensively, because the viability of bioprocesses is often questionable due to low yield and productivity for desired compounds. in order to make bioprocesses economically viable, it is essential to engineer microbial strains that offer enhanced yield of the desired product. synthetic biology provides the tools and techniques to design and construct artificial cells with minimal functionality containing a minimal genome, but with all the essential genes for survival in a defined environment and possessing replication capabilities. such minimal cells provide a platform for efficient production of desired chemicals and decontamination of waste streams. strains with reduced genomes have been created by deleting large number of non-essential genes. these strains have shown to have equal or better growth performance compared to their parent strains. in biotechnology applications, improved performance has been reported by the strains with minimal metabolism created by blocking handful of reactions that drive the metabolic flux through the predefined minimal metabolic reactions. burgard et al. proposed a mathematical programming approach to find the minimal reaction sets under different uptake environments. their study finds that minimal reaction sets are strongly dependent on medium constituents and cellular objectives. this approach does not provide any indication on what reactions have to be blocked in order to construct the cell with minimal metabolism besides its computational complexity is high. the approach used by trinh el at. identifies the reactions to be blocked to design the cell with minimal metabolism by considering the reduction in elementary flux modes achieved by removing reactions from the metabolic model. efms represent the various independent pathways available for the cell to achieve its cellular objectives. efms analysis has so far been employed only with small metabolic models representing the central metabolism but the computational complexity of efms analysis prevents its application to genome-scale models. we have previously proposed a graph theory based approach for identifying minimal reaction set in metabolic networks. the approach exploits the network structure of metabolic networks and uses math programming efficiently, to identify the minimal reaction set. significant reduction in the computational time has been achieved using the graph theory based approach compared to classical math programming. the presence of redundant pathways in metabolic networks results in alternate optimal solutions and consequently create mismatch between model predictions and experimental observations. redundant pathways also lead to multiple minimal reaction sets with different biological significance. several factors related to the physical and biological functioning of the designed cells including high substrate utilization, deregulated pathways, high tolerance to inhibitors, robust reproduction, predictable metabolic interactions, and physical robustness to sustain the stress and strain during fermentation have to be considered before creating minimal cells. though not all these factors may be equally important in designing minimal metabolic cells, some such as practically achievable metabolic fluxes, thermodynamically favourable pathways, and high substrate utilization should be incorporated. also, the number of reactions to be knocked-out in order to create a minimal metabolism cell is an important factor. each solution, minimal metabolism cell, found through computational analysis of metabolic network has different properties that may not be captured in the model used for computational analysis. identifying all the minimal reaction sets would enable us to evaluate such non-quantifiable properties of different minimal reaction sets and select the one most suitable for experimental development. in this paper, we propose a graph theory based recursive math programming approach to identify all the minimal reaction sets in the metabolic network. methods milp for finding minimal reaction set the metabolic network of a given microorganism with n metabolites and m reactions is mathematically represented as: ∑j=1msijvj= i=,…,n the stoichiometric matrix s captures interactions among reactions where sij is the stoichiometric coefficient of the ith metabolite in the jth reaction and vj is the flux of reaction j. the zero in the right hand side is due to the steady-state assumption generally considered in metabolic network analysis. the mathematical representation of metabolic networks enables analysis of the metabolism using optimization methods to identify internal flux distribution, metabolic capabilities, and strain improvement strategies through gene knock-out or insertion of non-native reactions. identification of minimal reaction set can be represented as an optimization problem given by: minimizez=∑j=1myjs.t∑j=1msijvj=0i=,…,nvjmin·yj≤vj≤vjmax·yjj=,…,myj∈ j=,…,mvbiomass≥vbiomassmaxvj∈r here, vjmin and vjmax represent the lower and upper bounds for the flux through reaction j. a binary variable yj is associated with each reaction with ‘1’ indicating the presence/activation of the reaction and ‘0’ its absence/deactivation. cellular objectives are incorporated as constraints, for example, the objective in eq. is to produce at least νmaxbiomass biomass. although the mixed-integer linear programming in eq. has been reported to be successfully solved in some cases, the computational time increases exponentially with number of reactions. we previously reported an efficient approach that combines graph theory with math programming to solve this problem. in our hybrid approach, a metabolic network is considered as an and-or graph where nodes represent metabolites and arcs represent reactions. reactions that require multiple metabolites to proceed are considered to be related by a and-logic, while reactions that can produce or consume a metabolite using independent routes are considered to be conjoined by an or logic. a depth is associated with each node and arc in the network starting with the extracellular metabolites and primary uptake reactions which are deemed to be of depth the depth of every other metabolite and reaction is assigned as an increment over its predecessor’s. there are two phases in the hybrid approach. based on the depth of reactions, phase decomposes the metabolic network into sub-networks which are then analyzed in isolation using small milps to classify reactions as essential, extraneous or indeterminate. essential reactions are required for the cell to meet biological objectives and hence they are the part of every minimal set. extraneous reactions are not necessary for the cell. indeterminate reactions primarily consist of substitutable reactions i.e. reactions which can be substituted with other reactions to achieve the cellular objectives. these irs are holistically analysed in the subsequent phase using a milp with the same structure as that in but smaller than the monolithic one. through this, a subset of irs called additional reactions necessary for the minimal metabolism cell are identified which together with srs identified in phase form the minimal reaction set. a substantial reduction in the computational time required to identify one solution could be achieved through the hybrid approach compared to solving eq. directly. in this paper, we extend the above hybrid approach to identify all minimal reaction sets. the theoretical basis of the proposed approach is discussed first. reaction dependency and grouping reactions in the metabolic network are dependent on each other since the network is an interconnected system designed to achieve the biological objectives of the cell. two different kinds of dependencies can be identified – linear and flux dependency. linear dependency arises between reactions due to the structure of the network where the product of a reaction feed into exactly one other reaction. when a set of reactions are all linearly dependent, they can be considered to form a linear pathway. instances of several reactions forming a linear pathway in metabolic networks are common. for example, in the sample metabolic network shown in figure 2a, two external metabolites a_ext and b_ext enter into the cell and biomass is produced from them through the reaction network. two linear pathways can be identified in this sample network {r r r r5} and {r r6}. in linear pathways, under steady-state assumption mentioned above, the flux through all the reactions has to be equal. hence, deletion of any reaction in the linear pathway would result in the deletion of the whole pathway. another kind of dependency, flux dependency, exists between reactions which are not structurally linear, but are required to co-exist to balance the fluxes. if a reaction produces two products which are then consumed by two different reactions, then these three reactions are dependent on each other, since at steady state, the metabolites produced by the first reaction have to be consumed by the two down-stream reactions. deletion of any one of these reactions would make all other reactions incapable of carrying flux at steady state. in the sample network, although reactions {r r r11} are not linearly dependent, they have flux dependencies, since reactions {r r11} consume the products of r and are dependent on each other through the flux balance requirement. the linear and flux dependency among reactions in a metabolic network can be exploited to assemble reactions into groups for analysis, rather than analyzing them individually. since deletion of any reaction would force all the dependent reactions to be excluded from the network, for network optimization using a milp, it is sufficient to associate a single binary variable with each group of reactions. reduction of the number of binary variables reduces the search space and consequently reduces the computational cost of finding solutions. identification of dependent reactions and simplification of metabolic networks using the reaction dependency has been reported in literature. groups of dependent reactions are generally identified by comparing the rows in the null-space matrix of the stoichiometric matrix, s. the null-space represents all the possible steady-state flux distributions that satisfy eq.  and the dependent reactions are the rows in this matrix with same values after normalization with no contradictions in the directionality for irreversible reactions. since this procedure depends on the directionality of reactions without considering the structural features of metabolic netwotk, it may not identify some dependent reaction groups due to imperfect assignment of reaction directionality in the metabolic networks. also, identification of dependent reaction groups strictly based on flux distributions results in groups of structurally unrelated reactions which hinders interpretation. we have developed a graph based algorithm that exploits the structure of metabolic network to identify groups of dependent reactions as described next. given a metabolic network where the depth has been assigned to reactions as described in milp for finding minimal reaction set, the algorithm first creates a list of reactions, sorted in the ascending order of their depth. reactions are grouped from this list using an iterative procedure. in each iteration, a new group is created starting with the first reaction in the list, i.e. reaction with the lowest depth. dependent reactions are then added to this group step-by-step by searching the metabolic network in a breath-first manner. in each step, all the reactions at the current depth + are collected and tested for linear or flux dependency. reactions are added to the group if they dependent on other reactions that are already present in the group. specifically, a single reaction that receives its reactants exclusively from another reaction in that group is deemed as linearly dependent. similarly, multiple reactions are deemed to have flux dependency if all their reactants originate from one reaction already in the group. the search continues until al the reactions are evaluated i.e., the highest depth in the network is reached or no dependent reaction can be found at a given depth. once a group of dependent reactions is identified, all these reactions are removed from the reaction list. in the subsequent iteration, the algorithm continues with the creation of a new group with the first reaction in the updated reaction list. the algorithm terminates when the reaction list becomes empty. we illustrate the algorithm using the sample network shown in figure 2a. in the first iteration, the algorithm starts a new group with reaction r, identifies reaction r as linearly dependent at depth in the first step, and then reaction r in the second step, and r in the third step. in the fourth step, r is added to the group since one of its two reactants is exclusively from r reactions {r r11} are then identified as having flux dependency since they exclusively receive their reactants from r the search stops at this step since there are no further reactions in the list at higher depths. thus the first group of dependent reactions is {r r r r r r r11}. the graph based algorithm thus groups reactions based on both linear and flux dependencies. the number of reactions in a group is called its norm. thus, the norm of this group is all these reactions are removed from the reaction list. the second iteration starts with reaction r, and identifies reaction r as its dependent. the search for dependent reactions stops here since two reactions, {r r8}, consume the product of r hence, the second group has reactions {r r6}. continuing in this fashion, two single reaction groups {r7} and {r8} are also identified. thus, in total, there are four different groups in the sample network as shown in figure 2b. once the groups of dependent reactions have been identified in the metabolic network, analysis can be carried out on these groups rather than on the individual reactions. if a reaction from a group is essential for the cell, all the reactions in that group become essential since they are dependent on each other. similarly, all the reactions in the group become extraneous or indeterminate if one of the reactions in the group is extraneous or indeterminate, respectively. hence, the minimal reaction set identification problem is reduced to identification essential reaction groups, extraneous reaction groups, and indeterminate reaction groups. recursive milp for finding all minimal reaction sets the proposed recursive milp approach for identifying all the minimal reaction sets in metabolic network is shown in figure  a given metabolic network is described using the groups of dependent reactions where a single binary variable is associated with each group. then, phase of the proposed approach classifies these groups into essential, extraneous and indeterminate groups using the algorithm described in jonnalagadda et al.. as described in milp for finding minimal reaction set, the essential reaction groups are necessary for the cell to meet its cellular objectives and hence these groups have to be present in all minimal reaction sets. extraneous reaction groups are unnecessary for the cell and will be absent in every minimal reaction set. indeterminate reaction groups comprise substitutable reactions which are the source of multiple optimal solutions hence, all minimal reaction sets can be identified by finding all the different additional reaction groups from the irgs. these multiple sets of args together with the srgs identified in phase forms all the minimal reaction sets. the algorithm for finding all args from irgs formulated as a recursive milp as shown in figure  the first set of arg is found by solving the milp with the same constraints as given in eq.  but considering only the irgs where binary variables have been associated for each group. the objective function for the optimization is the minimization of ∑lirgwl·yl where the wl is the norm of the group and yl is the binary variable associated with that group. the optimization procedure thus will identify the args such that the total number of reactions is minimal. the args together with the srg from phase forms the first minimal reaction set. once an optimal solution is found, a constraint is added to the model to exclude that solution from the search space. based on lee et al., the following constraint is added to eq. 2: ∑r∈nzwr·yr≤nz− where nz is the groups in the optimal solution, and yr is the binary variable associated with the groups in nz. eq. means that at least one of the non-zero binary variable in the optimal solution is set to zero. hence in the next recursion, nz is excluded from the search space and the optimizer is forced to find a new optimal solution. this recursive procedure terminates when the optimizer returns a sub-optimal solution, i.e. a solution with more reactions than that in the first solution. in principle, all the minimal reaction sets can be identified by recursively solving the milp with a new constraint added to the model in each recursion. however, the grouping of reactions offers insights which enable the computational cost to be reduced significantly by generating additional solutions without solving the milp, using group substitutability analysis. group substitutability analysis for identifying solutions in metabolic networks, some reactions share similar cellular functions such as producing, consuming a metabolite, or recycling co-factors. for example, in the sample network shown in figure 2a reactions r and r consume metabolite d and produce metabolite g. the presence of reactions with similar functions enables the cell to survive under different conditions, stress, and malfunction of genes through substitution of reaction for another inactive reaction. these reactions are considered substitutable since they result in alternate optima. in the minimal reaction set identification problem, substitutable reactions lead to multiple minimal reaction sets. the above recursive milp approach can be employed to identify all minimal reaction sets. alternatively, many candidate solutions can be generated more efficiently by simply substituting a reaction in an optimal solution. in this paper, we perform this substitution analysis on groups to efficiently identify alternate optimal solutions. two types of group substitution are possible – single and multi-group substitution. in single group substitution, a group is substituted with another group of the same norm. so, the total number of reactions in the optimal solution remains unchanged. for example, in figure 2b groups and are substitutable as both have the same metabolic function and have a norm of thus if group is present in an optimal solution, another candidate solution can be generated by replacing it with group substitutability for single groups can be identified easily since all the groups would produce and consume the same metabolites, and can hence be identified by or gates in the graph representation of the metabolic network. sets of groups could also be analysed for substitutability but this multi-group substitution is computationally complex and is beyond the scope of this paper. group substitutability analysis is conducted in the proposed approach following the identification of a solution by the milp and candidate solutions generated. not every candidate solution identified by the qualitative approach would meet the cellular objectives. therefore, it is essential to verify the candidate solutions to ensure that the predefined biological objectives are satisfied. this verification is conducted by solving a linear program with the objective of maximizing the cellular objective, which is computationally efficient. only candidate solutions that satisfy the objective are deemed as optimal solutions to the original milp and appended to the set of optimal solutions. other candidate solutions are discarded. a new constraint is also added to the model for each optimal solution thus identified to eliminate their re-identification in future recursions. the algorithm then continues with solving the milp to find other optimal solutions.  we illustrate the proposed method by identifying all minimal reaction sets that support predefined growth for two systems – escherichia coli and saccharomyces cerevisiae. case study 1: aerobic growth of escherichia coli on glucose here, we identify all the minimal reaction sets from the e. coli metabolic network so as to meet cellular objective νmaxbiomass≥ g/gdw∙h for a glucose uptake rate of 10 mmol /gdw∙h. the network contains metabolites and reactions. these reactions are first grouped based on dependency as described in reaction dependency and grouping. there are in total groups — groups with norm groups with norm and groups with norm i.e. single reaction groups. hence, the number of binary variables required for milp is reduced from to the proposed recursive milp method is then employed to identify all the minimal reaction sets from this network. phase of the proposed approach classified groups as essential reaction groups of which groups with norm the remaining groups of norm the method also denoted groups as extraneous. there were groups identified as indeterminate containing group with norm groups with norm and the remaining with norm hence, the number of binary variables defined in phase is reduced from to then the recursive milp is employed to identify irgs. the first optimal solution contains groups — groups of norm and group with norm the recursive milp found two more optimal solutions that meet the predefined cellular objective. in the fourth iteration, the optimizer found a sub-optimal solution with ars and hence is terminated. these three sets of additional reactions together with the srs from phase form the three different minimal reaction sets. to cross-validate the results, we also implemented the classical monolithic milp approach with binary variables. the monolithic milp also identified the same three minimal reaction sets thus confirming the accuracy of the proposed approach. the three minimal reaction sets identified in the e. coli metabolic network are shown in figure  the reactions in the minimal reaction set are shown by thick solid line. the three minimal reaction sets differ from each other by the presence of a single unique reaction while of reactions in the minimal reaction set are common to all three. this indicates that out of the reactions identified in phase by recursive milp are common to all three minimal sets. however, these reactions are deemed as indeterminate in phase since there exist alternative pathways. minimal reaction set has a unique reaction phosphoenolpyruvate carboxykinase while set has pyruvate kinase and set has transhydrogenase. the comparison of flux distributions from the different reaction sets reveals how the cell meets its biological objective while still staying minimal. for example, minimal reaction set contains ppck which converts oxaloacetate, produced from phosphoenol pyruvate through phosphoenolpyruvate corboxylase reaction, back to phosphoenol pyruvate forming a cycle. since such cycles may not generally be active at steady-state, considering thermodynamics, this minimal reaction set may not be suitable for developing minimal metabolism cell. similarly, minimal reaction set has large flux through the transhydrogenase reaction that regenerates cofactors nadh, nadp from nad and nadph. this set is also not desirable for developing minimal metabolism cell since such a high flux may not be practically possible in the organism. in comparison, set has a unique reaction pyk that converts phosphoenolpyruvate to pyruvate which is part of glycolysis pathway in aerobically growing e.coli and contains no coupled reactions; hence, it is a suitable reaction set for developing the minimal metabolism cell. we identify the number of reactions to be knocked-out from e.coli in order to develop minimal metabolism cell based on each minimal reaction set using the graph theory based approach. in brief, the procedure iteratively selects the reaction with lowest depth from the list of reactions not present in the minimal reaction set as the knock-out candidate. in each iteration, all the reactions dependent on the selected reaction are excluded from the list. the procedure continues until the list of reactions becomes empty. thus all the reactions selected in this procedure have to be removed from the strain to achieve the minimal metabolism based on the selected minimal reaction set. in this case study, all three minimal reaction sets require reactions to be knocked-out from escherichia coli; of these reactions are the same in all cases. hence, the minimal metabolism cell can be constructed by suitable blocking out the reactions corresponding to minimal reaction set in summary, finding multiple minimal sets enables us to develop the best minimal metabolism cell by selectively deleting the remaining two reactions. case study 2: aerobic growth of saccharomyces cerevisiae on glucose we now illustrate the computational efficiency of the proposed method by identifying all the minimal reaction sets for a genome-scale model of saccharomyces cerevisiae containing metabolites and reactions. the cellular objective is selected as νmaxbiomass of  g /gdw∙h for glucose uptake rate of 1 mmol /gdw∙h. the model is reduced to reactions by removing reactions that are not connected to the glycolysis pathway and reactions which differ in a cofactor. there are groups of dependent reactions in the norm range. phase of the proposed approach identified groups as essential, groups as extraneous, and groups as indeterminate. the extraneous reactions are removed from further analysis. unlike the e. coli model, the saccharomyces cerevisiae model has compartments. out of the indeterminate reactions, reactions are involved in transporting metabolites among the compartments and inter-converting cofactor metabolites. these are deemed to be essential reactions. the remaining groups containing indeterminate reactions are further analyzed in phase using recursive milp to find all additional reactions. the results are given in table  there are reactions in the first solution that together with the essential reactions from phase form the first minimal reaction set with reactions. saccharomyces cerevisiae case study the proposed approach identified minimal reaction sets within a reasonable amount of time. based on the first solution, other minimal reaction sets that meet the predefined cellular objective are identified through group substitutability analysis leading to a total of minimal reaction sets in the first iteration. these optimal solutions were excluded from the search space through addition of new constraints. in the seconds iteration, more minimal reaction sets were identified — from milp and from substitution analysis. the algorithm then continues with next iteration. the results for each iteration are shown in table  there are a total of minimal reaction sets for this metabolic network. the proposed recursive milp approach has to go through iterations to identify all these optimal solutions. further execution of milp resulted in a sub-optimal solution with reactions, hence it terminated. saccharomyces cerevisiae case study for all the iterations in phase 2 to quantify the improvement achieved, we executed the milp with binary variables for indeterminate reactions. the solver required, cplex iterations and 1038 seconds to find the optimal solution. the reduction of number of binary variables has resulted in a significant improvement with approximately 60% reduction in cplex iterations and 80% reduction in the time required for finding the optimal solution in phase we also compared the total time required for phase & to find the first solution by the proposed method with monolithic milp and graph theory based approach without grouping. the results are given in table  the time required by the proposed method is ~ 4% and 22% of the time required for the monolithic milp and graph theory based approach, accordingly. for all solutions, the proposed approach required 16311 seconds. the large computational time required for monolithic milp restrained its use for finding all minimal reaction sets. nonetheless, to validate the results monolithic milp was employed after excluding all the minimal reactions from search space. it found a sub-optimal solution with reactions in the minimal set. this guarantees that the proposed method identifies all minimal reaction sets. discussion and  development of cells with minimal metabolic functionality is increasingly gaining importance. the presence of redundant reactions in metabolic networks results in multiple minimal reactions sets that can meet the predefined cellular objectives. in this paper, we proposed a graph theory augmented recursive milp approach to identify all the minimal reaction sets in a metabolic network. the proposed method has been demonstrated by finding all the minimal reaction sets for escherichia coli and saccharomyces cerevisiae. the proposed approach correctly identified all the minimal reaction sets in both the cases. we also proposed the concept of grouping dependent reactions to reduce the number of binary variables for milp formulation. in the present study, several groups of dependent reactions are identified in escherichia coli and saccharomyces cerevisiae and exploited to reduce the number of binary variables and consequently the solution time. since the use of binary variables is very common in metabolic network analysis for identifying strain improvement strategies, the reaction group concept will benefit the other applications as well. here, we have developed a graph based algorithm that exploits the structure of the metabolic network to identify groups of dependent reactions. we now compare the groups of dependent reactions identified by the proposed method with the previously reported approach based on steady-state flux distribution. we used the metatool software to find the dependent reaction groups using steady-state flux distribution. while the proposed graph based approach found reaction groups with norm more than in the yeast model used in case study the flux based approach identified dependent reaction groups. the norm range of the reaction groups identified by the proposed approach is while that for the flux based approach is. the lower range of the proposed approach could be due to the strict use of the structure of the metabolic network in contrast to the flux based approach which can place structurally unconnected reactions in a group. nonetheless, the use of structure of metabolic network identified some dependent reaction groups that were not identified by the flux based approach due to the use of reaction directionality as described in reaction dependency and grouping. while there is a large overlap between the dependent reaction groups found by the two approaches, several groups are uniquely identified by each approach. hence, combining the two approaches could further simplify the metabolic network. generally, different minimal reaction sets have different structural and functional properties in terms of robustness, predictable metabolic interactions, practically achievable metabolic fluxes, and thermodynamically non-favourable pathways and cycles. we identified one such cycle in one of the solutions identified in e.coli model in case study and very high flux through co-factor recycling reaction in another solution. thus, of the three, only one solution has practically achievable fluxes without any coupled reactions and is suitable for developing the minimal metabolism cell. this clearly shows the importance of identifying all minimal reaction sets. in case study the proposed method identified solutions. while some of them use glycolysis pathway to produce dhap others use pentose phosphate pathway. other differences observed include use of different reactions for the consumption of dhap and formate, recycling of nad, nadh, nadp, and nadph, and production of ammonium. a preliminary analysis shows that some minimal reaction sets require complex set of reactions with impractical high fluxes compared to others. further analysis may reveal the details of these pathways that qualify them to developing into minimal metabolism cells. competing interests the authors declare that they have no competing interests. authors’ contributions both sj and rs contributed to the concept and methodology development. sj implemented the methodology and conducted the analysis of results and biological interpretation. rs supervised the study and assisted in implementation. sj drafted the manuscript. both authors read and approved the final manuscript. 
graph theoretical insights offer systematic methods to not only reduce the number of variables in math programming and increase its computational efficiency, but also provide efficient ways to find multiple optimal solutions. existing computational methods to identify minimal reaction sets in metabolic networks are computationally expensive. further, they identify only one of the several possible minimal reaction sets. in this paper, we propose an efficient graph theory based recursive optimization approach to identify all minimal reaction sets. in case study the proposed method identified three minimal reaction sets each containing reactions in escherichia coli central metabolic network with reactions. analysis of these three minimal reaction sets revealed that one of them is more suitable for developing minimal metabolism cell compared to other two due to practically achievable internal flux distribution. in case study the proposed method identified minimal reaction sets from the saccharomyces cerevisiae genome scale metabolic network with reactions. the proposed method required only  hours to identify all the minimal reaction sets and has shown a significant reduction in the solution time when compared to the existing methods for finding minimal reaction set. identification of all minimal reactions sets in metabolic networks is essential since different minimal reaction sets have different properties that effect the bioprocess development. the proposed method correctly identified all minimal reaction sets in a both the case studies. the proposed method is computationally efficient compared to other methods for finding minimal reaction sets and useful to employ with genome-scale metabolic networks. development of cells with minimal metabolic functionality is gaining importance due to their efficiency in producing chemicals and fuels. the efficacy of the proposed approach is demonstrated using case studies from escherichia coli and saccharomyces cerevisiae.
 development of cells with minimal metabolic functionality is gaining importance due to their efficiency in producing chemicals and fuels. existing computational methods to identify minimal reaction sets in metabolic networks are computationally expensive. further, they identify only one of the several possible minimal reaction sets.  in this paper, we propose an efficient graph theory based recursive optimization approach to identify all minimal reaction sets. graph theoretical insights offer systematic methods to not only reduce the number of variables in math programming and increase its computational efficiency, but also provide efficient ways to find multiple optimal solutions. the efficacy of the proposed approach is demonstrated using case studies from escherichia coli and saccharomyces cerevisiae. in case study the proposed method identified three minimal reaction sets each containing reactions in escherichia coli central metabolic network with reactions. analysis of these three minimal reaction sets revealed that one of them is more suitable for developing minimal metabolism cell compared to other two due to practically achievable internal flux distribution. in case study the proposed method identified minimal reaction sets from the saccharomyces cerevisiae genome scale metabolic network with reactions. the proposed method required only  hours to identify all the minimal reaction sets and has shown a significant reduction in the solution time when compared to the existing methods for finding minimal reaction set.  identification of all minimal reactions sets in metabolic networks is essential since different minimal reaction sets have different properties that effect the bioprocess development. the proposed method correctly identified all minimal reaction sets in a both the case studies. the proposed method is computationally efficient compared to other methods for finding minimal reaction sets and useful to employ with genome-scale metabolic networks. systems biotechnologystrain developmentminimal cellmixed-integer linear program multiple solutions 