BACKGROUND
metabolism is the process of synthesis and degradation of molecules occurring in living organisms.
metabolism is generally represented as a network where metabolites are interconnected by reactions.
in order to give a functional description of metabolism, metabolic networks are often decomposed into separated parts, called metabolic pathways.
the description of metabolism through metabolic pathways is useful, even though any division in pathways is arbitrary, because it helps in modeling and understanding the behavior of the full network.
a metabolic pathway can be defined as a coherent set of enzyme-catalyzed biochemical reactions by which a living organism transforms a set of source compounds into a set of target compounds.
by regulating enzyme and protein synthesis, living organisms can adapt to different environments.
this model of metabolism as composed by independent metabolic pathways is simplistic, since pathways are nested and interdependent.
in fact, metabolism is a complex system and pathways interact with each other.
the aim of the work presented here is to find all the viable sets of heterologous enzymes, which can produce a predefined target compound when added to the pool of endogenous enzymes of a given organism.
our method enables a metabolic engineer to find all heterologous metabolic pathways producing a target compound, for instance liquiritigenin , from the endogenous metabolites of e. coli k- <dig>  as shown in figure  <dig>
as depicted in figure  <dig>  our problem can be formulated as searching for all possible heterologous pathways linking a target compound to the endogenous metabolites of an organism.
to this purpose we provide software tools that enable the discovery of potential pathways producing a target chosen by the user  <cit> .
more precisely the user enters a target compound, a chassis organism, and our software tools return a ranked list of pathways  to be engineered into the chassis organism.
to achieve this task we have developed an approach composed of three steps.
in the first step, using a retrosynthesis software, reactions producing a target compound are iteratively searched backwards until the set of needed precursors only contains source metabolites.
this first step returns a retrosynthetic network connecting a target compound to the endogenous metabolites of an organism.
there may be several pathways in the retrosynthetic network linking the source metabolites to the target compounds and there is thus a need to enumerate all the possibilities.
pathway enumeration is performed by in the second step.
once the pathways have been enumerated, we evaluate in the third step the possibility to insert each pathway and its associated heterologous enzymes in the host organism.
this step consist of determining the catalytic efficacy of the enzymes, the toxicity of the products and the coproducts  <cit> , and the easiness of inserting the enzymes into the host.
the efficiency of the pathways can then be further estimated by flux models for the cell metabolism such as flux balance analysis  <cit> .
we have already discuss elsewhere the first and third step  <cit> , i.e.
methods to generate retrosynthetic networks and methods to rank pathway efficiency.
to apply these methods in the context of heterologous target production, we need a computationally fast method to enumerate all possible pathways.
we address the enumeration problem in the current paper.
different mathematical models that describe metabolism have been proposed .
we distinguish two main families of approaches: the ones computing steady states of the fluxes of reactions  and the ones based only the topology of the network.
typically, steady states are studied and simulated by generating the flux space.
of particular interest are the extreme pathways and the elementary modes, they both represent the smallest  generating set of the flux space and they both are composed of independent non-decomposable pathways in the network  <cit> .
the differences between extreme pathways and elementary modes have already been discussed in details  <cit>  and these differences arise when dealing with reversible reactions.
in the present paper we consider all reaction irreversible, and networks comprising reversible reactions are modeled by doubling each reversible reaction into a forward reaction and a reverse reaction.
algorithms have been developed to enumerate both extreme pathways  <cit>  and elementary modes  <cit>  and these algorithm are all variants of the double description method  <cit> , which enumerates all extreme rays of a polyhedral cone.
the algorithms use as input a stoichiometric matrix  representing the network  and output sets of fluxes  satisfying sv =  <dig>  one notices that extreme pathways and elementary modes while representing pathways  do not directly enumerate all pathways linking a source set to a target set of compounds.
however as shown in the subsection "enumerating pathways using the steady state approach" one can construct stoichiometric matrices where input fluxes are added to the set of source compounds and outgoing fluxes are associated to the target and heterologous coproducts such that the extreme pathways and elementary modes enumerated from these matrices do correspond to all pathways linking the source set to the target.
while as mentioned above, the problem of systematically enumerating pathways for heterologous production in chassis organisms has not yet been addressed, there are methods based on the steady state approach to search for heterologous pathways optimizing target productions  <cit> , and methods to search for shortest pathways between source and target sets of compounds  <cit>  and  <cit> .
all these methods are based on optimization and make use of integer linear programming.
precisely, the method of pharkya et al.
<cit> , is aimed at redesigning microbial chassis organisms through heterologous reaction addition and native reaction deletion for the overproduction of a target compound.
the addition and deletion are parameterized using binary variables attached to each reaction.
a mixed integer linear program  is then set up to maximize the target yield while minimizing the number of added reactions.
the methods of de figueiredo et al.
<cit> , and pey et al.
<cit>  are both aimed at searching for the k shortest pathways.
in de figueiredo et al.
<cit>  the k first shortest pathways are searched in entire metabolic networks, while in pey et al.
<cit>  the pathways are searched between a source metabolite and a target metabolite.
both methods solve the problem at steady state and search for fluxes, v, verifying sv =  <dig>  while minimizing the number of reactions turned on .
aside from the fact that integer linear programs suffer from computational complexity  all the above methods search for at most k optimized  pathways and do not guarantee a full enumeration of the possibilities.
in our methods the optimal pathways are computed in a post process by ranking the pathways that have been enumerated.
our approach allows one to decouple enumeration from optimization, and thus to plug any optimization criteria, including nonlinear functions and not only target yield or pathway length .
aside from using extreme pathways and elementary modes, we also present in this paper a topological model which directly enumerates all the possible heterologous pathways linking target compounds to a source set of compounds.
the main advantage of the topological approach compared to the stationary state approach is computational speed.
speed is in fact an important aspect when searching for the best pathways to engineer, as there are generally a combinatorial number of pathways between given source sets and target sets.
as an illustration of this combinatorial complexity, the work of hatzimanikatis et al.
<cit> , which provides a list of  <dig>  novel biochemical routes from chorismate to phenylalanine, and the work of cho et al.
<cit> , which enumerates  <dig>  reaction routes to produce isobutanol.
there exist standard graph-based methods to search and eventually enumerate pathways in metabolic networks, but these methods including pathfinding  <cit>  and pathway hunter tool  <cit>  are computing pathways and shortest pathways in graphs instead of hypergraphs.
the particularity of these techniques is that only main substrates and main products are taken into account when constructing pathways, and consequently these main compounds must be differentiated from the cofactors .
in the work of croes et al.
<cit>  cofactors are filtered out based on their connectivity in the network.
indeed, compounds highly connected such as atp, nadp, or h2o are cofactors of most reactions as they do not share carbon atoms with the products of the reactions.
in a more recent work  <cit> , the main compounds in the pathways linking source metabolites to target metabolites are detected using the kegg rpair annotation  <cit> , which enables one to follow the fate of atoms when going from a set of substrates to a set of products.
another approach to search for main substrates and main product is the one developed with the pathway hunter tool, which consists of mapping substrates to products using cheminformatics fingerprints.
while all the above techniques are computationally efficient, their main shortcoming is that they are not able to encompass reactions when a main product is formed from two main substrates.
there are plenty of such reactions in metabolic networks, consider for instance the formation of guanidinoacetate from arginine and glycine through a glycine amidinotransferase , or the formation of glutathione from γ-l-glutamyl-l-cysteine and glycine catalyzed by a glutathione synthase .
recently, some of the above topological methods have been benchmarked against the integer linear programming technique mentioned above  <cit>  to search for the shortest pathways linking various compounds, the recovery ratio for a set of  <dig> predefined reference pathways could not reach 100% with the graph based approach, exemplifying the shortcoming of that approach.
as reviewed above, while there are methods and theoretical results to enumerate elementary modes or extreme pathways and graph based techniques to search for pathways in a given metabolic network, to the best of our knowledge there is no known methods to directly enumerate pathways in the context of metabolic engineering, that is, to enumerate all the pathways encompassing all substrates and products necessary and sufficient to produce a given set of target compounds from a given set of source compounds.
in the present paper we address that specific problem and present two methods one based on elementary modes  and one based on a direct enumeration algorithm .
in order to address this problem we need, in addition, to consider the problem of determining supplement molecules, i.e.
metabolites that the organism cannot synthesize, but which can be added to the growth media in order to increase the number of viable pathways; and bootstrap molecules, i.e.
metabolites which are required fist in order to be produced  <cit> .
while in the general context of metabolic network analysis, finding elementary modes does not require to first search for bootstrap molecules, in the context of metabolic engineering however any heterologous pathway solution that comprises a compound that is first consumed before being produced is valid only when the compound is added to the growth medium.
therefore, in our study in the context of metabolic engineering, there is a need to first compute the bootstrap compounds prior to elementary modes.
the paper is divided as follows.
in the methods section we first provide some definitions, then outline our algorithms to solve the pathway enumeration problem with both the steady state approach and the topological approach.
the problem of finding and enumerating all the pathways going from a large source  to a target chosen by the user is considered.
all the algorithms presented for the topological approach  have polynomial worst-case running time, the algorithm for enumeration is a polynomial time per output algorithm on some classes of hypergraphs.
we also provide algorithms to determine supplements, which are metabolites that the organism cannot synthesize, but which can be added to the growth media in order to increase the number of viable pathways.
furthermore, an analysis of pathways containing supplements allows finding out pathways that contain bootstrap molecules.
in the results and discussion section we illustrate our algorithms with the enumeration of the possible pathways to synthesize more than  <dig> compounds in e. coli.
we "experimentally" probe the computational complexity of the steady state and topological approaches for a series of networks of growing sizes and discuss the theoretical complexity results of the topological approach, which are provided in appendices a and b.
while we illustrate our two methods for the production of heterologous compounds using as a source set all the endogenous metabolites of e. coli, our methods can be applied to any chassis organism and more generally to any source set .
methods
in the context of metabolic engineering, metabolic networks have been represented as directed graphs .
in such graphs, edges are directed and correspond to reactions connecting two compounds if one is the product of the other.
directed graphs can represent monomolecular reactions , but they are not well suited to capture more complex reactions.
as already discussed in the background section, when representing bimolecular reactions, one has to choose which molecules are connected by the edges of the graph and which ones have to be excluded from the graph because they are considered co-substrates or co-products.
additionally, one of the limitations of a model based on a graph representation is that depending on the criteria used to identify the co-substrates and co-products in the reactions, the networks obtained are different.
in the present paper to palliate the limitations of the directed graph model, we represent networks as directed hypergraphs.
the first examples of modeling through hypergraphs can be found in  <cit> .
in a hypergraph, each hyperarc connects a set of vertices, corresponding to reactants, to a disjoint set of vertices, representing the products.
in our model each hyperarc corresponds to a reaction that can be catalyzed by an enzyme.
it is worth noticing that hypergraph models have already been used to find minimal sets of metabolites sufficient to produce a set of target metabolites  <cit> .
unfortunately, the algorithms proposed in  <cit> , do not enumerate pathways and are therefore not directly applicable to our metabolic engineering problem.
definitions
definition  <dig> .
a directed hypergraph is a pair h= where v = {v <dig>  v <dig> ..., vn} is the set of vertices and e = {e <dig>  e <dig> ..., em} is the set of hyperarcs.
a hyperarc ei is an ordered pair ei =  of disjoint subsets of v.
the set xi is also called the tail of ei and the set yi is called the head, with reference to the graphical representation of arcs  and hyperarcs as arrows.
we denote by x:e→p the application that given an hyperarc ei returns its tail x ⊂ v. analogously we use y:e→p for the application that given a hyperarc returns its head.
definition  <dig> .
in a metabolic network each vertex corresponds to a metabolite and each hyperarc corresponds to a reaction.
a metabolic network of m metabolites and n reactions can be represented with a m × n stoichiometric matrix s, where the rows correspond to the m metabolites and the n columns to the reactions.
a reaction j is represented by the column vector sj = t where sij is the stoichiometric coefficient of metabolite i in reaction j. reactants have negative coefficients and products have positive coefficients.
examples of hypergraph, network, and stochiometric matrix are given in figure 2a.
we notice that the stoichiometric coefficients of the reactions are not taken into account in the hypergraph representation.
we also notice that the pair  is ordered so to make the distinction between reactants and products.
in this representation reactions are irreversible.
many biochemical reactions can be considered as irreversible, since in organisms the homeostatic equilibrium is often strongly polarized.
nonetheless, metabolic network may comprise reversible reactions, and we model these reactions by introducing both hyperarcs: , and .
hyperpaths, a generalization of simple paths in graphs where cycle free paths going from one vertex to another, are used to represent pathways.
a hyperpath connects a source set of vertices to a target set of nodes.
two examples of hyperpaths are given in figures 2a and 2c.
we remark that in a natural way a set e of hyperarcs defines a hypergraph ε =  ∪ ∪e∈ey , e).
by abuse of the terminology we denote by e the hypergraph corresponding to the set e of hyperarcs and all the heads and tails of the hyperarcs in e. the following definition for hyperpaths is borrowed from nielsen et al.
<cit> .
definition  <dig> .
a hyperpath p going from a source subset sh of v to a target subset tp of p in a hypergraph h= is a hypergraph hp= with vp ⊆ v, ep ⊆ e, such that there is an ordering f of the hyperarcs ep with the following properties.
• ∀k∈{ <dig> …,|f|},x⊆sh∪)
• tp⊆sh∪)
from the point of view of metabolism, the first condition corresponds to the requirement that reactants of reactions participating in the hyperpath can be produced without the presence of the reaction itself.
hyperpaths defined in this manner represent a metabolic route from the source to the target.
according to definition  the hypergraph of figure 2b with source a is not a hyperpath because neither reaction r <dig> nor r <dig> can happen until the other does not start.
the definition , though complex, is computationally tractable, meaning that the time required to determine if a hypergraph is a hyperpath is proportional to the number of reactions.
a polynomial time algorithm to determine if a hypergraph is a hyperpath is given in  <cit> , the algorithm findall presented below can also be used for that purpose.
in fact, as discussed below, if the set of reactions returned by findall  contains all the reactions in hp, then hp is a hyperpath.
the metabolic network described by a hypergraph has to be as comprehensive as possible, containing every known enzyme-catalyzed reaction occurring in organisms.
we say that a hyperpath produces a set of target metabolites if it contains all those target elements.
a set of target compounds is said to be reachable from a given source, or linked to the source, if there is at least one hyperpath producing the targets.
we are interested in the enumeration of pathways leading to the production of a desired compound.
hyperpaths do not generally give the best representation of pathways because hyperpaths can contain reactions not necessarily linking the target to the source.
minimal hyperpaths, cf.
definition , are an appropriate representation of pathways since they contain only the essential reactions linking the source to the target.
in the definition given below, we say that a hyperpath p is a subset of another hyperpath p′ if v ⊆ v' and e ⊆ e'.
for instance the hyperpath of figure 2c is a subset of the one of figure 2a.
definition  <dig> .
a hyperpath  with target tp is said to be minimal if it has no proper subsets with the same target.
the target is disconnected from the source if a reaction is removed from a minimal hyperpath.
in this sense minimal hyperpaths cannot be reduced.
from a metabolic engineering perspective the concept of minimal hyperpath is useful as it defines the minimum set of reactions necessary to produce a target heterologous compounds, and consequently the minimum set of enzymes needed to be inserted into the chassis organism where the compound is going to be produced.
in the following we define b to be the set of all molecules linked to the source for a given hypergraph h and source set sh.
the characterization of b is the first task to be solved before the enumeration.
once this set is known all the minimal hyperpaths can be enumerated for all the molecules associated to the vertices in b.
supplements
supplements for a target are molecules whose presence in the source set increases the number of pathways for target production.
finding supplements is an important improvement when exploring ways to produce the target, since they make possible new pathways.
for each target of interest one can look for vertices that once inserted in sh give place to pathways otherwise impassable.
in terms of metabolism we are looking for the "supplement" molecules, i.e., molecules that once introduced in the source set permit to find more pathways than those otherwise available.
we introduce below findsupp, an algorithm that returns the supplements.
an analysis of pathways containing supplements allows to find out pathways containing bootstrap molecules, i.e.
metabolites that are needed in reactions producing compounds afterwards used for the production of the bootstrap molecules.
as a matter of fact, many pathways can be made viable once bootstrap molecules become available in the metabolic network .
loosely speaking bootstrap molecules are molecules that cannot be produced by the reactions belonging to a hyperpath unless they are already present in the source.
cottret et al  <cit>  stated that given a source set the existence of a pathway making use of bootstrap molecules can be tested in polynomial time.
we provide later in this section an algorithm returning the bootstrap compounds, such algorithm can be used to determine if a target molecule is connected to the source through a pathway making use of bootstraps.
enumerating pathways using the steady state approach
in steady state, all possible pathways in a metabolic network are by definition stoichiometrically balanced, i.e.
all metabolites produced from the source set must be consumed except for those that are target products.
extreme pathways and elementary modes are two methods that compute the set of independent non-decomposable pathways in the network that generate all feasible steady state solutions in the flux space.
they do not directly enumerate all pathways linking a source set to a target set of compounds.
however, one can construct stoichiometric matrices where input fluxes are added to the set of source compounds and outgoing fluxes are associated to the target and heterologous co-products such that the extreme pathways and elementary modes enumerated from these matrices can be used to generate all pathways linking the source set to the target.
given a hyperpath hp= of a hypergraph h=, we can define a set of flux vectors vp for the hyperpath where components vpj corresponding to those reactions in the pathway ej∈hp are activated:
vpj= >0ej∈hp0ej∈h\hp
a hyperpath hp= of a hypergraph h= with input source subset sh and the target subset tp is defined as stoichiometrically balanced if the rows corresponding to each metabolite vi ∈ v that are obtained from the product of the stoichiometric matrix s and the associated flux vector vp verify:
svp=≤0vi∈sh≥0vi∈tp0vi∈v\sh,tp
a way to introduce the constraint on input and output metabolites in the previous equation is by adding to the stoichiometric matrix s additional columns corresponding to input reactions , and output reactions .
these auxiliary reactions, even if non-properly balanced in terms of the law of conservation of mass, are useful in order to define completely the problem in a compact manner:
sv=0v≥0v∈r
both extreme pathways and elementary modes make use of this formulation in order to compute the set of feasible solutions v. since in our hypergraph definition all reactions are irreversible, the set of pathways solving equation  <dig> computed by both extreme pathways and elementary modes are identical .
furthermore, solutions in v must contain only positive or null fluxes.
in order to determine all stoichiometrically balanced heterologous pathways hp that can be inserted into the chassis organism to produce a target set tp, we need to constrain the computation of elementary modes only to those that have non-zero fluxes for heterologous reactions.
efficient solutions to this problem have been considered in the divide-and-conquer approach  <cit>  by rearranging the constraints in an echelon form so that the constraints containing only the desired reactions appear at the bottom.
to define the constraints in our case, we consider first the hypergraph rt that is formed only by heterologous reactions.
this hypergraph rt is the subset of the hypergraph r formed by those hyperedges where at least one vertex v does not belong to the source set sr, i.e.
those metabolites endogenous to the chassis organisms.
by considering rt instead of the full hypergraph r, we are looking only for biosynthetic pathways involving heterologous reactions and therefore avoiding cycles internal to the chassis organism.
therefore, to compute all feasible steady state heterologous pathways, we reformulate equation  <dig> so that the stoichiometric matrix s is defined by reactions in rt; the input is given by all substrates in the source set sr∩x; and the output by all products of the reactions in the hypergraph y.
finally, from the computed set of solutions v for equation  <dig>  we are interested in enumerating all minimal hyperpaths from sr to the target set t on the hypergraph given by rt.
according to definition  <dig>  minimal hyperpaths for some target t are given by those cycle-free solutions in v containing only reactions linking the source to the target.
since any feasible flux pattern v is a superposition of elementary modes with non-negative coefficients  <cit> , the set of minimal hyperpaths for a given target t is a subset of the elementary modes producing t that are solution of equation  <dig>  namely, any feasible solution generated from the elementary modes will contain at least as many reactions as the ones that are in those elementary modes that form its basis.
therefore no additional minimal hyperpaths can be generated in this case by superposition of elementary modes.
enumerating pathways using the topological approach
the algorithm findall that allows to find b, the set of metabolites that can be linked to the source sh by a hyperpath.
findall, by explicitly constructing the ordered set f in definition , provides a proof of the tractability of the problem of checking if a hypergraph is a hyperpath.
moreover findallf permits to prune the original hypergraph enabling a faster enumeration algorithm.
as presented below the algorithm minimize, when called on the output of findall, returns, if exists, a minimal hyperpath linking a given target to the source.
these algorithms are the main components of the algorithm enumerating the pathways findpath described next.
then we present findsupp an algorithm to enumerate supplements.
finding one minimal hyperpath
let h= be the hypergraph representing the set of metabolic reactions, n = |v|, m = |e| and let sh be the set of source vertices representing the source metabolites.
the algorithm findall returns all the reactions that can contribute to the production of any element in b, i.e., the set of all compounds that can be connected to the source.
findall is a linear algorithm in the number of vertices, hyperarcs and in the total coordination; the complexity is o |+ |y-1|) that is bounded by o. therefore, such algorithm can be applied to the hypergraph h of all reactions in order to obtain a pruned sub-hypergraph h′= where the set of vertices v′:=sh∪b, and the set of edges e' is the set of reactions returned by findall.
in the context of metabolic engineering findall returns all the compounds that can be produced from a given set of source compounds and reactions.
for instance, using findall with all know metabolic reactions one can determine all the compounds that can be produced from the metabolites of e. coli.
algorithm findall
input:
h,sh
<dig>  for all r in h
<dig>     x ← x
<dig>  end for
<dig>  v←sh
<dig>  d←sh
<dig>  f ← {∅}
<dig>  while v ≠ {∅}
<dig>     let i be an element of v
<dig>     v ← v \ i
<dig>     d ← d ∪ i
<dig>     for all r ∈ h such that i ∈ x:
<dig>        x ← x \ i
<dig>        if x = {∅}
<dig>           f ← {f, r}
<dig>           for all j in y  and not in d:
<dig>              v ← v ∪ j
<dig>           end for
<dig>        end if
<dig>     end for
<dig>  end while
output:
f
let d be the union of the source set and of the heads of all the reactions output in f by findall.
the correctness of the algorithm above is given by the following claims: every element in d is the target of some hyperpath or is part of the source, and every vertex in h that can be reached from the source is in d. for the first claim we can give a constructive proof by using the output vector f, the second claim is proved by contradiction.
• the proof of the fact that every element in d is reachable from the source is given constructively by the ordered set f returned by the algorithm.
in fact at each step f is a hyperpath.
this claim can be proved by induction on the steps of the algorithm, each time a hyperarc r is appended to f  the tail x is contained in d  and if a vertex j is added to d  it means that it was previously added to v  and thus it was in the head y  of some hyperarc already in the hyperpath.
• the second claim can be proved by contradiction: if an element of b were not in d there would be a hyperpath linking it to the source.
in such hyperpath let consider the first  reaction r whose x belongs to d and such that one of the elements of y  does not.
consider among x the last one that has been inserted into the set v ; after its removal from x this set becomes empty and the elements of y  are inserted into v  and then in d , which is a contradiction.
from the above statements follows that each vertex appearing in a hyperpath having as source sh is an element of d and every hyperarc is an element of f. thus the algorithm findall provides an effective pruning of the original hypergraph: in h there is no minimal hypergraph with source sh containing hyperarcs not in ∪kfk or vertices not in b. the output hyperarcs are the only ones that can belong to a minimal hyperpath, and h′=),∪kfk) is the pruned hypergraph only containing reachable vertices and hyperarcs.
notice that findall algorithm as presented above returns in polynomial time a hyperpath valid for each target vertex in b. even though there are more efficient algorithms for finding a hyperpath for one single target, for the sake of simplicity we avoid to introduce here an additional algorithm and just remark that since findall is polynomial, the use of it does not affect the complexity analysis of the algorithms making use of its output.
remark that a minimal hyperpath going to a specific target can be easily extracted from the hyperpath output of findall.
namely, given a hyperpath p connecting s to t, it is always possible to find a minimal hyperpath p′ subset of p. moreover it can be done in polynomial time, for instance by using minimize , the algorithm introduced below.
minimize is an algorithm that takes as input a hypergraph p, a hyperpath rf subset of p, a target set of vertices t and a source s. if p does not link t to s the empty set is returned, otherwise a hyperpath contained in p, containing rf and linking t to s is returned.
in particular, if rf is empty, the output of minimize is the minimal hyperpath going from s to t, provided it exists.
minimize returns a hyperpath obtained by removing all inessential hyperarcs except for the ones in rf.
in the context of metabolic engineering, pathways containing a small number of heterologous reactions are generally preferred, since they are easier to engineer in the host organism.
therefore, given two pathways that produce the same target, where one is subset of the other, the one requiring the smaller number of heterologous reactions has to be selected.
this is the reason that makes relevant to obtain minimal hyperpaths from generic hyperpaths.
algorithm minimize
input:
p,rf,t,s
<dig>  f ← findall
<dig>  p' ← p
<dig>  if t ⊄ ∪k y
<dig>  p′←{∅}
<dig>  else
<dig>     for all r in p
<dig>        if r not in rf
<dig>           f ← findall
<dig>           if t ⊂ ∪k y
<dig>              p′←p′\r
<dig>           end if
<dig>        end if
<dig>     end for
<dig>  end if
output:
p′
the proof of correctness of this algorithm is simple and is based on the fact that p′⊆p implies findall⊆findall.
if a reaction in p has not been removed from p′, then any subset of p <dig> not containing r does not produce the target.
the worst-case time for this algorithm is o|+|y|)).
since x and y have bounded values, the algorithm has a quadratic complexity.
even though faster algorithms can be designed, here we presented this one because of its conceptual simplicity.
remark that, since minimize returns a minimal hyperpath p′ subset of p if it exists, then the minimality of a hyperpath p can be tested by checking whether p′=p or not.
a related problem to minimize is the minimal constrained hyperpath problem: the problem of finding if a minimal hyperpath from a given source to a given target, containing the hyperarcs in rf exists.
notice that minimize, although linked to this problem does not solve it.
in fact, if the output of minimize is an empty set then there are no minimal hyperpaths satisfying the constraints; however if the output is a minimal hyperpath then obviously a minimal constrained hyperpath exists; and finally, if the output is a non-minimal hyperpath then we do not know if a minimal hyperpath satisfying the constraints exists or not.
below we will discuss why we are interested in algorithms for the minimal constrained hyperpath problem, while in appendix a.
<dig> we show that in general the problem is np-complete .
pathways enumeration
the basic idea behind the enumeration algorithm presented below is to introduce an iterative refinement of partitions of the space of feasible solutions i.e.
of the space of hyperpaths and in each part to look for a solution.
in our implementation, a part is defined by two sets of reactions  of the original hypergraph.
these sets are used during the enumeration process, rf is a set of hyperarcs that must be present in the enumerated hyperpath and rn is the set of hyperarcs that must not be part of the enumerated hyperpath.
the problem of finding a solution in one of the parts is addressed at each iteration and if a solution is found the part is divided in finer parts.
this process is repeated until all the minimal hyperpaths have been found.
enumeration by means of the minimal constrained hyperpath problem
first we describe informally the enumeration algorithm through the toy example hypergraph in figure 2d and 3a, then we outline in figure 3b a typical run for a more involved example: liquiritigenin .
a minimal hyperpath p <dig> connecting the node v <dig> to the source nodes v <dig>  v <dig> on the hypergraph h of figure 2d can be obtained by calling minimize  on the hypergraph p′ obtained by findall.
the hypergraph p′ is represented in figure 3a.
once p1={r <dig> r3} has been obtained, the search space is divided into three parts:
• the hypergraphs which do not contain r <dig>
• the hypergraphs which do contain r <dig> and do not contain r <dig>
• the hypergraphs which do contain r <dig> and r <dig>
the first set does not contain hyperpaths connecting the target to the source: once the reaction r <dig> is removed, v <dig> is disconnected from the source.
the second set contains a solution and thus has to be partitioned.
the third set contains only one minimal pathway .
the minimal hyperpath in the second set is found by running findall on h\r <dig> and then minimize with constraint rf = {r4}.
the minimal hyperpath so obtained is the one only containing hyperarcs r <dig>  r <dig>  the set of the hypergraphs defined by  is partitioned in two parts defined by new sets of constraints.
the way the partition is done is explained in detail in algorithm findpath and gives two non overlapping sets:
• the hypergraphs which do contain r <dig> and do not contain r <dig> and r <dig>
• the hypergraphs which do contain r <dig> and r <dig> and do not contain r <dig>
the first of these sets does not contain hyperpaths going to v8: once r <dig> and r <dig> are removed, node v <dig> is disconnected from the source.
the second one only contains the second and last minimal hyperpath: the one consisting of hyperarcs r <dig>  r <dig>  the algorithm here sketched is based on the fact that all minimal hyperpaths are found once the problem of finding a minimal hyperpath has been solved for each part of the partition.
relaxed hyperpath minimization
the enumeration procedure is performed by the algorithm findpath, which enumerates all the minimal pathways and does not output duplicate hyperpaths.
precisely findpath  returns a set of hyperpaths containing all the minimal hyperpaths in h connecting t to sh and containing all the reactions in rf.
findpath  returns all the minimal hyperpaths from the source sh to the target t in h.
a schematic representation of how findpath works for the enumeration of the pathways of liquiritigenin is given in figure 3b where we represent each call with a box connected by an arrow to its parent process.
for each call of findpath either a new hyperpath is found and then findpath is executed with new constraints, or there are no new hyperpaths and the branching process is stopped.
the new constraints sets rf ', rn' for a new call of findpath are obtained by incrementing the sets rf, rn of the father process.
given an order for the hyperarcs of the hyperpath p found for the father process, the set rn' relative to the child process is constructed by incrementing rn by one element r belonging to p, the set rf ' is constructed by incrementing rf by all the hyperarcs coming before r. for each element in p not belonging to rn a child process is called.
findpath  returns all the minimal hyperpaths from the source sh to the target t in h. in the context of metabolic engineering findpath returns all the metabolic pathways for the production of the target compounds.
algorithm findpath
input:
h, rf, t, sh
<dig>  f ← findall
<dig>  p ← minimize
<dig>  en ← ∅
<dig>  if p≠∅
<dig>     en←p
<dig>        f←findall
<dig>        for all k in |f|,..., 1}
<dig>           r = fk
<dig>           if r not in rf :
<dig>              en←{en,findpath}
<dig>              rf ← rf ∪ r
<dig>           end if
<dig>        end for
<dig>     end if
output:
en
the loop at line  <dig> of findpath is done according to the order given by line  <dig> where the hyperarcs are ordered so that at least one of the head vertices of each hyperarc is a tail vertex of some previous reaction.
such an ordering is always possible since p is a hyperpath.
as said above and illustrated in figure 3b, findpath is an algorithm that iteratively calls itself, see line  <dig>  note that even if rn is not explicitly defined in findpath, it is constructed implicitly when at line  <dig> of findpath is called on the smaller graph h\r.
let us note that the output of the enumeration is not always composed of minimal hyperpaths.
this is due to the fact that the algorithm minimize while running in polynomial time can return a non-minimal hyperpath.
an algorithm always returning minimal hyperpaths cannot be polynomial since the problem of finding a minimal hyperpath containing a set rf of hyperarcs is an np-complete problem as showed in appendix a.
<dig>  however, in many practical instances , the algorithm minimize returns a minimal constrained hyperpath.
as a matter of fact, for all the enumeration studies we have so far carried out, we observed that the output obtained by minimize when called by the algorithm findpath introduced above was always a minimal hyperpath.
nonetheless, a characterization of hard instances of the minimal constrained hyperpath problem is given in appendix.
supplements enumeration
provided a given metabolic network and a set of source compounds  it may not be possible to link all the metabolites of the network to the source set.
when a target compound is not accessible from the source set, one can consider the possibility of inserting into the metabolism of the organism some precursors so that the target becomes reachable.
in practice such a task can be carried out through the enrichment of the growth media.
more generally, the insertion of supplements can be used even when the target compound is reachable in order to access to new pathways for the production of the target.
let a supplement for a target t be any compound i∉b that is involved as reactant in at least one minimal hyperpath going from a superset of sh∪i to the target t. below we give the algorithm findsupp finding the supplements for the production of a given target.
in figure 2d supplements are highlighted in red.
therefore, the process of finding supplements is use-ful as a general strategy in metabolic engineering in order to determine which metabolites might be part of the metabolism that produces a given target.
algorithm findsupp
input:
h, sh, t
<dig>  wishlist ← t
<dig>  d ← {∅}
<dig>  while wishlist \ d ≠ {∅}
<dig>     let i be an element of wishlist \ d
<dig>     d ← d ∪ i
<dig>     aux ← {∅}
<dig>     for all reactions r with i ∈ y
<dig>        aux←aux∪\)
<dig>     end for
<dig>        wishlist ← wishlist ∪ aux
<dig>  end while
<dig>  f←findall
<dig>  d←d\sh∪)
output:
d
bootstraps
bootstrap molecules relative to a source sh are the molecules that cannot be produced by a hyperpath with source sh unless they are already present in the media.
an example of bootstrap nodes are nodes v <dig>  v <dig> of figure 2d.
in this section we give an algorithm finding in polynomial time all the bootstraps of a hypergraph h with source vertices sh.
bootstraps are special kind of supplement, if at any step of a pathway, a heterologous metabolite is needed as a substrate but has not yet been produced from the source set, then this metabolite is a bootstrap and must be added in the growth media for the reaction to take place, and for the pathway to be a valid pathway.
the algorithm given below enables one to detect bootstraps prior enumerating pathways running the findpath algorithm.
algorithm findbootstraps
input:
h,sh
<dig>  f←findall
<dig>  d←sh∪)
<dig>  h′←{∅}
<dig>  for all r in h
<dig>     r' ←  \ d, y  \ d)
<dig>     if y  ≠ {∅}:
<dig>        h′←h′∪r′
<dig>     end if
<dig>  end for
<dig>  while exists v in ∪r∈h′y\∪r∈h′x
<dig>     for all r' containing v:
<dig>        r' ← , y \ v)
<dig>        if  = {∅}) or ):
<dig>           h′←h′\r′
<dig>        end if
<dig>     end for
<dig>  end while
<dig>  b=∪r∈h′y
output:
b
the findbootstraps algorithm is linear in the number of vertices, hyperarcs and in the total coordination.
remark that the set ∪r∈h′y obtained in line  <dig> is equal to ∪r∈h′x.
in fact the bootstrap vertices b constitute the largest set of vertices not reachable from the source and such that each element of the set belongs to the head of at least one reaction whose tail only contains vertices in b or in b. notice that the set of bootstrap vertices in a hypergraph h only depends on the source vertices and does not depend on the target.
RESULTS
to evaluate the performance of our topological approach , we have compared running times of this approach with the running times of the steady state approach presented in the methods section.
all our tests were run on a mac pro server  <dig> ×  <dig>  ghz quad-core intel xeon,  <dig> gb.
all the algorithms of the topological approach were implemented in python.
for the steady state approach we used two software products, one computing elementary modes and the other extreme pathways.
elementary modes were computed by using the matlab interface to the java implementation of efmtool version  <dig> .
<dig>  <cit> .
extreme pathways were computed by using the mac osx version of the expa program  <cit> .
the running time comparison test was performed for different random samples of the hypergraph h constructed from the kegg database  <cit> .
e. coli was chosen as source organism.
the hypergraph h was composed of  <dig> metabolites connected by  <dig> reactions including  <dig> metabolites endogenous to e. coli and  <dig> heterologous compounds.
each sampled hypergraph hs was built by randomly sampling a given fraction of the total reactions in the hypergraph h. tests for each sampling fraction were repeated  <dig> times.
in the case of elementary modes, we were able to run the test only up to 50% sampling of the full metabolic network, due to memory constraints in matlab.
for extreme pathways, the test was run up to 85% due again to memory constraints as well.
prior to enumerating pathways with elementary modes, extreme pathways or the direct algorithm, bootstraps were identified by using the algorithm findbootstraps defined in methods and added to the growth media.
the observed trend is consistent in the three cases, as it is shown in figure  <dig>  for the network composed of all  <dig> metabolites and  <dig> reactions, there are  <dig> pathways connecting  <dig> out of  <dig> heterologous metabolites to the source.
this relative low number of pathways producing heterologous compounds is related to the fact that heterologous maps usually show a degree of hierarchy higher than the one observed in native metabolic networks such as central metabolism .
it is indeed not surprising to find the same number of extreme pathways and elementary modes as all reactions in our networks are irreversible.
in the case of findpath, the fact that the topological approach does not take stoichiometry into account might produce some inconsistent pathways.
checking that a given pathway is stoichiometrically balanced can be done using linear programming  <cit> .
in the context of metabolic engineering  we found only few inconsistent pathways .
these pathways are eliminated after enumeration by our ranking function as one of the criteria to rank pathways is to solve the steady state equations  <cit> .
once unbalanced pathways were removed from the list of enumerated pathways, we obtained the same pathways as through the calculation of elementary modes.
these results are consistent with the given definitions of the algorithms.
in a forthcoming paper, we plan to generalize them into a formal proof.
we also find that findpath has the fastest execution time, the algorithm required less than  <dig> min to compute the full network.
elementary modes has an execution time approximately 10-fold slower than findpath.
the computation of extreme pathways, in turn, is the slowest one, being between  <dig> ~ 103-fold slower than findpath.
the execution times are on average of  <dig>  ±  <dig>  s per output with the findpath algorithm, of  <dig>  ±  <dig>  s with elementary modes, and of  <dig>  ±  <dig>  s with extreme pathways.
the observed execution times mean that for the compound with the highest number of pathways in our tests, the anti-diabetic drug acarbose containing  <dig> pathways, enumerating the pathways with elementary modes takes more than  <dig>  seconds, while it can be computed in  <dig>  seconds with findpath.
we have measured running times per output and memory usage as function of input size and output size .
in all cases we found linear growth o for both input size and output size.
in the case of running times, findpath and efmtool have running times per output approximately constant in function of the size of the input stoichiometric matrix s of  <dig>  × 10- <dig> and  <dig>  × 10- <dig> seconds, respectively, being expa the less efficient code with a constant time of  <dig>  × 10- <dig> seconds and a linear growth of  <dig>  × 10- <dig> seconds per input.
similar values were obtained depending on the size of the output, although the scaling in this case for expa was of  <dig>  × 10- <dig> seconds per output.
regarding memory usage, expa and findpath are less demanding, especially for output size, with linear growths of  <dig>  and  <dig> kb per output, respectively, while efmtool required significantly more memory allocation,  <dig>  ×  <dig> kb per output.
performance of the codes findpath, efmtool and expa in terms of run time per output and memory use in function of size of input stoichiometric matrix s and size of output .
to the best of our knowledge, the computational complexity of enumerating elementary modes on networks comprising irreversible reactions is up-to-date unknown  <cit> .
in appendix a.
<dig> we prove that enumerating minimal pathway is an np-complete problem  but as examplified in appendix b the hard instances are rather rare and obviously none of these hard instances were found in our running tests as all pathways returned by findpath were indeed minimal pathways.
to further probe the computational complexity of findpath we computed distribution of the run time for each of the  <dig> heterologous targets that are linked to e. coli.
as shown in figure  <dig>  the distribution is exponentially distributed, the average run time is thus finite.
examples of pathway enumeration for heterologous compounds with therapeutical applications were provided in our previous methodology study about pathway ranking  <cit>  for penicillin n and taxol .
these two examples contained a relative low level of combinatorial complexity, with a maximum number of  <dig> different pathways producing penicillin n. the combinatorial complexity of enumerating all putative biosynthetic routes producing a target compound, however, might become considerably higher than in these two previous examples .
further examples are those involved in the biosynthetic pathways for plant steroids leading to brassinolide, where  <dig> pathways can be used to produce campesterol, one of the initial precursors.
the number of pathways grows as we proceed downstream going up to  <dig> for teasterone, and to  <dig> for typhasterol.
finally, the number of pathways for the end products, castasterone and brassinolide is of  <dig>  this example illustrates how the complexity of pathways enumeration can grow with the number of intermediates involved in the synthesis of the final product.
since findpath algorithm has been designed with the aim to enumerate pathways in metabolic engineering applications, it does not directly address other general metabolic network analysis problems like finding shortest paths between two compounds of endogenous pathways in central metabolisms.
indeed, findpath does not enumerate pathways comprising compounds  that are consumed before being produced unless they were added as supplement in the growth medium .
another aspect of pathway design is the definition of a cost function that estimates limiting effects on production efficiency by multiple factors associated with genetic and metabolic engineering of the pathway, such as gene heterogeneity, metabolite toxicity, or steady-state fluxes.
in order to facilitate the designer in the selection of the best synthetic pathways to be implemented each hyperpath enumerated by our findpath algorithm is ranked depending on this cost estimation.
once a cost function to minimize is introduced, the search for an optimal pathway can be formulated as a shortest path problem.
the shortest path problem for weighted graphs consists of finding a path going from a given source vertex to a given target vertex while minimizing a cost function given by the sum of the costs of the arcs involved in the path.
in order to define a shortest-hyperpath problem we need a definition of cost for hyperpaths in a weighted hypergraph .
a natural generalization of the cost function for paths in graphs to hypergraphs is the sum of the costs of the hyperarcs contained in the hyperpath.
if on the one hand this generalization seems to be natural, on the other hand the two problems have, however, not the same complexity.
in fact the shortest-hyperpath problem with this cost function is known to be an np-hard problem .
the reason why the algorithms commonly used for graphs cannot be easily adapted for hypergraphs is that given the costs w for the hyperpaths ai⊂p leading to the vertices ai, tails of a hyperarc ej whose cost is wj, then the cost of the hyperpath containing all the hyperarcs in ∪iai and ej is not equal to wj+ ∑iw if the overlap of the hyperpaths ai is non trivial.
provided an additive cost function as defined in  <cit> , finding the shortest hyperpath can be done in polynomial time by an algorithm of complexity o) given in  <cit> , which finds the shortest hyperpath in b-hyperpaths.
in fact, this algorithm does apply to minimal hyperpaths given in definition .
generally speaking, the strategy we used in the enumeration algorithm  is often used to solve the problem of finding the first k solutions of a combinatorial problem, but it cannot be extended to the k-shortest problem on hypergraphs not even with an additive cost function , since when splitting the problem into a constrained part and a free one  we get a problem that is known to be np-hard .
therefore, the k-shortest path problem might be solvable efficiently if one is able to develop a specific cost function making the problem tractable.
however such an artifact cost function may not be necessarily appropriate to practical problems, such as the production of heterologous targets using metabolic engineering.
as an alternative we have proposed in this paper an enumeration algorithm, which systematically enumerates all pathways linking source compounds to target compounds.
in all the cases that we have so far processed for metabolic networks, our algorithm runs in polynomial time per output, and therefore pathways can be ranked by the designer based on their own user-defined cost functions.
CONCLUSIONS
in summary, the methods presented in this paper provide metabolic engineers with powerful tools that extend the toolbox for heterologous biosynthetic pathway design.
besides pathway enumeration of biosynthetic routes for a given target product, our methods have several other possible applications.
for instance, they can be used in combination with gene deletion strategies in order to determine pathway manipulations leading to overproduction of the target compound.
another application is in biodegradation and bioremediation, where our algorithms would need to be slightly modified in order to reverse the pathway search so that it can identify degradation routes for a given compound, while the underlying structure of the algorithms remains still valid.
finally beyond metabolism, our algorithms could also be utilized in the context of chemical synthesis to enumerate all the possible routes linking a target molecule to a source set of starting reactants, enabling the search for the best routes in terms of production costs.
availability and requirements
a web server is available: http://bioretrosynth.issb.genopole.fr/tools/metahype see details in additional file  <dig>
authors' contributions
jlf and pc designed the algorithms and the experiments.
df implemented the algorithms, and wrote the proof of the np-completeness for the minimal constrained hyperpath problem.
sbp wrote the example of pathway enumeration.
pc implemented the web service prototype, performed the experiments and analyzed the results.
all authors wrote the manuscript.
all authors read and approved the final manuscript.
appendix a reduction proofs for the shortest hyperpath problem and the minimal constrained hyperpath problem
a.
<dig> shortest hyperpath problem
in  <cit>  a reduction of the shortest-hyperpath problem to minimum set cover  is given.
we have to adapt the proof to our case for two reasons: the definition of directed hypergraph that was used is more restrictive | = 1) and the hyperpath was ill-defined.
in fact the given definition by these authors does not permit to say if some hypergraphs  are also hyperpaths or not.
in other words their definition is ambiguous: does not permit to determine the nature of all the directed hypergraphs and thus can be completed in several ways.
nonetheless the proof of np-hardness they gave is valid for our more general hypergraphs and minimal hyperpaths because the set of directed hypergraphs employed in the reduction proof in  is a sub-ensemble of the directed hypergraphs we defined above in the definitions section and all the hypergraphs employed for the reduction are well defined as hyperpaths, independently of the way the incomplete definition they gave is completed.
since our definition is a way to complete the definition in , then the two definitions agree on the set of hypergraphs employed for the reduction.
from these facts follows that the shortest hyperpath problem is an np-hard problem.
and, in particular, if the weights on hyperarc are non-negative, since hyperpaths always contain at least one minimal hyperpath, the shortest minimal hyperpath problem is np-hard too.
additive cost functions
the reason why the algorithms commonly used for the shortest path problem on graphs cannot be easily adapted for hypergraphs is that given the costs w for the hyperpaths ai⊂p leading to the vertices ai, tails of a hyperarc ej whose cost is wj, then the cost of the hyperpath containing all the hyperarcs in ∪iai and ej is not equal to wj+ ∑iw if the overlap of the hyperpaths ai is non trivial.
in order to define a shortest path problem that can be solved polynomially by a variant of dijkstra algorithm, the additive cost functions have been introduced in  <cit>  for the b-hyperpaths.
we adapt below the notion of "additive" cost function for hyperpaths.
a cost function w returning the cost for reaching the vertex x with the hyperpath p starting from a source s whose elements s ∈ s have w:=  <dig> is additive if w  is the minimum over all the arcs ex∈p whose head contains x of ex + f where wi:=w are the costs for reaching the tail vertices i of ex in p, f is an increasing monotone function of its argument and f ≥ wi ∀i.
remark that the cost of a hyperpath determined with an additive cost function in general is not given by the sum of the costs of the hyperarcs.
a.
<dig> minimal constrained hyperpath problem
consider a 3-sat instance concerning n variables σi and m clauses xj consisting of the problem of deciding if there exists an assignment of boolean values to the σi such that all the clauses are satisfied.
for each boolean variable σi contained in at least one clause introduce one hyperarc εi with the head of each εi having two vertices y = {vi+, vi-}.
for each clause xj consider a vertex νj and seven hyperarcs .
a boolean assignment is a triple a <dig>  a <dig>  a <dig> of boolean values.
let these hyperarcs be μj <dig> ..., μj <dig> and let the head y of a hyperarc μjk corresponding to the combination a <dig>  a <dig>  a <dig> of the boolean variables σj <dig>  σj <dig>  σj <dig> be y = {vj1a <dig>  vj2a <dig>  vj3a <dig>  νj}.
now let the tails of each hyperarc introduced be connected to the source nodes.
and let consider a node t being the product of the reaction r having as substrates the vertices νj and the heads {vi+, vi-} of the hyperarcs εi.
given the hypergraph described above  consider the minimal constrained hyperpath problem where all the hyperarcs εi are mandatory, and the target is t.
a solution of this problem gives in linear time a solution for the underlying sat problem, which makes the problem of minimal constrained hyperpath an np-complete problem.
in fact, given a minimal hyperpath m, solution of this problem, for each i consider vi* the only one of the two head vertices {vi+, vi-} belonging to the head of one or more of the μ arc in m .
the boolean assignments σi=vi* are a solution of the 3-sat problem and, inversely if a solution of the 3-sat problem exists then a minimal pathway satisfying the constraints exists and is the one obtained using only one of the hyperarcs for each xj among the ones whose head only contains νj and vertices vi*.
in additional file  <dig>  figure s <dig> for simplicity we consider the reduction of a single-clause satisfaction problem to finding if a minimal hyperpath satisfying the constraints exists.
there exist seven minimal hyperpaths connecting the target vertices to the source and satisfying the constraint that ε <dig>  ε <dig>  ε <dig> are parts of the hyperpath.
each solution corresponds to a valid boolean assignment of the variables σ <dig>  σ <dig>  σ <dig>
appendix b hard instances of minimal constrained hyperpath problem
on many hypergraphs the algorithm enumerating the pathways only returns minimal hyperpaths, this is the case for the metabolic networks that we analyzed in the main sections of this paper.
in this section we give a characterization of the hypergraphs where the algorithm minimize solves the minimal constrained hyperpath problem, characterizing these instances helps to individuate which hypergraphs are expected to give an output only containing minimal hyperpaths.
let y be the set of all the metabolites produced by reactions in rf, the mandatory reactions: y:=∪r∈rfy.
given a hypergraph h and the sets rf, rn we say that the well-separation condition holds if for every reaction r∈h\ the set y of products of r either is a subset of y or does not contain elements of y. if the well-separation condition holds for a hypergraph h with constrained reactions rf, the algorithm minimize returns a minimal hyperpath solving the minimal constrained hyperpath problem if a solution exists.
the well-separation condition holds for every choice of rf in a hypergraph whose reactions have one only product, as the hypergraphs defined in  <cit> .
if on the one hand, this condition can appear too constraining, on the other hand it can be generalized to larger sets of hypergraphs.
for instance, the algorithm minimize returns a minimal hyperpath solving the minimal constrained hyperpath problem even if the well-separation condition holds on the pruned graphs obtained by keeping from the original hypergraph only the reactions belonging at least to one hyperpath linking the target to the source and only the nodes being tail of such reactions.
examples of hard instances can be found among the ones used for the proof of np-completeness.
in general, hard instances h of the enumeration problem have to violate the condition of well-separation for same choice of rf⊂h, in order that the corresponding minimal constrained hyperpath problem becomes hard.
this happens if several compounds are products of more than one reaction producing more than one compound.
this is the case for nested networks as the one in additional file  <dig>  figure s <dig>  while the given network is small enough to be solvable by hand, it contains nevertheless the principal ingredients of complexity that would asymptotically make harder the problem as the size of the instances grows.
finding one minimal hyperpaths leading to the production of v <dig> is a simple problem, but finding new ones gets more and more involved.
this is a consequence of the fact that the nodes v <dig>  v <dig>  v <dig>  v <dig>  v <dig> can be produced by different choices of the reactions r <dig>  r <dig>  r <dig>  r <dig>  r <dig> and each of these reaction has more of one product susceptible to participate to the production of the target.
supplementary material
additional file 1
figure s <dig>  distribution of graph hierarchies  in heterologous metabolic networks  in comparison with the graph hierarchy of central , nucleotide , lipid , and amino acid  metabolic networks in e. coli.
click here for file
additional file 2
figure s <dig>  performance comparisons for findpath, expa and efmtools for run time per output and memory use in function of size of input and size of output.
click here for file
additional file 3
an example of pathway enumeration with bootstraps.
an example of pathways enumeration with findpath using bootstraps.
click here for file
additional file 4
availability and requirements.
description of the available metahype web server for running the algorithms for kegg compounds.
click here for file
additional file 5
figure s <dig>  minimal constrained hyperpath problem reduction of a 3-sat formula.
click here for file
additional file 6
figure s <dig>  instance of minimal constrained hyperpath problem.
click here for file
acknowledgements
funding: genopole® through an atige grant; anr through a chair of excellence.

 metabolism is the process of synthesis and degradation of molecules occurring in living organisms. metabolism is generally represented as a network where metabolites are interconnected by reactions. in order to give a functional description of metabolism, metabolic networks are often decomposed into separated parts, called metabolic pathways. the description of metabolism through metabolic pathways is useful, even though any division in pathways is arbitrary, because it helps in modeling and understanding the behavior of the full network. a metabolic pathway can be defined as a coherent set of enzyme-catalyzed biochemical reactions by which a living organism transforms a set of source compounds into a set of target compounds. by regulating enzyme and protein synthesis, living organisms can adapt to different environments. this model of metabolism as composed by independent metabolic pathways is simplistic, since pathways are nested and interdependent. in fact, metabolism is a complex system and pathways interact with each other. the aim of the work presented here is to find all the viable sets of heterologous enzymes, which can produce a predefined target compound when added to the pool of endogenous enzymes of a given organism. our method enables a metabolic engineer to find all heterologous metabolic pathways producing a target compound, for instance liquiritigenin, from the endogenous metabolites of e. coli k- as shown in figure as depicted in figure our problem can be formulated as searching for all possible heterologous pathways linking a target compound to the endogenous metabolites of an organism. to this purpose we provide software tools that enable the discovery of potential pathways producing a target chosen by the user. more precisely the user enters a target compound, a chassis organism, and our software tools return a ranked list of pathways to be engineered into the chassis organism. to achieve this task we have developed an approach composed of three steps. in the first step, using a retrosynthesis software, reactions producing a target compound are iteratively searched backwards until the set of needed precursors only contains source metabolites. this first step returns a retrosynthetic network connecting a target compound to the endogenous metabolites of an organism. there may be several pathways in the retrosynthetic network linking the source metabolites to the target compounds and there is thus a need to enumerate all the possibilities. pathway enumeration is performed by in the second step. once the pathways have been enumerated, we evaluate in the third step the possibility to insert each pathway and its associated heterologous enzymes in the host organism. this step consist of determining the catalytic efficacy of the enzymes, the toxicity of the products and the coproducts, and the easiness of inserting the enzymes into the host. the efficiency of the pathways can then be further estimated by flux models for the cell metabolism such as flux balance analysis. we have already discuss elsewhere the first and third step, i.e. methods to generate retrosynthetic networks and methods to rank pathway efficiency. to apply these methods in the context of heterologous target production, we need a computationally fast method to enumerate all possible pathways. we address the enumeration problem in the current paper. different mathematical models that describe metabolism have been proposed. we distinguish two main families of approaches: the ones computing steady states of the fluxes of reactions and the ones based only the topology of the network. typically, steady states are studied and simulated by generating the flux space. of particular interest are the extreme pathways and the elementary modes, they both represent the smallest generating set of the flux space and they both are composed of independent non-decomposable pathways in the network. the differences between extreme pathways and elementary modes have already been discussed in details and these differences arise when dealing with reversible reactions. in the present paper we consider all reaction irreversible, and networks comprising reversible reactions are modeled by doubling each reversible reaction into a forward reaction and a reverse reaction. algorithms have been developed to enumerate both extreme pathways and elementary modes and these algorithm are all variants of the double description method, which enumerates all extreme rays of a polyhedral cone. the algorithms use as input a stoichiometric matrix representing the network and output sets of fluxes satisfying sv = one notices that extreme pathways and elementary modes while representing pathways do not directly enumerate all pathways linking a source set to a target set of compounds. however as shown in the subsection "enumerating pathways using the steady state approach" one can construct stoichiometric matrices where input fluxes are added to the set of source compounds and outgoing fluxes are associated to the target and heterologous coproducts such that the extreme pathways and elementary modes enumerated from these matrices do correspond to all pathways linking the source set to the target. while as mentioned above, the problem of systematically enumerating pathways for heterologous production in chassis organisms has not yet been addressed, there are methods based on the steady state approach to search for heterologous pathways optimizing target productions, and methods to search for shortest pathways between source and target sets of compounds and. all these methods are based on optimization and make use of integer linear programming. precisely, the method of pharkya et al., is aimed at redesigning microbial chassis organisms through heterologous reaction addition and native reaction deletion for the overproduction of a target compound. the addition and deletion are parameterized using binary variables attached to each reaction. a mixed integer linear program is then set up to maximize the target yield while minimizing the number of added reactions. the methods of de figueiredo et al., and pey et al. are both aimed at searching for the k shortest pathways. in de figueiredo et al. the k first shortest pathways are searched in entire metabolic networks, while in pey et al. the pathways are searched between a source metabolite and a target metabolite. both methods solve the problem at steady state and search for fluxes, v, verifying sv = while minimizing the number of reactions turned on. aside from the fact that integer linear programs suffer from computational complexity all the above methods search for at most k optimized pathways and do not guarantee a full enumeration of the possibilities. in our methods the optimal pathways are computed in a post process by ranking the pathways that have been enumerated. our approach allows one to decouple enumeration from optimization, and thus to plug any optimization criteria, including nonlinear functions and not only target yield or pathway length. aside from using extreme pathways and elementary modes, we also present in this paper a topological model which directly enumerates all the possible heterologous pathways linking target compounds to a source set of compounds. the main advantage of the topological approach compared to the stationary state approach is computational speed. speed is in fact an important aspect when searching for the best pathways to engineer, as there are generally a combinatorial number of pathways between given source sets and target sets. as an illustration of this combinatorial complexity, the work of hatzimanikatis et al., which provides a list of novel biochemical routes from chorismate to phenylalanine, and the work of cho et al., which enumerates reaction routes to produce isobutanol. there exist standard graph-based methods to search and eventually enumerate pathways in metabolic networks, but these methods including pathfinding and pathway hunter tool are computing pathways and shortest pathways in graphs instead of hypergraphs. the particularity of these techniques is that only main substrates and main products are taken into account when constructing pathways, and consequently these main compounds must be differentiated from the cofactors. in the work of croes et al. cofactors are filtered out based on their connectivity in the network. indeed, compounds highly connected such as atp, nadp, or h2o are cofactors of most reactions as they do not share carbon atoms with the products of the reactions. in a more recent work, the main compounds in the pathways linking source metabolites to target metabolites are detected using the kegg rpair annotation, which enables one to follow the fate of atoms when going from a set of substrates to a set of products. another approach to search for main substrates and main product is the one developed with the pathway hunter tool, which consists of mapping substrates to products using cheminformatics fingerprints. while all the above techniques are computationally efficient, their main shortcoming is that they are not able to encompass reactions when a main product is formed from two main substrates. there are plenty of such reactions in metabolic networks, consider for instance the formation of guanidinoacetate from arginine and glycine through a glycine amidinotransferase, or the formation of glutathione from γ-l-glutamyl-l-cysteine and glycine catalyzed by a glutathione synthase. recently, some of the above topological methods have been benchmarked against the integer linear programming technique mentioned above to search for the shortest pathways linking various compounds, the recovery ratio for a set of predefined reference pathways could not reach 100% with the graph based approach, exemplifying the shortcoming of that approach. as reviewed above, while there are methods and theoretical results to enumerate elementary modes or extreme pathways and graph based techniques to search for pathways in a given metabolic network, to the best of our knowledge there is no known methods to directly enumerate pathways in the context of metabolic engineering, that is, to enumerate all the pathways encompassing all substrates and products necessary and sufficient to produce a given set of target compounds from a given set of source compounds. in the present paper we address that specific problem and present two methods one based on elementary modes and one based on a direct enumeration algorithm. in order to address this problem we need, in addition, to consider the problem of determining supplement molecules, i.e. metabolites that the organism cannot synthesize, but which can be added to the growth media in order to increase the number of viable pathways; and bootstrap molecules, i.e. metabolites which are required fist in order to be produced. while in the general context of metabolic network analysis, finding elementary modes does not require to first search for bootstrap molecules, in the context of metabolic engineering however any heterologous pathway solution that comprises a compound that is first consumed before being produced is valid only when the compound is added to the growth medium. therefore, in our study in the context of metabolic engineering, there is a need to first compute the bootstrap compounds prior to elementary modes. the paper is divided as follows. in the methods section we first provide some definitions, then outline our algorithms to solve the pathway enumeration problem with both the steady state approach and the topological approach. the problem of finding and enumerating all the pathways going from a large source to a target chosen by the user is considered. all the algorithms presented for the topological approach have polynomial worst-case running time, the algorithm for enumeration is a polynomial time per output algorithm on some classes of hypergraphs. we also provide algorithms to determine supplements, which are metabolites that the organism cannot synthesize, but which can be added to the growth media in order to increase the number of viable pathways. furthermore, an analysis of pathways containing supplements allows finding out pathways that contain bootstrap molecules. in the results and discussion section we illustrate our algorithms with the enumeration of the possible pathways to synthesize more than compounds in e. coli. we "experimentally" probe the computational complexity of the steady state and topological approaches for a series of networks of growing sizes and discuss the theoretical complexity results of the topological approach, which are provided in appendices a and b. while we illustrate our two methods for the production of heterologous compounds using as a source set all the endogenous metabolites of e. coli, our methods can be applied to any chassis organism and more generally to any source set. methods in the context of metabolic engineering, metabolic networks have been represented as directed graphs. in such graphs, edges are directed and correspond to reactions connecting two compounds if one is the product of the other. directed graphs can represent monomolecular reactions, but they are not well suited to capture more complex reactions. as already discussed in the background section, when representing bimolecular reactions, one has to choose which molecules are connected by the edges of the graph and which ones have to be excluded from the graph because they are considered co-substrates or co-products. additionally, one of the limitations of a model based on a graph representation is that depending on the criteria used to identify the co-substrates and co-products in the reactions, the networks obtained are different. in the present paper to palliate the limitations of the directed graph model, we represent networks as directed hypergraphs. the first examples of modeling through hypergraphs can be found in. in a hypergraph, each hyperarc connects a set of vertices, corresponding to reactants, to a disjoint set of vertices, representing the products. in our model each hyperarc corresponds to a reaction that can be catalyzed by an enzyme. it is worth noticing that hypergraph models have already been used to find minimal sets of metabolites sufficient to produce a set of target metabolites. unfortunately, the algorithms proposed in, do not enumerate pathways and are therefore not directly applicable to our metabolic engineering problem. definitions definition. a directed hypergraph is a pair h= where v = {v v..., vn} is the set of vertices and e = {e e..., em} is the set of hyperarcs. a hyperarc ei is an ordered pair ei = of disjoint subsets of v. the set xi is also called the tail of ei and the set yi is called the head, with reference to the graphical representation of arcs and hyperarcs as arrows. we denote by x:e→p the application that given an hyperarc ei returns its tail x ⊂ v. analogously we use y:e→p for the application that given a hyperarc returns its head. definition. in a metabolic network each vertex corresponds to a metabolite and each hyperarc corresponds to a reaction. a metabolic network of m metabolites and n reactions can be represented with a m × n stoichiometric matrix s, where the rows correspond to the m metabolites and the n columns to the reactions. a reaction j is represented by the column vector sj = t where sij is the stoichiometric coefficient of metabolite i in reaction j. reactants have negative coefficients and products have positive coefficients. examples of hypergraph, network, and stochiometric matrix are given in figure 2a. we notice that the stoichiometric coefficients of the reactions are not taken into account in the hypergraph representation. we also notice that the pair is ordered so to make the distinction between reactants and products. in this representation reactions are irreversible. many biochemical reactions can be considered as irreversible, since in organisms the homeostatic equilibrium is often strongly polarized. nonetheless, metabolic network may comprise reversible reactions, and we model these reactions by introducing both hyperarcs:, and. hyperpaths, a generalization of simple paths in graphs where cycle free paths going from one vertex to another, are used to represent pathways. a hyperpath connects a source set of vertices to a target set of nodes. two examples of hyperpaths are given in figures 2a and 2c. we remark that in a natural way a set e of hyperarcs defines a hypergraph ε = ∪ ∪e∈ey, e). by abuse of the terminology we denote by e the hypergraph corresponding to the set e of hyperarcs and all the heads and tails of the hyperarcs in e. the following definition for hyperpaths is borrowed from nielsen et al.. definition. a hyperpath p going from a source subset sh of v to a target subset tp of p in a hypergraph h= is a hypergraph hp= with vp ⊆ v, ep ⊆ e, such that there is an ordering f of the hyperarcs ep with the following properties. • ∀k∈{ …,|f|},x⊆sh∪) • tp⊆sh∪) from the point of view of metabolism, the first condition corresponds to the requirement that reactants of reactions participating in the hyperpath can be produced without the presence of the reaction itself. hyperpaths defined in this manner represent a metabolic route from the source to the target. according to definition the hypergraph of figure 2b with source a is not a hyperpath because neither reaction r nor r can happen until the other does not start. the definition, though complex, is computationally tractable, meaning that the time required to determine if a hypergraph is a hyperpath is proportional to the number of reactions. a polynomial time algorithm to determine if a hypergraph is a hyperpath is given in, the algorithm findall presented below can also be used for that purpose. in fact, as discussed below, if the set of reactions returned by findall contains all the reactions in hp, then hp is a hyperpath. the metabolic network described by a hypergraph has to be as comprehensive as possible, containing every known enzyme-catalyzed reaction occurring in organisms. we say that a hyperpath produces a set of target metabolites if it contains all those target elements. a set of target compounds is said to be reachable from a given source, or linked to the source, if there is at least one hyperpath producing the targets. we are interested in the enumeration of pathways leading to the production of a desired compound. hyperpaths do not generally give the best representation of pathways because hyperpaths can contain reactions not necessarily linking the target to the source. minimal hyperpaths, cf. definition, are an appropriate representation of pathways since they contain only the essential reactions linking the source to the target. in the definition given below, we say that a hyperpath p is a subset of another hyperpath p′ if v ⊆ v' and e ⊆ e'. for instance the hyperpath of figure 2c is a subset of the one of figure 2a. definition. a hyperpath with target tp is said to be minimal if it has no proper subsets with the same target. the target is disconnected from the source if a reaction is removed from a minimal hyperpath. in this sense minimal hyperpaths cannot be reduced. from a metabolic engineering perspective the concept of minimal hyperpath is useful as it defines the minimum set of reactions necessary to produce a target heterologous compounds, and consequently the minimum set of enzymes needed to be inserted into the chassis organism where the compound is going to be produced. in the following we define b to be the set of all molecules linked to the source for a given hypergraph h and source set sh. the characterization of b is the first task to be solved before the enumeration. once this set is known all the minimal hyperpaths can be enumerated for all the molecules associated to the vertices in b. supplements supplements for a target are molecules whose presence in the source set increases the number of pathways for target production. finding supplements is an important improvement when exploring ways to produce the target, since they make possible new pathways. for each target of interest one can look for vertices that once inserted in sh give place to pathways otherwise impassable. in terms of metabolism we are looking for the "supplement" molecules, i.e., molecules that once introduced in the source set permit to find more pathways than those otherwise available. we introduce below findsupp, an algorithm that returns the supplements. an analysis of pathways containing supplements allows to find out pathways containing bootstrap molecules, i.e. metabolites that are needed in reactions producing compounds afterwards used for the production of the bootstrap molecules. as a matter of fact, many pathways can be made viable once bootstrap molecules become available in the metabolic network. loosely speaking bootstrap molecules are molecules that cannot be produced by the reactions belonging to a hyperpath unless they are already present in the source. cottret et al stated that given a source set the existence of a pathway making use of bootstrap molecules can be tested in polynomial time. we provide later in this section an algorithm returning the bootstrap compounds, such algorithm can be used to determine if a target molecule is connected to the source through a pathway making use of bootstraps. enumerating pathways using the steady state approach in steady state, all possible pathways in a metabolic network are by definition stoichiometrically balanced, i.e. all metabolites produced from the source set must be consumed except for those that are target products. extreme pathways and elementary modes are two methods that compute the set of independent non-decomposable pathways in the network that generate all feasible steady state solutions in the flux space. they do not directly enumerate all pathways linking a source set to a target set of compounds. however, one can construct stoichiometric matrices where input fluxes are added to the set of source compounds and outgoing fluxes are associated to the target and heterologous co-products such that the extreme pathways and elementary modes enumerated from these matrices can be used to generate all pathways linking the source set to the target. given a hyperpath hp= of a hypergraph h=, we can define a set of flux vectors vp for the hyperpath where components vpj corresponding to those reactions in the pathway ej∈hp are activated: vpj= >0ej∈hp0ej∈h\hp a hyperpath hp= of a hypergraph h= with input source subset sh and the target subset tp is defined as stoichiometrically balanced if the rows corresponding to each metabolite vi ∈ v that are obtained from the product of the stoichiometric matrix s and the associated flux vector vp verify: svp=≤0vi∈sh≥0vi∈tp0vi∈v\sh,tp a way to introduce the constraint on input and output metabolites in the previous equation is by adding to the stoichiometric matrix s additional columns corresponding to input reactions, and output reactions. these auxiliary reactions, even if non-properly balanced in terms of the law of conservation of mass, are useful in order to define completely the problem in a compact manner: sv=0v≥0v∈r both extreme pathways and elementary modes make use of this formulation in order to compute the set of feasible solutions v. since in our hypergraph definition all reactions are irreversible, the set of pathways solving equation computed by both extreme pathways and elementary modes are identical. furthermore, solutions in v must contain only positive or null fluxes. in order to determine all stoichiometrically balanced heterologous pathways hp that can be inserted into the chassis organism to produce a target set tp, we need to constrain the computation of elementary modes only to those that have non-zero fluxes for heterologous reactions. efficient solutions to this problem have been considered in the divide-and-conquer approach by rearranging the constraints in an echelon form so that the constraints containing only the desired reactions appear at the bottom. to define the constraints in our case, we consider first the hypergraph rt that is formed only by heterologous reactions. this hypergraph rt is the subset of the hypergraph r formed by those hyperedges where at least one vertex v does not belong to the source set sr, i.e. those metabolites endogenous to the chassis organisms. by considering rt instead of the full hypergraph r, we are looking only for biosynthetic pathways involving heterologous reactions and therefore avoiding cycles internal to the chassis organism. therefore, to compute all feasible steady state heterologous pathways, we reformulate equation so that the stoichiometric matrix s is defined by reactions in rt; the input is given by all substrates in the source set sr∩x; and the output by all products of the reactions in the hypergraph y. finally, from the computed set of solutions v for equation we are interested in enumerating all minimal hyperpaths from sr to the target set t on the hypergraph given by rt. according to definition minimal hyperpaths for some target t are given by those cycle-free solutions in v containing only reactions linking the source to the target. since any feasible flux pattern v is a superposition of elementary modes with non-negative coefficients, the set of minimal hyperpaths for a given target t is a subset of the elementary modes producing t that are solution of equation namely, any feasible solution generated from the elementary modes will contain at least as many reactions as the ones that are in those elementary modes that form its basis. therefore no additional minimal hyperpaths can be generated in this case by superposition of elementary modes. enumerating pathways using the topological approach the algorithm findall that allows to find b, the set of metabolites that can be linked to the source sh by a hyperpath. findall, by explicitly constructing the ordered set f in definition, provides a proof of the tractability of the problem of checking if a hypergraph is a hyperpath. moreover findallf permits to prune the original hypergraph enabling a faster enumeration algorithm. as presented below the algorithm minimize, when called on the output of findall, returns, if exists, a minimal hyperpath linking a given target to the source. these algorithms are the main components of the algorithm enumerating the pathways findpath described next. then we present findsupp an algorithm to enumerate supplements. finding one minimal hyperpath let h= be the hypergraph representing the set of metabolic reactions, n = |v|, m = |e| and let sh be the set of source vertices representing the source metabolites. the algorithm findall returns all the reactions that can contribute to the production of any element in b, i.e., the set of all compounds that can be connected to the source. findall is a linear algorithm in the number of vertices, hyperarcs and in the total coordination; the complexity is o |+ |y-1|) that is bounded by o. therefore, such algorithm can be applied to the hypergraph h of all reactions in order to obtain a pruned sub-hypergraph h′= where the set of vertices v′:=sh∪b, and the set of edges e' is the set of reactions returned by findall. in the context of metabolic engineering findall returns all the compounds that can be produced from a given set of source compounds and reactions. for instance, using findall with all know metabolic reactions one can determine all the compounds that can be produced from the metabolites of e. coli. algorithm findall input: h,sh for all r in h    x ← x end for v←sh d←sh f ← {∅} while v ≠ {∅}    let i be an element of v    v ← v \ i    d ← d ∪ i    for all r ∈ h such that i ∈ x:       x ← x \ i       if x = {∅}          f ← {f, r}          for all j in y and not in d:             v ← v ∪ j          end for       end if    end for end while output:    f let d be the union of the source set and of the heads of all the reactions output in f by findall. the correctness of the algorithm above is given by the following claims: every element in d is the target of some hyperpath or is part of the source, and every vertex in h that can be reached from the source is in d. for the first claim we can give a constructive proof by using the output vector f, the second claim is proved by contradiction. • the proof of the fact that every element in d is reachable from the source is given constructively by the ordered set f returned by the algorithm. in fact at each step f is a hyperpath. this claim can be proved by induction on the steps of the algorithm, each time a hyperarc r is appended to f the tail x is contained in d and if a vertex j is added to d it means that it was previously added to v and thus it was in the head y of some hyperarc already in the hyperpath. • the second claim can be proved by contradiction: if an element of b were not in d there would be a hyperpath linking it to the source. in such hyperpath let consider the first reaction r whose x belongs to d and such that one of the elements of y does not. consider among x the last one that has been inserted into the set v; after its removal from x this set becomes empty and the elements of y are inserted into v and then in d, which is a contradiction. from the above statements follows that each vertex appearing in a hyperpath having as source sh is an element of d and every hyperarc is an element of f. thus the algorithm findall provides an effective pruning of the original hypergraph: in h there is no minimal hypergraph with source sh containing hyperarcs not in ∪kfk or vertices not in b. the output hyperarcs are the only ones that can belong to a minimal hyperpath, and h′=),∪kfk) is the pruned hypergraph only containing reachable vertices and hyperarcs. notice that findall algorithm as presented above returns in polynomial time a hyperpath valid for each target vertex in b. even though there are more efficient algorithms for finding a hyperpath for one single target, for the sake of simplicity we avoid to introduce here an additional algorithm and just remark that since findall is polynomial, the use of it does not affect the complexity analysis of the algorithms making use of its output. remark that a minimal hyperpath going to a specific target can be easily extracted from the hyperpath output of findall. namely, given a hyperpath p connecting s to t, it is always possible to find a minimal hyperpath p′ subset of p. moreover it can be done in polynomial time, for instance by using minimize, the algorithm introduced below. minimize is an algorithm that takes as input a hypergraph p, a hyperpath rf subset of p, a target set of vertices t and a source s. if p does not link t to s the empty set is returned, otherwise a hyperpath contained in p, containing rf and linking t to s is returned. in particular, if rf is empty, the output of minimize is the minimal hyperpath going from s to t, provided it exists. minimize returns a hyperpath obtained by removing all inessential hyperarcs except for the ones in rf. in the context of metabolic engineering, pathways containing a small number of heterologous reactions are generally preferred, since they are easier to engineer in the host organism. therefore, given two pathways that produce the same target, where one is subset of the other, the one requiring the smaller number of heterologous reactions has to be selected. this is the reason that makes relevant to obtain minimal hyperpaths from generic hyperpaths. algorithm minimize input: p,rf,t,s f ← findall p' ← p if t ⊄ ∪k y p′←{∅} else    for all r in p       if r not in rf          f ← findall          if t ⊂ ∪k y             p′←p′\r          end if       end if    end for end if output:    p′ the proof of correctness of this algorithm is simple and is based on the fact that p′⊆p implies findall⊆findall. if a reaction in p has not been removed from p′, then any subset of p not containing r does not produce the target. the worst-case time for this algorithm is o|+|y|)). since x and y have bounded values, the algorithm has a quadratic complexity. even though faster algorithms can be designed, here we presented this one because of its conceptual simplicity. remark that, since minimize returns a minimal hyperpath p′ subset of p if it exists, then the minimality of a hyperpath p can be tested by checking whether p′=p or not. a related problem to minimize is the minimal constrained hyperpath problem: the problem of finding if a minimal hyperpath from a given source to a given target, containing the hyperarcs in rf exists. notice that minimize, although linked to this problem does not solve it. in fact, if the output of minimize is an empty set then there are no minimal hyperpaths satisfying the constraints; however if the output is a minimal hyperpath then obviously a minimal constrained hyperpath exists; and finally, if the output is a non-minimal hyperpath then we do not know if a minimal hyperpath satisfying the constraints exists or not. below we will discuss why we are interested in algorithms for the minimal constrained hyperpath problem, while in appendix a. we show that in general the problem is np-complete. pathways enumeration the basic idea behind the enumeration algorithm presented below is to introduce an iterative refinement of partitions of the space of feasible solutions i.e. of the space of hyperpaths and in each part to look for a solution. in our implementation, a part is defined by two sets of reactions of the original hypergraph. these sets are used during the enumeration process, rf is a set of hyperarcs that must be present in the enumerated hyperpath and rn is the set of hyperarcs that must not be part of the enumerated hyperpath. the problem of finding a solution in one of the parts is addressed at each iteration and if a solution is found the part is divided in finer parts. this process is repeated until all the minimal hyperpaths have been found. enumeration by means of the minimal constrained hyperpath problem first we describe informally the enumeration algorithm through the toy example hypergraph in figure 2d and 3a, then we outline in figure 3b a typical run for a more involved example: liquiritigenin. a minimal hyperpath p connecting the node v to the source nodes v v on the hypergraph h of figure 2d can be obtained by calling minimize on the hypergraph p′ obtained by findall. the hypergraph p′ is represented in figure 3a. once p1={r r3} has been obtained, the search space is divided into three parts: • the hypergraphs which do not contain r • the hypergraphs which do contain r and do not contain r • the hypergraphs which do contain r and r the first set does not contain hyperpaths connecting the target to the source: once the reaction r is removed, v is disconnected from the source. the second set contains a solution and thus has to be partitioned. the third set contains only one minimal pathway. the minimal hyperpath in the second set is found by running findall on h\r and then minimize with constraint rf = {r4}. the minimal hyperpath so obtained is the one only containing hyperarcs r r the set of the hypergraphs defined by is partitioned in two parts defined by new sets of constraints. the way the partition is done is explained in detail in algorithm findpath and gives two non overlapping sets: • the hypergraphs which do contain r and do not contain r and r • the hypergraphs which do contain r and r and do not contain r the first of these sets does not contain hyperpaths going to v8: once r and r are removed, node v is disconnected from the source. the second one only contains the second and last minimal hyperpath: the one consisting of hyperarcs r r the algorithm here sketched is based on the fact that all minimal hyperpaths are found once the problem of finding a minimal hyperpath has been solved for each part of the partition. relaxed hyperpath minimization the enumeration procedure is performed by the algorithm findpath, which enumerates all the minimal pathways and does not output duplicate hyperpaths. precisely findpath returns a set of hyperpaths containing all the minimal hyperpaths in h connecting t to sh and containing all the reactions in rf. findpath returns all the minimal hyperpaths from the source sh to the target t in h. a schematic representation of how findpath works for the enumeration of the pathways of liquiritigenin is given in figure 3b where we represent each call with a box connected by an arrow to its parent process. for each call of findpath either a new hyperpath is found and then findpath is executed with new constraints, or there are no new hyperpaths and the branching process is stopped. the new constraints sets rf', rn' for a new call of findpath are obtained by incrementing the sets rf, rn of the father process. given an order for the hyperarcs of the hyperpath p found for the father process, the set rn' relative to the child process is constructed by incrementing rn by one element r belonging to p, the set rf' is constructed by incrementing rf by all the hyperarcs coming before r. for each element in p not belonging to rn a child process is called. findpath returns all the minimal hyperpaths from the source sh to the target t in h. in the context of metabolic engineering findpath returns all the metabolic pathways for the production of the target compounds. algorithm findpath input: h, rf, t, sh f ← findall p ← minimize en ← ∅ if p≠∅    en←p       f←findall       for all k in |f|,..., 1}          r = fk          if r not in rf:             en←{en,findpath}             rf ← rf ∪ r          end if       end for    end if output:    en the loop at line of findpath is done according to the order given by line where the hyperarcs are ordered so that at least one of the head vertices of each hyperarc is a tail vertex of some previous reaction. such an ordering is always possible since p is a hyperpath. as said above and illustrated in figure 3b, findpath is an algorithm that iteratively calls itself, see line note that even if rn is not explicitly defined in findpath, it is constructed implicitly when at line of findpath is called on the smaller graph h\r. let us note that the output of the enumeration is not always composed of minimal hyperpaths. this is due to the fact that the algorithm minimize while running in polynomial time can return a non-minimal hyperpath. an algorithm always returning minimal hyperpaths cannot be polynomial since the problem of finding a minimal hyperpath containing a set rf of hyperarcs is an np-complete problem as showed in appendix a. however, in many practical instances, the algorithm minimize returns a minimal constrained hyperpath. as a matter of fact, for all the enumeration studies we have so far carried out, we observed that the output obtained by minimize when called by the algorithm findpath introduced above was always a minimal hyperpath. nonetheless, a characterization of hard instances of the minimal constrained hyperpath problem is given in appendix. supplements enumeration provided a given metabolic network and a set of source compounds it may not be possible to link all the metabolites of the network to the source set. when a target compound is not accessible from the source set, one can consider the possibility of inserting into the metabolism of the organism some precursors so that the target becomes reachable. in practice such a task can be carried out through the enrichment of the growth media. more generally, the insertion of supplements can be used even when the target compound is reachable in order to access to new pathways for the production of the target. let a supplement for a target t be any compound i∉b that is involved as reactant in at least one minimal hyperpath going from a superset of sh∪i to the target t. below we give the algorithm findsupp finding the supplements for the production of a given target. in figure 2d supplements are highlighted in red. therefore, the process of finding supplements is use-ful as a general strategy in metabolic engineering in order to determine which metabolites might be part of the metabolism that produces a given target. algorithm findsupp input: h, sh, t wishlist ← t d ← {∅} while wishlist \ d ≠ {∅}    let i be an element of wishlist \ d    d ← d ∪ i    aux ← {∅}    for all reactions r with i ∈ y       aux←aux∪\)    end for       wishlist ← wishlist ∪ aux end while f←findall d←d\sh∪) output: d bootstraps bootstrap molecules relative to a source sh are the molecules that cannot be produced by a hyperpath with source sh unless they are already present in the media. an example of bootstrap nodes are nodes v v of figure 2d. in this section we give an algorithm finding in polynomial time all the bootstraps of a hypergraph h with source vertices sh. bootstraps are special kind of supplement, if at any step of a pathway, a heterologous metabolite is needed as a substrate but has not yet been produced from the source set, then this metabolite is a bootstrap and must be added in the growth media for the reaction to take place, and for the pathway to be a valid pathway. the algorithm given below enables one to detect bootstraps prior enumerating pathways running the findpath algorithm. algorithm findbootstraps input: h,sh f←findall d←sh∪) h′←{∅} for all r in h    r' ← \ d, y \ d)    if y ≠ {∅}:       h′←h′∪r′    end if end for while exists v in ∪r∈h′y\∪r∈h′x    for all r' containing v:       r' ←, y \ v)       if = {∅}) or ):          h′←h′\r′       end if    end for end while b=∪r∈h′y output:    b the findbootstraps algorithm is linear in the number of vertices, hyperarcs and in the total coordination. remark that the set ∪r∈h′y obtained in line is equal to ∪r∈h′x. in fact the bootstrap vertices b constitute the largest set of vertices not reachable from the source and such that each element of the set belongs to the head of at least one reaction whose tail only contains vertices in b or in b. notice that the set of bootstrap vertices in a hypergraph h only depends on the source vertices and does not depend on the target.  to evaluate the performance of our topological approach, we have compared running times of this approach with the running times of the steady state approach presented in the methods section. all our tests were run on a mac pro server × ghz quad-core intel xeon, gb. all the algorithms of the topological approach were implemented in python. for the steady state approach we used two software products, one computing elementary modes and the other extreme pathways. elementary modes were computed by using the matlab interface to the java implementation of efmtool version.. extreme pathways were computed by using the mac osx version of the expa program. the running time comparison test was performed for different random samples of the hypergraph h constructed from the kegg database. e. coli was chosen as source organism. the hypergraph h was composed of metabolites connected by reactions including metabolites endogenous to e. coli and heterologous compounds. each sampled hypergraph hs was built by randomly sampling a given fraction of the total reactions in the hypergraph h. tests for each sampling fraction were repeated times. in the case of elementary modes, we were able to run the test only up to 50% sampling of the full metabolic network, due to memory constraints in matlab. for extreme pathways, the test was run up to 85% due again to memory constraints as well. prior to enumerating pathways with elementary modes, extreme pathways or the direct algorithm, bootstraps were identified by using the algorithm findbootstraps defined in methods and added to the growth media. the observed trend is consistent in the three cases, as it is shown in figure for the network composed of all metabolites and reactions, there are pathways connecting out of heterologous metabolites to the source. this relative low number of pathways producing heterologous compounds is related to the fact that heterologous maps usually show a degree of hierarchy higher than the one observed in native metabolic networks such as central metabolism. it is indeed not surprising to find the same number of extreme pathways and elementary modes as all reactions in our networks are irreversible. in the case of findpath, the fact that the topological approach does not take stoichiometry into account might produce some inconsistent pathways. checking that a given pathway is stoichiometrically balanced can be done using linear programming. in the context of metabolic engineering we found only few inconsistent pathways. these pathways are eliminated after enumeration by our ranking function as one of the criteria to rank pathways is to solve the steady state equations. once unbalanced pathways were removed from the list of enumerated pathways, we obtained the same pathways as through the calculation of elementary modes. these results are consistent with the given definitions of the algorithms. in a forthcoming paper, we plan to generalize them into a formal proof. we also find that findpath has the fastest execution time, the algorithm required less than min to compute the full network. elementary modes has an execution time approximately 10-fold slower than findpath. the computation of extreme pathways, in turn, is the slowest one, being between ~ 103-fold slower than findpath. the execution times are on average of ± s per output with the findpath algorithm, of ± s with elementary modes, and of ± s with extreme pathways. the observed execution times mean that for the compound with the highest number of pathways in our tests, the anti-diabetic drug acarbose containing pathways, enumerating the pathways with elementary modes takes more than seconds, while it can be computed in seconds with findpath. we have measured running times per output and memory usage as function of input size and output size. in all cases we found linear growth o for both input size and output size. in the case of running times, findpath and efmtool have running times per output approximately constant in function of the size of the input stoichiometric matrix s of × 10- and × 10- seconds, respectively, being expa the less efficient code with a constant time of × 10- seconds and a linear growth of × 10- seconds per input. similar values were obtained depending on the size of the output, although the scaling in this case for expa was of × 10- seconds per output. regarding memory usage, expa and findpath are less demanding, especially for output size, with linear growths of and kb per output, respectively, while efmtool required significantly more memory allocation, × kb per output. performance of the codes findpath, efmtool and expa in terms of run time per output and memory use in function of size of input stoichiometric matrix s and size of output. to the best of our knowledge, the computational complexity of enumerating elementary modes on networks comprising irreversible reactions is up-to-date unknown. in appendix a. we prove that enumerating minimal pathway is an np-complete problem but as examplified in appendix b the hard instances are rather rare and obviously none of these hard instances were found in our running tests as all pathways returned by findpath were indeed minimal pathways. to further probe the computational complexity of findpath we computed distribution of the run time for each of the heterologous targets that are linked to e. coli. as shown in figure the distribution is exponentially distributed, the average run time is thus finite. examples of pathway enumeration for heterologous compounds with therapeutical applications were provided in our previous methodology study about pathway ranking for penicillin n and taxol. these two examples contained a relative low level of combinatorial complexity, with a maximum number of different pathways producing penicillin n. the combinatorial complexity of enumerating all putative biosynthetic routes producing a target compound, however, might become considerably higher than in these two previous examples. further examples are those involved in the biosynthetic pathways for plant steroids leading to brassinolide, where pathways can be used to produce campesterol, one of the initial precursors. the number of pathways grows as we proceed downstream going up to for teasterone, and to for typhasterol. finally, the number of pathways for the end products, castasterone and brassinolide is of this example illustrates how the complexity of pathways enumeration can grow with the number of intermediates involved in the synthesis of the final product. since findpath algorithm has been designed with the aim to enumerate pathways in metabolic engineering applications, it does not directly address other general metabolic network analysis problems like finding shortest paths between two compounds of endogenous pathways in central metabolisms. indeed, findpath does not enumerate pathways comprising compounds that are consumed before being produced unless they were added as supplement in the growth medium. another aspect of pathway design is the definition of a cost function that estimates limiting effects on production efficiency by multiple factors associated with genetic and metabolic engineering of the pathway, such as gene heterogeneity, metabolite toxicity, or steady-state fluxes. in order to facilitate the designer in the selection of the best synthetic pathways to be implemented each hyperpath enumerated by our findpath algorithm is ranked depending on this cost estimation. once a cost function to minimize is introduced, the search for an optimal pathway can be formulated as a shortest path problem. the shortest path problem for weighted graphs consists of finding a path going from a given source vertex to a given target vertex while minimizing a cost function given by the sum of the costs of the arcs involved in the path. in order to define a shortest-hyperpath problem we need a definition of cost for hyperpaths in a weighted hypergraph. a natural generalization of the cost function for paths in graphs to hypergraphs is the sum of the costs of the hyperarcs contained in the hyperpath. if on the one hand this generalization seems to be natural, on the other hand the two problems have, however, not the same complexity. in fact the shortest-hyperpath problem with this cost function is known to be an np-hard problem. the reason why the algorithms commonly used for graphs cannot be easily adapted for hypergraphs is that given the costs w for the hyperpaths ai⊂p leading to the vertices ai, tails of a hyperarc ej whose cost is wj, then the cost of the hyperpath containing all the hyperarcs in ∪iai and ej is not equal to wj+ ∑iw if the overlap of the hyperpaths ai is non trivial. provided an additive cost function as defined in, finding the shortest hyperpath can be done in polynomial time by an algorithm of complexity o) given in, which finds the shortest hyperpath in b-hyperpaths. in fact, this algorithm does apply to minimal hyperpaths given in definition. generally speaking, the strategy we used in the enumeration algorithm is often used to solve the problem of finding the first k solutions of a combinatorial problem, but it cannot be extended to the k-shortest problem on hypergraphs not even with an additive cost function, since when splitting the problem into a constrained part and a free one we get a problem that is known to be np-hard. therefore, the k-shortest path problem might be solvable efficiently if one is able to develop a specific cost function making the problem tractable. however such an artifact cost function may not be necessarily appropriate to practical problems, such as the production of heterologous targets using metabolic engineering. as an alternative we have proposed in this paper an enumeration algorithm, which systematically enumerates all pathways linking source compounds to target compounds. in all the cases that we have so far processed for metabolic networks, our algorithm runs in polynomial time per output, and therefore pathways can be ranked by the designer based on their own user-defined cost functions.  in summary, the methods presented in this paper provide metabolic engineers with powerful tools that extend the toolbox for heterologous biosynthetic pathway design. besides pathway enumeration of biosynthetic routes for a given target product, our methods have several other possible applications. for instance, they can be used in combination with gene deletion strategies in order to determine pathway manipulations leading to overproduction of the target compound. another application is in biodegradation and bioremediation, where our algorithms would need to be slightly modified in order to reverse the pathway search so that it can identify degradation routes for a given compound, while the underlying structure of the algorithms remains still valid. finally beyond metabolism, our algorithms could also be utilized in the context of chemical synthesis to enumerate all the possible routes linking a target molecule to a source set of starting reactants, enabling the search for the best routes in terms of production costs. availability and requirements a web server is available: http://bioretrosynth.issb.genopole.fr/tools/metahype see details in additional file authors' contributions jlf and pc designed the algorithms and the experiments. df implemented the algorithms, and wrote the proof of the np-completeness for the minimal constrained hyperpath problem. sbp wrote the example of pathway enumeration. pc implemented the web service prototype, performed the experiments and analyzed the results. all authors wrote the manuscript. all authors read and approved the final manuscript. appendix a reduction proofs for the shortest hyperpath problem and the minimal constrained hyperpath problem a. shortest hyperpath problem in a reduction of the shortest-hyperpath problem to minimum set cover is given. we have to adapt the proof to our case for two reasons: the definition of directed hypergraph that was used is more restrictive | = 1) and the hyperpath was ill-defined. in fact the given definition by these authors does not permit to say if some hypergraphs are also hyperpaths or not. in other words their definition is ambiguous: does not permit to determine the nature of all the directed hypergraphs and thus can be completed in several ways. nonetheless the proof of np-hardness they gave is valid for our more general hypergraphs and minimal hyperpaths because the set of directed hypergraphs employed in the reduction proof in is a sub-ensemble of the directed hypergraphs we defined above in the definitions section and all the hypergraphs employed for the reduction are well defined as hyperpaths, independently of the way the incomplete definition they gave is completed. since our definition is a way to complete the definition in, then the two definitions agree on the set of hypergraphs employed for the reduction. from these facts follows that the shortest hyperpath problem is an np-hard problem. and, in particular, if the weights on hyperarc are non-negative, since hyperpaths always contain at least one minimal hyperpath, the shortest minimal hyperpath problem is np-hard too. additive cost functions the reason why the algorithms commonly used for the shortest path problem on graphs cannot be easily adapted for hypergraphs is that given the costs w for the hyperpaths ai⊂p leading to the vertices ai, tails of a hyperarc ej whose cost is wj, then the cost of the hyperpath containing all the hyperarcs in ∪iai and ej is not equal to wj+ ∑iw if the overlap of the hyperpaths ai is non trivial. in order to define a shortest path problem that can be solved polynomially by a variant of dijkstra algorithm, the additive cost functions have been introduced in for the b-hyperpaths. we adapt below the notion of "additive" cost function for hyperpaths. a cost function w returning the cost for reaching the vertex x with the hyperpath p starting from a source s whose elements s ∈ s have w:= is additive if w is the minimum over all the arcs ex∈p whose head contains x of ex + f where wi:=w are the costs for reaching the tail vertices i of ex in p, f is an increasing monotone function of its argument and f ≥ wi ∀i. remark that the cost of a hyperpath determined with an additive cost function in general is not given by the sum of the costs of the hyperarcs. a. minimal constrained hyperpath problem consider a 3-sat instance concerning n variables σi and m clauses xj consisting of the problem of deciding if there exists an assignment of boolean values to the σi such that all the clauses are satisfied. for each boolean variable σi contained in at least one clause introduce one hyperarc εi with the head of each εi having two vertices y = {vi+, vi-}. for each clause xj consider a vertex νj and seven hyperarcs. a boolean assignment is a triple a a a of boolean values. let these hyperarcs be μj..., μj and let the head y of a hyperarc μjk corresponding to the combination a a a of the boolean variables σj σj σj be y = {vj1a vj2a vj3a νj}. now let the tails of each hyperarc introduced be connected to the source nodes. and let consider a node t being the product of the reaction r having as substrates the vertices νj and the heads {vi+, vi-} of the hyperarcs εi. given the hypergraph described above consider the minimal constrained hyperpath problem where all the hyperarcs εi are mandatory, and the target is t. a solution of this problem gives in linear time a solution for the underlying sat problem, which makes the problem of minimal constrained hyperpath an np-complete problem. in fact, given a minimal hyperpath m, solution of this problem, for each i consider vi* the only one of the two head vertices {vi+, vi-} belonging to the head of one or more of the μ arc in m. the boolean assignments σi=vi* are a solution of the 3-sat problem and, inversely if a solution of the 3-sat problem exists then a minimal pathway satisfying the constraints exists and is the one obtained using only one of the hyperarcs for each xj among the ones whose head only contains νj and vertices vi*. in additional file figure s for simplicity we consider the reduction of a single-clause satisfaction problem to finding if a minimal hyperpath satisfying the constraints exists. there exist seven minimal hyperpaths connecting the target vertices to the source and satisfying the constraint that ε ε ε are parts of the hyperpath. each solution corresponds to a valid boolean assignment of the variables σ σ σ appendix b hard instances of minimal constrained hyperpath problem on many hypergraphs the algorithm enumerating the pathways only returns minimal hyperpaths, this is the case for the metabolic networks that we analyzed in the main sections of this paper. in this section we give a characterization of the hypergraphs where the algorithm minimize solves the minimal constrained hyperpath problem, characterizing these instances helps to individuate which hypergraphs are expected to give an output only containing minimal hyperpaths. let y be the set of all the metabolites produced by reactions in rf, the mandatory reactions: y:=∪r∈rfy. given a hypergraph h and the sets rf, rn we say that the well-separation condition holds if for every reaction r∈h\ the set y of products of r either is a subset of y or does not contain elements of y. if the well-separation condition holds for a hypergraph h with constrained reactions rf, the algorithm minimize returns a minimal hyperpath solving the minimal constrained hyperpath problem if a solution exists. the well-separation condition holds for every choice of rf in a hypergraph whose reactions have one only product, as the hypergraphs defined in. if on the one hand, this condition can appear too constraining, on the other hand it can be generalized to larger sets of hypergraphs. for instance, the algorithm minimize returns a minimal hyperpath solving the minimal constrained hyperpath problem even if the well-separation condition holds on the pruned graphs obtained by keeping from the original hypergraph only the reactions belonging at least to one hyperpath linking the target to the source and only the nodes being tail of such reactions. examples of hard instances can be found among the ones used for the proof of np-completeness. in general, hard instances h of the enumeration problem have to violate the condition of well-separation for same choice of rf⊂h, in order that the corresponding minimal constrained hyperpath problem becomes hard. this happens if several compounds are products of more than one reaction producing more than one compound. this is the case for nested networks as the one in additional file figure s while the given network is small enough to be solvable by hand, it contains nevertheless the principal ingredients of complexity that would asymptotically make harder the problem as the size of the instances grows. finding one minimal hyperpaths leading to the production of v is a simple problem, but finding new ones gets more and more involved. this is a consequence of the fact that the nodes v v v v v can be produced by different choices of the reactions r r r r r and each of these reaction has more of one product susceptible to participate to the production of the target. supplementary material additional file 1 figure s distribution of graph hierarchies in heterologous metabolic networks in comparison with the graph hierarchy of central, nucleotide, lipid, and amino acid metabolic networks in e. coli. click here for file additional file 2 figure s performance comparisons for findpath, expa and efmtools for run time per output and memory use in function of size of input and size of output. click here for file additional file 3 an example of pathway enumeration with bootstraps. an example of pathways enumeration with findpath using bootstraps. click here for file additional file 4 availability and requirements. description of the available metahype web server for running the algorithms for kegg compounds. click here for file additional file 5 figure s minimal constrained hyperpath problem reduction of a 3-sat formula. click here for file additional file 6 figure s instance of minimal constrained hyperpath problem. click here for file acknowledgements funding: genopole® through an atige grant; anr through a chair of excellence.
we consider the possibility of engineering metabolic pathways in a chassis organism in order to synthesize novel target compounds that are heterologous to the chassis. to select the best pathways to engineer, we have developed an objective function that computes the cost of inserting a heterologous pathway in a given chassis organism. for this purpose, we model metabolic networks through hypergraphs where reactions are represented by hyperarcs. yet both methods are generally scalable in time with the number of pathways in the metabolic network. we follow a retrosynthetic approach in order to search in the metabolic space for pathways linking the target compounds to a source set of compounds. in order to find minimum-cost pathways, we propose in this paper two methods based on steady state analysis and network topology that are to the best of our knowledge, the first to enumerate all possible heterologous pathways linking a target compounds to a source set of compounds. we also provide an algorithm for identifying precursors which can be supplied to the growth media in order to increase the number of ways to synthesize specific target compounds. in the context of metabolic engineering, the source set is composed of all naturally produced chassis compounds and the target set can be any compound of the chemical space. we find the topological approach to be faster by several orders of magnitude than the steady state approach. therefore this work provides a powerful tool for pathway enumeration with direct application to biosynthetic pathway design.
 we consider the possibility of engineering metabolic pathways in a chassis organism in order to synthesize novel target compounds that are heterologous to the chassis. for this purpose, we model metabolic networks through hypergraphs where reactions are represented by hyperarcs. each hyperarc represents an enzyme-catalyzed reaction that transforms set of substrates compounds into product compounds. we follow a retrosynthetic approach in order to search in the metabolic space for pathways linking the target compounds to a source set of compounds.  to select the best pathways to engineer, we have developed an objective function that computes the cost of inserting a heterologous pathway in a given chassis organism. in order to find minimum-cost pathways, we propose in this paper two methods based on steady state analysis and network topology that are to the best of our knowledge, the first to enumerate all possible heterologous pathways linking a target compounds to a source set of compounds. in the context of metabolic engineering, the source set is composed of all naturally produced chassis compounds and the target set can be any compound of the chemical space. we also provide an algorithm for identifying precursors which can be supplied to the growth media in order to increase the number of ways to synthesize specific target compounds.  we find the topological approach to be faster by several orders of magnitude than the steady state approach. yet both methods are generally scalable in time with the number of pathways in the metabolic network. therefore this work provides a powerful tool for pathway enumeration with direct application to biosynthetic pathway design. 