Efficiently gap-filling reaction networks

Flux Balance Analysis (FBA) is a method for analyzing genome-scale steady-state models of metabolic networks [1, 2]. Constructing an FBA model often requires completing, or gap-filling, the reaction network. This completion is needed when the model does not grow under some specific growth condition (i.e., given sets of nutrients and secretions). Gap-filling is a computational technique to complete a reaction network based on FBA without referring to the genome. Indeed, a complete knowledge of the functionality of the genome provides a complete reaction network and no gap-filling would be needed. Essentially, given a set of candidate reactions, gap-filling suggests to add some of these reactions to the model so that the FBA model grows, but without guaranteeing that the enzymes for the added reactions exist in the organism. This computational approach to complete a network is viable if it is unknown how to complete the reaction network based on the genome.

To the best of our knowledge, all the current approaches to gap-fill [2–8] use Mixed-Integer Linear Programming (MILP) and compute a minimum set of reactions to add. But, searching for a minimum set of reactions using MILP is computationally expensive if the set of candidate reactions is large. For example, the MetaCyc [9] database contains more than 12,000 reactions. With such a set of candidate reactions, and searching for a minimum set of reactions, MILP typically takes from several minutes to hours, because the number of integer variables increases proportionally to the number of candidate reactions. Frequently, more than an hour is needed for optimally solving such MILPs. Section “Results” presents such concrete cases.

For another example, [8] reports that GapFind with GapFill took about 14 hours when gap-filling a specific FBA model of C. neoformans. GapFind/GapFill uses MILP to find metabolites that cannot be produced and gap-fill the network of reactions of the model to produce them.

Although MILP can find a minimum number of reactions, those reactions are still only plausible suggestions to complete a reaction network. The reactions still need to be analyzed to confirm their validity for the organism. This analysis can be done by searching the scientific literature for references to pathways based on the suggested reactions to add. FastGapFilling also needs such an analysis, but it also provides several sets of reactions from which a broader analysis can be done.

Moreover, these MILP computational techniques focus on finding the smallest set of candidate reactions to complete a model, but that set is not necessarily the most appropriate solution. For example, if only one reaction is needed to produce an essential compound, it might be the case that there is a solution with two reactions that better correspond to the taxonomy of the organism. Knowing the existence of this other non minimum solution might help discover the right missing reactions, which is a solution that might be given by FastGapFilling.

In this paper we present a computational technique using linear programming (LP), avoiding MILP, for reaction gap-filling. Solving a linear program with thousands and even tens of thousands of variables is typically about 1 second. Therefore, LP is a much faster computational technique, compared to MILP, that can give almost immediate feedback on what is needed to complete a reaction network. The gap-filling solutions found by FastGapFilling can be the same or different from those with the MILP approach, as the applications presented in Section “Results” demonstrate. In fact, several sets of candidate reactions are suggested by FastGapFilling, all offering positive growth condition, which may add more insight to complete a reaction network.

The next section presents the FastGapFilling algorithm based on LP.

To the best of our knowledge, all computational techniques published to date for gap-filling a reaction network are based on MILP, including the one used by MetaFlux [6]. MetaFlux can solve FBA models, simulate knockout experiments, and gap-fill the reaction network, the secretions, the nutrients, and the list of biomass metabolites. In MetaFlux, the gap-filling process is based on try-sets and numerical value parameters, called weights, provided by the user. A try-set is a set of candidates to fill the incomplete model. For example, the try-set of reactions is a set, typically large, of candidate reactions for completing the reaction network. A try-set of nutrients would be a set of candidate compounds that could be transported into the organism to enable growth. Each type of candidates (reactions, nutrients, secretions, and biomass metabolites) has its own weight. For each candidate added to the model, its corresponding weight is added to a global objective function and MILP maximizes that objective. Therefore, a positive weight indicates a desire to include as many candidates as possible in the model, and a negative weight indicates a desire to include a minimum number. A typical scenario is assigning a large positive weight for the biomass metabolites and assigning relatively smaller negative weights to all other try-sets. Such a scenario includes as many as possible candidate biomass metabolites by adding a minimum number of reactions, secretions, and nutrients. In MILP, the inclusion/exclusion of each candidate is controlled by an integer variable taking only one of two values, 0 (do not include) or 1 (include). Thousands of such binary variables can exist because — besides the nutrients, secretions, and biomass metabolites — thousands of candidate reactions are typical.

LP and MILP differ in one main respect: all variables in LP can take only fractional values, whereas some variables in MILP must have integer values (for example the binary variables). This difference alone makes MILP computationally more demanding. Indeed, a typical MILP solver must iteratively apply a series of LP problems to find the right integer values to satisfy all constraints. This process can be demanding when thousands or even hundred of thousands of LP cases need solving. In contrast, our technique typically requires solving only a few LP problems.

The FastGapFilling algorithm is presented in Figure 1 and can be summarized in the following way:

1. 1.

   When creating the LP formulation, all candidate reactions M are included with the actual reactions N of the model to complete.

    

The algorithm does not try to gap-fill nutrients or secretions, because the model is assumed to have the correct growth environment specified before attempting to complete the reaction networks. Accordingly, our description of the algorithm does not discuss the growth environment.

The FastGapFilling algorithm accurately represents the model, because the entire biomass reaction is used as if solving the FBA model. The algorithm’s primary use is finding a small set of candidate reactions for adding to the model to produce greater biomass. In particular, if the biomass is zero (i.e., no growth is currently possible), the algorithm tries to obtain growth by activating candidate reactions. That particular case implements reaction gap-filling. Typically, M is a database of reactions such as MetaCyc, which includes approximately 12,000 reactions.

The generic reactions of MetaCyc are instantiated. That is, reactions that operate on classes of compounds are transformed to one or several reactions, where the classes are replaced by specific compounds that exist in the model. Any reversible reaction is transformed into two reactions, one for each direction, a standard practice in constructing FBA models. This transformation simplifies the mathematical modeling as the flux values of all reactions become non-negative.

The algorithm’s output is a list of reactions sets. An implementation may output only the smallest set that would suggest the reactions for model completion. But other implementations may present all sets found and let the user select the most biologically appropriate set. The algorithm does not necessarily find the smallest possible set of candidate reactions to complete the given model. It does try to reduce the size of that set, but it does not compute the absolute minimum set of candidate reactions to complete the model. In other words, FastGapFilling uses an heuristic searching for the smallest set of candidate reactions to complete a model, but it may not find the smallest set.

The search for a small set of candidate reactions to add to the model is based on one weight δ in the objective function. At each iteration of the search, that weight is adjusted either to help the solver increase the biomass flux f _b or to reduce the number of candidate reactions added. The number of iterations done is bounded by ⌈log2|M|⌉.

The costs c _r given as input to the algorithm are typically based on the kind (e.g., enzymatic vs. nonenzymatic) or taxonomic range (e.g., bacteria) of the reactions. The cost of a reaction is proportional to its probability of occurring in a model. For example, a reaction that is not known to occur in the taxonomic range of the model (e.g., a reaction only known in plants, but the model is for a bacteria) should have a higher cost than a reaction known to occur in the taxonomic range of the model. Examples of costs are presented in Section “Results”.

The mathematical formulation of the LP problem of this algorithm follows the standard approach of solving an FBA model, except for the objective function described below and specified in Figure 1. We refer the reader to the papers [2, 4, 6] that describe the generation of such LP. The statement “Generate LP with reactions N and M” in the algorithm assumes that sets of nutrients and secretions are given. Actually, the implementation guides the details of this generation, as it may include many other parameters, such as lower and upper bounds on the flux of some reactions, constraints to avoid ineffective high fluxes, and so on. Essentially, the generated LP must be solvable, must satisfy the steady-state constraints, and must restrict the fluxes to non-negative values for all reactions of N and M, including the biomass reaction.

The main novelty of the LP formulation is the objective function that changes at each iteration of the binary search. The objective function to maximize is $\delta f_b-\sum _{r\in M}{c}_r{f}_r$, where δ is a calculated weight by the binary search, f _b is the flux of the biomass reaction, and the summation of the terms c _r f _r is done over all candidate reactions of M. This function tends to increase the flux of the biomass, as a maximization is applied and the δ f _b is the only positive term in it. This function also tends to decrease the number of candidate reactions that have a non-zero flux (that is, these are active candidate reactions) because the term $-\sum _{r\in M}{c}_r{f}_r$ is negative (note that flux f _r and cost c _r are non-negative values). Indeed, this is the objective of the algorithm: increasing the flux of the biomass while decreasing the number of candidate reactions to add to do so. The essential point is to find the right value for δ such that f _b is non-zero (that is, growth is obtained, and the number of active candidate reactions is small). This is the objective of the binary search done by modifying δ such that growth is maintained and the number of active candidate reactions is reduced. Note that the initial value of δ is |M| (that is, the number of candidate reactions in M). This value enables many candidate reactions to be active simultaneously at the beginning of the binary search.

Note that, in general, this algorithm actually suggests adding candidate reactions to the model to increase the biomass. It could also be used for engineering more efficient metabolic pathways, although this aspect is not covered in this paper.

The main shortcoming of the algorithm is that it does not find a solution when no set of candidate reactions can be found to produce the entire set of biomass metabolites. In such cases, we recommend reducing the number of metabolites in the biomass reaction.

The following section demonstrates the application of FastGapFilling on an E. coli and yeast models.

We have applied the FastGapFilling algorithm to four incomplete E. coli models and to one incomplete yeast model. We started with a working E. coli model for which no added reaction are needed to obtain growth and then removed from 1 to 14 reactions in order to prohibit growth. A similar procedure was done to create the incomplete yeast model. See the following subsections for more details on the creation of these incomplete models.

The original E. coli model is based on the EcoCyc database (version 17.5). The model has 1,763 enzymatic reactions, including the instantiated generic enzymatic reactions. The biomass reaction has 69 metabolites. The growth media is composed of glucose, oxygen, ammonium, phosphate, diphosphate, sulfate, and iron. The secretions are carbon-dioxide, acetate, and water. An upper bound of 16 mmol/g/h was constraining the intake of glucose.

For each incomplete model, we also ran MetaFlux [6] in development mode to gap-fill the reaction network. The development mode of MetaFlux uses MILP with the SCIP [10] linear solver. The MILP gap-fill solutions were compared with the gap-fill solutions found by using the FastGapFilling algorithm.

The relative values of these weights (costs and gain), as opposed to their absolute values, are important. For example, the cost of a spontaneous reaction is smaller than the cost of an enzymatic reaction, because a reaction that does not require an enzyme can occur in any organism. Similarly, the cost of a reaction that is outside the taxonomic range of the organism is set higher than the cost of a reaction in the taxonomic range of the organism because we prefer reactions in the taxonomic range of the organism selected. The candidate reactions are all coming from the MetaCyc database, and their taxonomic ranges are based on the curated pathways of that database.

Notice that the gain of one biomass metabolite is high compared to the cost of any reactions. For example, the ratio between the gain and the cost of a reaction in the taxonomic range of the organism is 100. This means that up to 100 such reactions could be added to the model if one more biomass metabolite can be produced. The intention is to produce a maximum number of biomass metabolites even if a large, but still bounded, number of candidate reactions is needed to produce them. As shown in the following subsections, this general result was obtained for all cases presented because a few reactions produced all metabolites.

FastGapFilling is an efficient technique to gap-fill reaction networks allowing interactive completion of reaction networks of FBA models, whereas MILP fails to offer such an efficiency.

The FastGapFilling algorithm increases the biomass flux by minimizing a weighted minimum sum of fluxes from new reactions, whereas the gap-filling approach using MILP is based on adding a minimum number of new reactions to obtain a non-zero biomass flux.

Using concrete examples, we show that the FastGapFilling algorithm is efficient for performing reaction gapfilling requiring much less computation time than MILP-based approaches. This result is not surprising because MILP may have to solve hundreds or even thousands of LP problems, whereas the FastGapFilling approach requires solving only a few LP problems.