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 (a plant secondary metabolite with therapeutic applications), from the endogenous metabolites of E. coli K-12, as shown in Figure ​Figure11.

As depicted in Figure ​Figure1,1, 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 [1]. More precisely the user enters a target compound, a chassis organism, and our software tools return a ranked list of pathways (each list being composed of enzymes) 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 [2], 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 [3].

We have already discuss elsewhere the first and third step [1,2], 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 (cf. [4] for a review of the different models). We distinguish two main families of approaches: the ones computing steady states of the fluxes of reactions (one well-known application being the flux balance analysis) 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 (minimal) generating set of the flux space and they both are composed of independent non-decomposable pathways in the network [4]. The differences between extreme pathways and elementary modes have already been discussed in details [5] 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 [6] and elementary modes [7] and these algorithm are all variants of the double description method [8], which enumerates all extreme rays of a polyhedral cone. The algorithms use as input a stoichiometric matrix (S) representing the network (cf. [3] for definition of stoichiometric matrix) and output sets of fluxes (v) satisfying Sv = 0. One notices that extreme pathways and elementary modes while representing pathways (to each flux verifying Sv = 0 correspond a stoichiometrically balanced pathway) 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 [9], and methods to search for shortest pathways between source and target sets of compounds [10] and [11]. All these methods are based on optimization and make use of integer linear programming. Precisely, the method of Pharkya et al. [9], 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 (MILP) is then set up to maximize the target yield while minimizing the number of added reactions. The methods of de Figueiredo et al. [10], and Pey et al. [11] are both aimed at searching for the k shortest pathways. In de Figueiredo et al. [10] the k first shortest pathways are searched in entire metabolic networks, while in Pey et al. [11] 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 = 0, while minimizing the number of reactions turned on (using a binary variable). Aside from the fact that integer linear programs suffer from computational complexity (MILP is an NP-hard problem) all the above methods search for at most k optimized (shortest) 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 (cf. page 3 and Carbonell et al. [1] for a list of criteria entering our metabolic engineering optimization problem).

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. [12], which provides a list of 75,000 novel biochemical routes from chorismate to phenylalanine, and the work of Cho et al. [13], which enumerates 107,272 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 [14-16] and Pathway Hunter Tool [17] 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 (i.e. co-substrates and co-products). In the work of Croes et al. [14,15] cofactors are filtered out based on their connectivity in the network. Indeed, compounds highly connected such as ATP, NADP, or H₂O are cofactors of most reactions as they do not share carbon atoms with the products of the reactions. In a more recent work [16], the main compounds in the pathways linking source metabolites to target metabolites are detected using the Kegg RPAIR annotation [18,19], 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 (EC 2.1.4.1), or the formation of glutathione from γ-L-glutamyl-L-cysteine and glycine catalyzed by a glutathione synthase (EC 6.3.2.3). Recently, some of the above topological methods have been benchmarked against the integer linear programming technique mentioned above [11] to search for the shortest pathways linking various compounds, the recovery ratio for a set of 40 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 (steady state approach) and one based on a direct enumeration algorithm (topological approach). 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 [20]. 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 (as for instance al the metabolites of an organism) to a target chosen by the user is considered. All the algorithms presented for the topological approach (with the exception of the algorithm for enumeration) 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 5000 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.

A directed hypergraph is a pair where V = {v₁, v₂ ..., v_n} is the set of vertices and E = {e₁, e₂,..., e_m} is the set of hyperarcs. A hyperarc e_i is an ordered pair e_i = (X_i, Y_i) of disjoint subsets of V.

The set X_i is also called the tail of e_i and the set Y_i is called the head, with reference to the graphical representation of arcs (directed edges) and hyperarcs as arrows.

We denote by the application that given an hyperarc e_i returns its tail X(e_i) V. Analogously we use for the application that given a hyperarc returns its head.

Definition 2 (Reactions and networks).

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 S_j = (s_1j,..., s_mj)^T where s_ij 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 ​Figure2A.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 (X, Y) 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: (X, Y), and (Y, X).

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 ​Figures2A2A and ​and2C.2C. We remark that in a natural way a set E of hyperarcs defines a hypergraph ε = (_eEX(e) _eEY (e), 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. [22].

Definition 3 (Hyperpaths).

A hyperpath P going from a source subset of V to a target subset T_P of P in a hypergraph is a hypergraph with VP V, EP E, such that there is an ordering F of the hyperarcs EP with the following properties.

•

•

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 (3) the hypergraph of Figure ​Figure2B2B with source a is not a hyperpath because neither reaction R₁ nor R₂ can happen until the other does not start. The definition (3), 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 [23], 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 , then 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 (4), 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 is a subset of another hyperpath if V V' and E E'. For instance the hyperpath of Figure ​Figure2C2C is a subset of the one of Figure ​Figure2A2A.

Definition 4 (Minimal Hyperpaths).

A hyperpath (V_P, E_P) 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 to be the set of all molecules linked to the source for a given hypergraph and source set . The characterization of 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 .

The observed trend is consistent in the three cases, as it is shown in Figure ​Figure4.4. For the network composed of all 6542 metabolites and 8392 reactions, there are 23564 pathways connecting 2338 out of 5571 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 (a comparison is provided in supplementary Additional file 1, Figure S1). 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 [28]. In the context of metabolic engineering (i.e. enumerating pathways between native metabolites and heterologous targets) we found only few inconsistent pathways (less than 1% of the 23564 enumerated 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 [1]. 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 5 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 10² ~ 10³-fold slower than FindPath. The execution times are on average of 0.0136 ± 0.0051 s per output with the FindPath algorithm, of 0.218 ± 0.042 s with elementary modes, and of 2.84 ± 1.24 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 1513 pathways, enumerating the pathways with elementary modes takes more than 329.83 seconds, while it can be computed in 20.58 seconds with FindPath. We have measured running times per output and memory usage as function of input size and output size (see Table ​Table11 and Additional file 2, Figure S2). In all cases we found linear growth O(n) 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 3.1 × 10^-2 and 1.6 × 10^-1 seconds, respectively, being ExPA the less efficient code with a constant time of 4.2 × 10^-1 seconds and a linear growth of 1.6 × 10^-7 seconds per input. Similar values were obtained depending on the size of the output, although the scaling in this case for ExPA was of 5.7 × 10^-3 seconds per output. Regarding memory usage, ExPA and FindPath are less demanding, especially for output size, with linear growths of 3.9 and 54 kB per output, respectively, while efmtool required significantly more memory allocation, 3.3 × 10⁴ kB per output.

To the best of our knowledge, the computational complexity of enumerating elementary modes on networks comprising irreversible reactions is up-to-date unknown [28]. In Appendix A.2 we prove that enumerating minimal pathway is an NP-complete problem (reduction to 3-SAT) 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 2338 heterologous targets that are linked to E. coli. As shown in Figure ​Figure5,5, 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 [1] for penicillin N and taxol (see also Additional file 4, Figure S5 and S6). These two examples contained a relative low level of combinatorial complexity, with a maximum number of 14 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 (cf. Figure ​Figure22 in [1] where more than 10000 pathways can be found between tyrosine and chorismate). Further examples are those involved in the biosynthetic pathways for plant steroids leading to brassinolide, where 8 pathways can be used to produce campesterol, one of the initial precursors. The number of pathways grows as we proceed downstream going up to 40 for teasterone, and to 224 for typhasterol. Finally, the number of pathways for the end products, castasterone and brassinolide is of 328. 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 (bootstraps) that are consumed before being produced unless they were added as supplement in the growth medium (cf. in Additional file 3 an example of pathway enumeration with FindPath using bootstrap compounds).

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 (i.e. a hypergraph whose hyperarcs have associated a non-negative real number representing their cost). 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 (see Appendix A.1).

The reason why the algorithms commonly used for graphs cannot be easily adapted for hypergraphs is that given the costs for the hyperpaths leading to the vertices a_i, tails of a hyperarc e_j whose cost is w_j, then the cost of the hyperpath containing all the hyperarcs in _iA_i and e_j is not equal to if the overlap of the hyperpaths A_i is non trivial.

Provided an additive cost function as defined in [23], finding the shortest hyperpath can be done in polynomial time by an algorithm of complexity O(m · (n + ln m)) given in [23], which finds the shortest hyperpath in B-hyperpaths. In fact, this algorithm does apply to minimal hyperpaths given in definition (4).

Generally speaking, the strategy we used in the enumeration algorithm (iterative partition of the space of feasible solutions) 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 (cf. definition in Appendix A.1), since when splitting the problem into a constrained part and a free one (by using the sets R_n and R_f ) we get a problem that is known to be NP-hard (see the proof in [29] where an algorithm is given for the k-shortest hyperpath on hypergraphs whose hyperarcs have only one head vertex).

Let Y(R_f) be the set of all the metabolites produced by reactions in R_f, the mandatory reactions: . Given a hypergraph and the sets R_f, R_n we say that the well-separation condition holds if for every reaction the set Y(r) of products of r either is a subset of Y(R_f) or does not contain elements of Y(R_f). If the well-separation condition holds for a hypergraph with constrained reactions R_f, 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 R_f in a hypergraph whose reactions have one only product, as the hypergraphs defined in [29]. 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 of the enumeration problem have to violate the condition of well-separation for same choice of , 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 6, Figure S9. 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.