BACKGROUND
genome-wide metabolic reconstructions provide a powerful platform for the analysis of metabolic pathways.
reconstructions for at least sixty-five different species spanning fifty-one different genera, including humans  <cit> , are available today  <cit> .
such reconstructions have been used successfully in several fields of research, including metabolic engineering, evolutionary analysis and metabolic network properties  <cit> .
the mathematical formulation of these reconstructions are called constraint based models.
these models are defined at the core by a sparse stoichiometric matrix, where each column represents a reaction, each row a metabolite, and each entry the corresponding stoichiometric coefficient.
a vast array of computational methods for the analysis of constraint based models is available, and the number continues to grow  <cit> .
perhaps the most widely used of these methods is flux balance analysis .
fba returns a flux distribution through the entire metabolism which optimizes  a given objective function or reaction, such as atp or biomass production.
this flux distribution satisfies a steady-state assumption, such that there is no net production or consumption of any metabolite.
a pre-defined set of constraints, such as upper and lower bounds for each reaction and substrate availability, is also defined  <cit> .
this technique is most commonly applied to metabolic reconstructions of unicellular organisms with the optimization of a theoretical biomass reaction.
this reaction consumes resources such as amino acids, nucleotides and atp at the rate the given organism would need in order to grow and multiply  <cit> .
fba alone has predicted in e. coli the uptake and release rates of certain metabolites  <cit> , cell growth rate under different environmental conditions  <cit>  and gene essentiality  <cit>  with great success.
however, the prediction of internal cell fluxes remains a challenge  <cit> , mainly due to four reasons: the fba solution is not unique.
there are several fluxes within the model that are not well defined, but can exist within certain bounds while the objective function is being optimized  <cit> .
these fluctuations then define an fba solution space, rather than a single, unique output.
organisms might not be operating at maximum capacity  <cit> .
in this case, the objective function might not be fully optimized, but instead be in a near-optimal or sub-optimal state.
furthermore, flux variability analysis shows that the fba solution space increases drastically when considering a near-optimal to optimal state  <cit> , exacerbating the possibility of multiple fba solutions.
the observed metabolic state is also not unique.
that is, a single flux distribution cannot explain all the flux states observed in the experiments, and variation occurs within the bacterial population.
wintermute et al.
<cit>  proposes a cloud theory for metabolic regulation, where bacteria are allowed to fluctuate within a near-optimal solution space.
the study also demonstrates that the variability of fluxes within this region matches the observed variability within the data.
thermodynamically infeasible loops can appear in the fba output.
these are sets of pathways termed “type iii pathways”  <cit> , which are composed of internal reactions that can be combined to yield a steady-state with no net input or output.
these cycles can clutter the fba output, hindering any subsequent analysis  <cit> .
to better predict internal fluxes, some studies have relied on additional thermodynamic constraints  <cit> .
several of these studies successfully reduce the fba solution space, but still leave a variety of states that, although thermodynamically feasible, are physiologically improbable  <cit> .
other studies rely on the optimization of a thermodynamic cost, and may miss other key biological factors that govern metabolism  <cit> .
significant efforts have also been directed toward integrating omics data and fba  <cit> , but the generally low correlation between gene expression and the associated reaction flux makes this a challenging field  <cit> .
furthermore, the predictive power of these integrated omics-fba models is limited, since experimental data is needed to predict the fluxes.
more recent studies have relied on the theory that e. coli chooses its pathways based on the minimization of protein cost  <cit>  and number of metabolic steps  <cit> .
several successful methods based on these assumptions have been proposed, which include minimization of net metabolic flux  <cit> , minimization of the number of steps in the metabolism  <cit>  and enzymatic level constraints  <cit> .
one alternative method utilizes elementary modes to pre-define the directionality of reactions, thus reducing the fba solution space  <cit> .
these methods have been very successful at improving predictions of growth rate, substrate usage and internal cell fluxes in unicellular organisms.
several studies have suggested, however, that unicellular organisms in reality grow at rates lower than those predicted by fba alone  <cit> .
although some of the above-mentioned methods do predict growth rates lower than the standard fba, they all rely on implementing additional constraints upon optimization of the objective function.
furthermore, although this approach has been successfully implemented in the study of cancer cells  <cit> , this objective function will most likely not hold true for the analysis of healthy multicellular organisms, like mammals, as healthy cells in these systems have more complex objectives than to simply grow and multiply.
additionally, these one-step optimization methods with additional constraints return a single flux distribution, and they are unable to explore the near-optimal solution space.
this limitation is significant given a recent study has suggested the flux distribution of e. coli can vary freely within a near-optimal space  <cit> .
therefore, in order to further improve fba predictions, especially as the field expands to include multicellular organisms, new techniques which explore the sub-optimal solution space need to be developed.
to address this need, we propose a two-step optimization fba method for predicting internal fluxes, termed cost reduced sub-optimal fba , which is suitable for sub-optimal objectives.
rather than imposing additional constraints when optimizing the objective function, we fix this objective at a predefined value.
an estimated protein cost throughout the metabolism is then minimized in order to predict the internal cell fluxes.
by varying the biomass production value, we are able to predict how changes in pathway usage depend on changes in other conditions, particularly glucose availability.
although this method is not well suited for predicting growth rates, we are able to predict metabolic flux distributions within a near-optimal solution space.
furthermore, by decomposing the model using an adapted version of the extreme pathway analysis, we are able to break down the energy production pathways and further understand the model behavior as glucose concentrations change.
we validate our methods using e. coli as the model organism.
methods
cost calculation and implementation
some of the methods most successful at predicting internal cell fluxes have been based on the theory that e. coli is subject to a protein cost constraint  <cit> .
five of such methods available in the literature, along with corsofba, are summarized in table  <dig>  here we explore a similar principle, and calculate the protein cost of a reaction based on the net flux through the reaction , the enzyme molecular weight , and a thermodynamic penalty for reversible reactions.
the cost term used in corsofba is defined as:  j·mw·expα·Δrg′or·t table  <dig>
optimization methods with additional constraints
method
enzymatic cost
constraint
∑aiji
∑gi·mwi
∑mi+δ·∑gi
∑mwijikcati
∑ji++keqi·ji−
∑ji·mwi·expα·Δrg′ior·t
variables
k
k
these methods calculate an overall cost through the metabolism by summing individual costs associated with each reaction.
while previous methods have been used to constrain the fba solution, corsofba optimizes this cost for any given growth rate.
where Δrg′o is the standard gibbs free energy of the given reaction, r is the ideal gas constant and t is the temperature.
the molecular weight term is included to represent the amount of resources needed in order to produce a sufficient amount of enzymes to maintain the associated reaction flux.
the thermodynamic penalty, applied only to reversible reactions, represents the cost associated with the change in the concentration of metabolites required for the reaction to flow in the desired direction.
that is, a reaction with positive Δrg′o would require a large disparity in the concentration of the associated metabolites, so the cost term is multiplied by a positive number larger than one.
on the other hand, if Δrg′o is negative, and the reaction is therefore favorable, the cost term is multiplied by a positive number between zero and one, given that little or no metabolic adjustment is necessary.
it is worth noting that this thermodynamic cost does not directly represent a change in enzyme concentration, but rather an increase or decrease in the overall cost associated with the reaction based on its standard gibbs free energy.
furthermore, this term is applied separately to forward and backward directions of reversible reactions, favoring the direction with negative standard gibbs free energy.
a similar thermodynamic cost, used without molecular weight considerations, has been applied by holzhütter  <cit> , weighing only the fluxes directed in the thermodynamically unfavorable direction.
finally, the parameter α is chosen to be  <dig> rtmolkcal, making the final thermodynamic cost simply:  exp <dig> ·Δrg′omolkcal.
this value is chosen in order to balance the contribution of molecular weight and thermodynamic penalty, yielding a prediction considering both costs .
without this parameter, due to the exponential nature of the thermodynamic penalty, this term would quickly approach extremely large numbers  or zero , and only the thermodynamic cost would essentially be considered.
in that case, reaction directionality would essentially be predetermined based solely on the standard gibbs free energy, which could be problematic.
the gibbs free energy of reactions was estimated from the metacyc database  <cit>  as the change in gibbs energy of formation between reaction compounds.
molecular weight values were obtained for  <dig> of the  <dig>  enzyme catalyzed reactions in the e. coli ijr <dig> reconstruction  <cit> .
these values were extracted from the brenda database  <cit>  whenever available.
if e. coli measurements were not available, we extracted the values for the organism most closely related to e. coli.
a small number of molecular weights was estimated from the ecocyc database  <cit> , and the remaining  <dig> values were defined as the median of the calculated values.
additional information concerning the standard gibbs free energy of reactions, protein molecular weights and protein cost calculations can be found in the supplemental information .
since the protein cost term scales linearly with the flux through the associated reaction, we incorporated the molecular weight and thermodynamic cost directly into the model.
all enzyme associated reactions were split into forward and backward parts, and the cost was added as a produced metabolite:  initialreactionrxn:a+b⇔c+dreactionswithcostrxnf:a+b⇒c+d+mwrxn·expα·Δrg′rxnor·trxnb:c+d⇒a+b+mwrxn·exp−α·Δrg′rxnor·t
a reaction which consumes this produced metabolite was then added to the model.
while every original reaction in the model produces this cost at a rate relative to its absolute flux, this added reaction becomes the only one to consume this cost.
minimizing the flux through this reaction will then return the flux distribution with the lowest enzyme associated cost to achieve the predefined biomass production value.
reactions with no associated enzyme were given a cost of zero.
these methods were applied to the e. coli ijr <dig> reconstruction  <cit> .
during all simulations the glucose uptake bound was set to an arbitrary value.
uptake bounds for co <dig>  fe <dig>  h2o, h, na, nh <dig>  o <dig>  phosphate, and so <dig> were set to arbitrarily large values, much larger than the glucose uptake bound, and were considered to be present in excess.
all exchange reactions release bounds were set to arbitrarily large values except for glucose, formate and o <dig>  which were set to zero.
predicted flux values were then normalized by either the flux of glucose to glucose-6-phosphate or overall flux through enzyme associated reactions.
it is worth noting that this protein cost term does not account for the enzyme turnover rates or substrate affinity.
we have opted for this cost term since the kinetic parameters of enzymes are extremely difficult to obtain experimentally, and vary greatly based on ph and temperature values.
furthermore, the simplification of a constant turnover rate and substrate affinity for all enzymes has been shown to yield significant results  <cit> .
fundamental pathways analysis
different methods have been proposed to characterize and decompose the fba solutions space  <cit> , including elementary modes  <cit> , extreme pathways  <cit>  and minimal generators  <cit> .
however, these tools have mostly been applied to small networks.
the number of pathways obtained for each of these methods, as well as the computational cost of these analyses, increase exponentially with the size of the reconstruction  <cit> .
with large-scale metabolic reconstructions, these decompositions quickly become impractical tools.
here we propose a modified version of extreme pathways which yields a significantly lower number of pathways during the decomposition.
each of these pathways is associated with a protein cost and an atp production potential, which helps elucidate the model behavior as the value of the objective function varies.
we find that this large number of pathways stems partially from the need to balance each pathway to a steady-state, such that no metabolite has a net production or consumption rate.
certain metabolites, however, known as currency metabolites, are responsible for basic cell functions and occur in numerous reactions in the model.
standard reactions involving these metabolites, such as atp hydrolysis, can be obtained through different combinations of other model reactions.
these combinations, or loops, can then be used multiple times with the same purpose, increasing the number of pathways which have essentially the same core reactions.
this phenomenon can be illustrated by extreme pathways analysis applied to the pathway and set of reactions presented in figure  <dig>  the pathway depicted in figure 1a is the production of energy through the conversion of glucose into lactate, which produces two mols of atp for every mol of glucose.
reaction names are extracted from the e. coli reconstruction by orth et al.
<cit> .
figure 1b illustrates loops and reactions present in the same reconstruction that hydrolyze atp, any of which can be combined with the reactions in figure 1a to produce a pathway satisfying the atp steady-state condition.
furthermore, atps4r also introduces a proton imbalance when hydrolyzing atp.
this imbalance can then be corrected by any of the proton transport loops also presented in figure 1b.
the combination of this pathway with the given loops then give us multiple extreme pathways with essentially the same core reactions; and they represent the same phenomenon - the conversion of glucose to lactate through glycolysis.figure  <dig>
metabolic loops increase pathway multiplicity.
pathway converting  <dig> mol of glucose to  <dig> mols of lactate and producing  <dig> mols of atp.
reaction fluxes are represented in parentheses.
extracellular metabolites are defined by .
all other metabolites are in the cytosolic space.
network loops and reactions that can hydrolyze atp and transport protons across the cell membrane.
these reactions and loops can be used to balance the reactions in , yielding multiple balanced pathways with the same set of backbone reactions
in order to extract only the core reactions in each pathway, we then remove a selected set of currency metabolites from the reconstruction.
this removal reduces the number of pathways we obtain from extreme pathway analysis.
this metabolite removal was performed under certain considerations: metabolites removed should be completely deleted from all reactions and compartments in the model.
removing a metabolite from one reaction and not others could block certain reactions, since fewer options to balance that metabolite would be present in the reconstruction.
no reaction should become a sink or source of metabolites.
allowing reactions to become sinks or sources would introduce exchange reactions for certain metabolites, affecting the model decomposition.
reactions removed from the reconstruction should not contain any metabolite which has not been removed from the reconstruction.
the removal of these reactions would remove any pathway which uses these reactions to balance the metabolite still in the model.
while we should be aware of these considerations when modifying the model, breaking any of these specifications should not invalidate the fundamental pathway analysis.
the model decomposition could still be performed, but some limitations might have to be considered.
we have illustrated this fundamental pathways method using the e. coli model by orth et al.
<cit> , and have demonstrated that the original pathways can be recovered by adding back the identified loops.
this analysis can be found in the supplemental information .
while we obtain fewer pathways in this tailored reconstruction, each set of reactions is associated with an imbalance in the original model.
for example, instead of multiple extreme pathways, the set of reactions in figure 1a would yield a single fundamental pathway after the removal of currency metabolites from the reconstruction.
these reactions alone, however, would yield atp, adp, phosphate, proton and water imbalances in the original model.
in the fundamental pathways analysis, each resulting fundamental pathway is then associated with this imbalance, which can be used to further analyze each pathway.
in this study, this imbalance is used to analyze the atp production capacity of each pathway obtained.
RESULTS
comparison between different cost functions
cost-optimal flux simulations were performed using four different types of cost:  molecular weights alone,  thermodynamic penalties alone,  the previously mentioned combination of molecular weights and thermodynamic penalties and  uniform costs.
the use of uniform costs gives a simple minimization of the overall flux through the enzyme associated reactions.
this minimization strategy is a modified version of a widely regarded two-step optimization method called pfba  <cit> .
pfba has been previously used to predict unique flux distributions at the predicted fba optima.
simulations comparing these four costs have been performed for growth rates ranging from  <dig> to 100% of the predicted optima.
flux distributions were normalized by the flux of glucose to glucose-6-phosphate, and results for several central carbon metabolism reactions can be found in the supplemental information .
while simulations using molecular weights, thermodynamic costs, or a combination of both present considerably different flux distributions at different sub-optimal values, simulations using uniform costs quickly converge to a unique flux distribution as the objective value is lowered.
the same set of simulations was also compared to several data points from two metabolic flux analysis  experiments: ishii et al.
<cit>  and yao et al.
<cit> .
correlation and sum of squared error between all simulated and experimental flux distributions for several central carbon metabolism reactions were calculated.
these results are presented in the supplemental information .
a similar method to that used by holzhütter  <cit>  was also included in additional file 8: figure s <dig>  while pfba simulations yield the highest correlation values in the suboptimal space when compared to the ishii et al.
dataset, these same simulations also yield the highest sum of squared error.
furthermore, when the objective function was lowered from 100% down to 65% of optima, corsofba yields the lowest sum of squared error for all data points except one , and the highest correlation values for the higher dilution rates in the yao et al.
data .
comparison between metabolic flux analysis and simulated fluxes
in order to visualize the trend of how different fluxes change in the sub-optimal space, simulations utilizing molecular weight costs, thermodynamic costs and a combination of both, normalized again by the flux of glucose to glucose-6-phosphate, were plotted alongside the experimental values from the two mfa experiments.
fluxes were plotted starting at 100% down to 75% of the predicted optima, and the results for several central carbon metabolism reactions are presented in figure  <dig>  the threshold of 75% was chosen in order to match the trend observed in our simulations to the experimental data  <cit> .figure  <dig>
comparison between simulations and mfa experiments.
comparison between selected simulated fluxes using molecular weights, thermodynamic penalties, and a combination of the two costs, and metabolic flux analysis  data.
fluxes are normalized by glucose to glucose-6-phosphate conversion rates.
x-axes values are kept constant for all plots.
reaction names are taken from the ecoli ijr <dig> model  <cit> .
yellow boxes indicate pentose phosphate pathway , blue indicates glycolysis, green indicates tca cycle and red indicates other reactions
as the objective function value decreases, corsofba flux distributions show close agreement with mfa data for higher growth rates , particularly for simulations considering both thermodynamic and molecular weight costs.
fluxes through the pentose phosphate pathway  and glycolysis remain relatively constant, while the flux through the tca cycle gradually decreases.
the decrease in tca cycle usage leads to higher levels of acetate release.
quantitatively, the fluxes through glycolysis and the tca cycle are slightly under-predicted by our simulations, while the flux through the ppp is generally over-predicted.
similar patterns have been observed in other fba method approaches  <cit> , and this may arise due to the optimization of the biomass function alone.
the optimization of atp production alongside biomass has been shown to yield better predictions  <cit> .
atp production optimization could also lead to higher fluxes through the tca cycle and lower ppp fluxes.
it is also worth noting that, although fluxes through the tca cycle are under-predicted according to figure  <dig>  mfa experimental results can be inconsistent, and our predicted tca fluxes show good agreement with two other mfa studies  <cit> .
these simulations mirror a behavior known as overflow metabolism  <cit> , where e. coli, at high growth rates, moves away from the full oxidation of glucose through the tca cycle and uses less atp efficient pathways, releasing acetate instead of co  <dig>  our simulations support the hypothesis that this behavior stems from a tradeoff between enzyme cost and energy yield  <cit> .
that is, when more glucose is available, e. coli uses a higher flux through pathways that are less enzymatically costly, but which produce fewer atp per mmol of glucose.
predicted fluxes using molecular weights only and a combined cost also found excellent agreement with the experimental values in the reactions ppc  and icl  <cit> .
both experiments show a decrease in icl flux and an increase in ppc flux, from near-optimal conditions to a growth rate of  <dig>  h- <dig>  this result is further supported by the experimental results from nanchen et al.
<cit> , that found a lower flux through icl during a growth rate of approximately  <dig>  h- <dig> when compared to a growth rate of  <dig>  h- <dig>  suggesting a transient response through this pathway as glucose concentration increases.
the same transient response was found in our simulations.
these results suggest that e. coli utilizes these anaplerotic reactions  to relieve enzymatic cost, and that considering these costs during fba may increase internal flux predictions when comparing to the minimization of metabolic steps alone.
simulation results also suggest that the optimization of the biomass function may yield fluxes that are fundamentally different from those at near-optimal conditions.
at optimal growth, there is no predicted flux through the glyoxylate shunt, considerably lower fluxes through the tca cycle and ppp, and a flux through the ppc reaction not present at near-optimal conditions.
although comparisons between experimental data and simulated flux distributions show that the highest correlation and lowest error are found at optimal growth, our simulations also indicate that the implementation of these costs yield comparable results in the sub-optimal space .
moreover, the suboptimal corsofba approach can better predict fluxes through certain reactions, such as ppc and icl.
while previous studies have shown that the fba solution space increases drastically when considering the objective function at near-optimal to optimal values  <cit> , here we show that the optimization of the enzymatic cost at near-optimal conditions yields results that are more consistent with experimental data for certain reactions.
this result strongly suggests that exploring the fba solution space at near-optimal to optimal objective values may increase the predictive accuracy of internal cell fluxes.
the near optimal solution space was also analyzed for knockout strains.
flux distributions for six e. coli knockout strains reported by ishii et al.
<cit>  were compared to cost-optimal simulations using the combined cost function.
results show that cost-optimal flux distribution from 95% to optimal yield similar or better predictions according to both correlation and sum of square error .
predictions using uniform costs quickly diverge from experimental data.
these simulations further support the “cloud theory” proposed by wintermute et al.
<cit> , and suggest that the experimental data can also be in good agreement with near-optimal flux distributions.
comparison between gene expression data and simulated fluxes
to further validate this analysis, simulations using the combined cost function were also compared to gene expression values reported in the mfa studies  <cit> .
while the overall correlation between reaction flux and gene expression is moderate at best  <cit> , increased expression of all genes participating in a particular pathway can be considered a good indication of increased flux  <cit> .
simulations were performed under the same conditions as before, utilizing the combination of molecular weights and thermodynamic penalty, but this time the flux distribution was normalized by the overall flux through the enzyme associated reactions.
comparisons between simulated values and gene expression data are shown in figure  <dig> figure  <dig>
comparison between simulations and gene expression data.
comparison between selected simulated fluxes, using a combination of molecular weights and thermodynamic penalties, and associated gene expression.
plot axes are the same as the ones defined in the ptsg plot unless otherwise specified.
multiple genes associated with the same reaction are included in the same box.
reaction names are taken from the e. coli ijr <dig> reconstruction  <cit> .
yellow boxes indicate pentose phosphate pathway, blue indicates glycolysis, green indicates tca cycle and red indicates other reactions
although gene expression data can also be inconsistent between experiments, and some reactions are associated with multiple genes, this comparison further supports the qualitative results presented in the previous section.
when normalizing the predicted fluxes by the overall flux through the enzyme associated reactions, an increased relative uptake of glucose is observed as the objective value decreases.
as a result, an increased relative flux through glycolysis is also predicted.
furthermore, the relative flux through the tca cycle still decreases as the growth rate increases.
the increase in relative glucose uptake and flux through glycolysis, as well as the decreased flux through the tca cycle, are supported by the relative gene expression associated with these pathways  <cit> .
simulated values also find good agreement with fluxes through the entner-doudoroff  pathway.
genes edd and eda, associated with this pathway, have long been known to exist in e. coli  <cit> , but the activity of their associated reactions has been observed mainly under growth on gluconate, glucuronate, and methyl-beta-d-glucuronide; phosphate limitation; and carbon starvation  <cit> .
due to this fact, these reactions are generally not included in mfa experiment networks.
in contrast with murray et al., our simulations using molecular weight and combined costs predict the use of the ed pathway under excess glucose conditions, not starvation .
these predictions are supported both by the genetic data presented by the mfa studies  <cit> , which show an increase in edd and eda expression at high glucose concentrations, and mfa experimental results by harcombe et al.
<cit> , where fluxes through the ed pathway were measured in bacteria growing at growth rates above  <dig> h- <dig>
fundamental pathways analysis
to better understand the transition between metabolic states as the value of the objective function is decreased, the energy producing pathways of the e. coli ijr <dig> reconstruction were decomposed using the fundamental pathways analysis described in the methods section.
details on how these reactions were calculated can be found in the supplemental information .
briefly, starting with the atp imbalance associated with each fundamental pathway, the total atp potential of each pathway was estimated based on the imbalance of other energy generating metabolites, such as nadh and ubiquinol- <dig>  the lowest cost of converting these metabolites into atp was then added to the total cost associated with the pathway.
the final atp production by enzymatic cost was then compared to the potential atp production by mmol of glucose .figure  <dig>
fundamental pathway analysis.
atp production potential of energy producing fundamental pathways, compared to the associated protein cost of each mmol of atp.
optimal  and near optimal  pathways are highlighted in the scatter plot.
ops are plotted in blue and nops in red in the central carbon metabolism diagrams
in accordance with general e. coli metabolism knowledge, the fundamental pathway analysis found two optimal pathways , according to atp production by both mmol of glucose and protein cost.
the most energy productive pathway is the full oxidation of glucose through the tca cycle .
the most cost efficient pathway was found to be the production of acetate through glycolysis .
although this pathway potentially produces less than half the amount of atp as op <dig>  the total enzymatic cost per mmol of atp is lower.
also in agreement with general knowledge, this analysis predicted the release of acetate to be more efficient than the release of both lactate and ethanol.
in the simulations presented in previous sections, a gradual transition from op <dig> to op <dig> is observed as the objective function value decreases.
this trend also supports the idea that overflow metabolism takes place in order to alleviate enzymatic cost.
under low glucose conditions, e. coli would need to extract from glucose the largest possible amount of energy in order to sustain growth.
as more glucose becomes available, however, this bacteria can afford to consume more substrate and produce energy through cheaper pathways.
one interesting observation from this analysis is the existence of near optimal pathways , which combine the previously described optimal pathways with the ppp or ed pathways.
nop <dig> and nop <dig> both produce more atp per mmol of glucose than op <dig> but less than op <dig>  while remaining less cost efficient than op <dig> but slightly more so than op <dig>  nop <dig> and nop <dig> reduce the flux through upper glycolysis by partially re-routing the flux through the ppp, and they exhibit no flux through phosphoglucose isomerase .
a third near optimal pathway, termed nop <dig>  skips the upper glycolysis through the ed pathway.
this pathway is both less cost and energy effective than op <dig>
the use of these near optimal pathways could explain the inconsistent ppp fluxes reported by the mfa experiments considered here.
a short recovery in the simulated ppp flux near a growth rate of  <dig>  h − <dig>  demonstrates that this pathway is a viable option at certain dilution rates, due to its coupling with the ppp and biomass production.
it has also been demonstrated that Δpfka deficient e. coli strains reduce flux through phosphofructokinase  <dig> by diverting fluxes through the ppp  <cit> , hence using nop <dig> and nop <dig> instead of op <dig> or op <dig> when the cost through upper glycolysis is increased.
this fundamental pathway analysis elucidates how we would expect the model to behave in terms of energy generation as we optimize the protein cost at sub-optimal states.
based on this analysis, in order to fulfill its catabolic needs, the model transitions from op <dig> to op <dig> as we move away from optimal growth.
although this is in fact the general trend we observe, anaplerotic needs also need to be considered.
these needs are addressed most efficiently by reactions not in op <dig> or op <dig>  such as pep carboxylase and the glyoxylate shunt.
interestingly, the model also predicts the use of the ed pathway with anaplerotic purposes at high growth rates.
while this pathway has been studied mostly for its catabolic activity in z. mobilis and several pseudomonas species, simulations here predict this pathway to be the cheapest way to get glucose shuttled to the tca cycle for the production of building blocks, while giving up little efficiency in energy pathways through nop <dig>
CONCLUSIONS
here we present corsofba, a method for analyzing and predicting metabolic flux distributions at sub-optimal states based on protein and thermodynamic cost minimization.
corsofba demonstrates that the optimization of protein cost at near-optimal states can produce significantly different results from those at the optimal objective.
furthermore, although correlation and error calculations indicate better predictions at the optimal solution, protein-optimal results at sub-optimal objectives show better agreement with mfa experiments and gene expression profiles for anaplerotic reactions.
these results suggest it is important to explore the sub-optimal fba space in order to better predict internal cell fluxes.
although the method described here is not suited for predicting growth rates, it provides a platform for analyzing internal cell fluxes in the sub-optimal space.
we believe corsofba will be particularly useful given recent studies have suggested e. coli can exist freely within a near-optimal space  <cit> .
we also believe this method can be useful in predicting fluxes in healthy, multicellular organisms, which have more complex objectives than the production of biomass.
future studies are merited to explore whether the predictive power of methods which optimize the objective function under enzymatic level constraints, such as moment  <cit>  and fbawmc  <cit> , would also benefit from exploring the sub-optimal solution space.
our results suggest that significant conclusions could be drawn by adapting these methods to optimize the enzymatic cost they implement in a near-optimal space, rather than using these simply as a model constraint.
additional files
additional file  <dig>
figure s <dig>  sensitivity analysis of parameter α. model simulations ranging from  <dig> to 75% of optima with varying values of α. simulations are performed using both thermodynamic costs only and combined costs.
additional file  <dig>
gibbs free energy change distribution.
file showing plot and distribution of the gibbs free energy change of reactions considered.
additional file  <dig>
molecular weight and cost calculation.
details on how molecular weight values were obtained and final costs calculated.
molecular weights and final costs distribution are also plotted.
additional file  <dig>
table s <dig>  thermodynamic, molecular weight and final cost values.
table includes all values used during simulations: gibbs energy of formation for all metabolites, calculated gibbs energy change for all reactions based on metabolite energy of formation and molecular weights obtained from databases.
additional file  <dig>
fundamental pathway analysis.
example demonstrating the recovery of extreme pathways by adding reactions loops back into the fundamental pathways.
details on how the fundamental pathway analysis presented in the manuscript was performed are also available.
additional file  <dig>
table s <dig>  model decomposition.
extreme pathways and fundamental pathways calculated both for the analysis presented in the manuscript and the pathway recovery analysis in additional file  <dig>
additional file  <dig>
figure s <dig>  cost values comparison.
cost-optimal simulations using molecular weights alone, thermodynamic penalties alone, a combination of the two and uniform costs are presented for major central carbon metabolism reactions in the range of  <dig> to 50% of maximal growth.
additional file  <dig>
figure s <dig>  error and correlation calculations.
correlation and sum of squared error are calculated between all simulated and experimental flux distributions for major central carbon metabolism reactions.
additional file  <dig>
figure s <dig>  error and correlation calculations for knockout strains.
correlation and sum of squared error are calculated between mfa experimental data and simulated flux distributions using the combined cost function for six knockout strains.
abbreviations
tcatricarboxylic acid
fbaflux balance analysis
ppppentose phosphate pathway
mfametabolic flux analysis
edentner-doudoroff
glcglucose
g6pglucose-6-phosphate
f6pfructose-6-phosphate
fdpfructose- <dig> -biphosphate
dhapdihydroxyacetone phosphate
g3pglyceraldehyde 3-phosphate
13dpg3-phospho-d-glyceroyl phosphate
3pg3-phospho-d-glycerate
2pgd-glycerate 2-phosphate
pepphosphoenolpyruvate
pyrpyruvate
accoaacetyl-coa
actpacetyl phosphate
acacetate
citcitrate
icitisocitrate
akg2-oxoglutarate
succoasuccinyl-coa
succsuccinate
fumfumarate
mal-ll-malate
oaaoxaloacetate
6pgl6-phospho-d-glucono- <dig> -lactone
6pgc6-phospho-d-gluconate
ru5p-dd-ribulose 5-phosphate
r5palpha-d-ribose 5-phosphate
xu5p-dd-xylulose 5-phosphate
s7psedoheptulose 7-phosphate
e4pd-erythrose 4-phosphate
2ddg6p2-dehydro-3-deoxy-d-gluconate 6-phosphate
competing interests
the authors declare that they have no competing interests.
authors’ contributions
as and aq conceived the methods and revised the manuscript.
as drafted the manuscript and generated the results.
both authors have read the manuscript and approved the final version.

 genome-wide metabolic reconstructions provide a powerful platform for the analysis of metabolic pathways. reconstructions for at least sixty-five different species spanning fifty-one different genera, including humans, are available today. such reconstructions have been used successfully in several fields of research, including metabolic engineering, evolutionary analysis and metabolic network properties. the mathematical formulation of these reconstructions are called constraint based models. these models are defined at the core by a sparse stoichiometric matrix, where each column represents a reaction, each row a metabolite, and each entry the corresponding stoichiometric coefficient. a vast array of computational methods for the analysis of constraint based models is available, and the number continues to grow. perhaps the most widely used of these methods is flux balance analysis. fba returns a flux distribution through the entire metabolism which optimizes a given objective function or reaction, such as atp or biomass production. this flux distribution satisfies a steady-state assumption, such that there is no net production or consumption of any metabolite. a pre-defined set of constraints, such as upper and lower bounds for each reaction and substrate availability, is also defined. this technique is most commonly applied to metabolic reconstructions of unicellular organisms with the optimization of a theoretical biomass reaction. this reaction consumes resources such as amino acids, nucleotides and atp at the rate the given organism would need in order to grow and multiply. fba alone has predicted in e. coli the uptake and release rates of certain metabolites, cell growth rate under different environmental conditions and gene essentiality with great success. however, the prediction of internal cell fluxes remains a challenge, mainly due to four reasons: the fba solution is not unique. there are several fluxes within the model that are not well defined, but can exist within certain bounds while the objective function is being optimized. these fluctuations then define an fba solution space, rather than a single, unique output. organisms might not be operating at maximum capacity. in this case, the objective function might not be fully optimized, but instead be in a near-optimal or sub-optimal state. furthermore, flux variability analysis shows that the fba solution space increases drastically when considering a near-optimal to optimal state, exacerbating the possibility of multiple fba solutions. the observed metabolic state is also not unique. that is, a single flux distribution cannot explain all the flux states observed in the experiments, and variation occurs within the bacterial population. wintermute et al. proposes a cloud theory for metabolic regulation, where bacteria are allowed to fluctuate within a near-optimal solution space. the study also demonstrates that the variability of fluxes within this region matches the observed variability within the data. thermodynamically infeasible loops can appear in the fba output. these are sets of pathways termed “type iii pathways”, which are composed of internal reactions that can be combined to yield a steady-state with no net input or output. these cycles can clutter the fba output, hindering any subsequent analysis. to better predict internal fluxes, some studies have relied on additional thermodynamic constraints. several of these studies successfully reduce the fba solution space, but still leave a variety of states that, although thermodynamically feasible, are physiologically improbable. other studies rely on the optimization of a thermodynamic cost, and may miss other key biological factors that govern metabolism. significant efforts have also been directed toward integrating omics data and fba, but the generally low correlation between gene expression and the associated reaction flux makes this a challenging field. furthermore, the predictive power of these integrated omics-fba models is limited, since experimental data is needed to predict the fluxes. more recent studies have relied on the theory that e. coli chooses its pathways based on the minimization of protein cost and number of metabolic steps. several successful methods based on these assumptions have been proposed, which include minimization of net metabolic flux, minimization of the number of steps in the metabolism and enzymatic level constraints. one alternative method utilizes elementary modes to pre-define the directionality of reactions, thus reducing the fba solution space. these methods have been very successful at improving predictions of growth rate, substrate usage and internal cell fluxes in unicellular organisms. several studies have suggested, however, that unicellular organisms in reality grow at rates lower than those predicted by fba alone. although some of the above-mentioned methods do predict growth rates lower than the standard fba, they all rely on implementing additional constraints upon optimization of the objective function. furthermore, although this approach has been successfully implemented in the study of cancer cells, this objective function will most likely not hold true for the analysis of healthy multicellular organisms, like mammals, as healthy cells in these systems have more complex objectives than to simply grow and multiply. additionally, these one-step optimization methods with additional constraints return a single flux distribution, and they are unable to explore the near-optimal solution space. this limitation is significant given a recent study has suggested the flux distribution of e. coli can vary freely within a near-optimal space. therefore, in order to further improve fba predictions, especially as the field expands to include multicellular organisms, new techniques which explore the sub-optimal solution space need to be developed. to address this need, we propose a two-step optimization fba method for predicting internal fluxes, termed cost reduced sub-optimal fba, which is suitable for sub-optimal objectives. rather than imposing additional constraints when optimizing the objective function, we fix this objective at a predefined value. an estimated protein cost throughout the metabolism is then minimized in order to predict the internal cell fluxes. by varying the biomass production value, we are able to predict how changes in pathway usage depend on changes in other conditions, particularly glucose availability. although this method is not well suited for predicting growth rates, we are able to predict metabolic flux distributions within a near-optimal solution space. furthermore, by decomposing the model using an adapted version of the extreme pathway analysis, we are able to break down the energy production pathways and further understand the model behavior as glucose concentrations change. we validate our methods using e. coli as the model organism. methods cost calculation and implementation some of the methods most successful at predicting internal cell fluxes have been based on the theory that e. coli is subject to a protein cost constraint. five of such methods available in the literature, along with corsofba, are summarized in table here we explore a similar principle, and calculate the protein cost of a reaction based on the net flux through the reaction, the enzyme molecular weight, and a thermodynamic penalty for reversible reactions. the cost term used in corsofba is defined as: j·mw·expα·Δrg′or·t table optimization methods with additional constraints method enzymatic cost constraint ∑aiji ∑gi·mwi ∑mi+δ·∑gi ∑mwijikcati ∑ji++keqi·ji− ∑ji·mwi·expα·Δrg′ior·t variables k k these methods calculate an overall cost through the metabolism by summing individual costs associated with each reaction. while previous methods have been used to constrain the fba solution, corsofba optimizes this cost for any given growth rate. where Δrg′o is the standard gibbs free energy of the given reaction, r is the ideal gas constant and t is the temperature. the molecular weight term is included to represent the amount of resources needed in order to produce a sufficient amount of enzymes to maintain the associated reaction flux. the thermodynamic penalty, applied only to reversible reactions, represents the cost associated with the change in the concentration of metabolites required for the reaction to flow in the desired direction. that is, a reaction with positive Δrg′o would require a large disparity in the concentration of the associated metabolites, so the cost term is multiplied by a positive number larger than one. on the other hand, if Δrg′o is negative, and the reaction is therefore favorable, the cost term is multiplied by a positive number between zero and one, given that little or no metabolic adjustment is necessary. it is worth noting that this thermodynamic cost does not directly represent a change in enzyme concentration, but rather an increase or decrease in the overall cost associated with the reaction based on its standard gibbs free energy. furthermore, this term is applied separately to forward and backward directions of reversible reactions, favoring the direction with negative standard gibbs free energy. a similar thermodynamic cost, used without molecular weight considerations, has been applied by holzhütter, weighing only the fluxes directed in the thermodynamically unfavorable direction. finally, the parameter α is chosen to be rtmolkcal, making the final thermodynamic cost simply: exp ·Δrg′omolkcal. this value is chosen in order to balance the contribution of molecular weight and thermodynamic penalty, yielding a prediction considering both costs. without this parameter, due to the exponential nature of the thermodynamic penalty, this term would quickly approach extremely large numbers or zero, and only the thermodynamic cost would essentially be considered. in that case, reaction directionality would essentially be predetermined based solely on the standard gibbs free energy, which could be problematic. the gibbs free energy of reactions was estimated from the metacyc database as the change in gibbs energy of formation between reaction compounds. molecular weight values were obtained for of the enzyme catalyzed reactions in the e. coli ijr reconstruction. these values were extracted from the brenda database whenever available. if e. coli measurements were not available, we extracted the values for the organism most closely related to e. coli. a small number of molecular weights was estimated from the ecocyc database, and the remaining values were defined as the median of the calculated values. additional information concerning the standard gibbs free energy of reactions, protein molecular weights and protein cost calculations can be found in the supplemental information. since the protein cost term scales linearly with the flux through the associated reaction, we incorporated the molecular weight and thermodynamic cost directly into the model. all enzyme associated reactions were split into forward and backward parts, and the cost was added as a produced metabolite: initialreactionrxn:a+b⇔c+dreactionswithcostrxnf:a+b⇒c+d+mwrxn·expα·Δrg′rxnor·trxnb:c+d⇒a+b+mwrxn·exp−α·Δrg′rxnor·t a reaction which consumes this produced metabolite was then added to the model. while every original reaction in the model produces this cost at a rate relative to its absolute flux, this added reaction becomes the only one to consume this cost. minimizing the flux through this reaction will then return the flux distribution with the lowest enzyme associated cost to achieve the predefined biomass production value. reactions with no associated enzyme were given a cost of zero. these methods were applied to the e. coli ijr reconstruction. during all simulations the glucose uptake bound was set to an arbitrary value. uptake bounds for co fe h2o, h, na, nh o phosphate, and so were set to arbitrarily large values, much larger than the glucose uptake bound, and were considered to be present in excess. all exchange reactions release bounds were set to arbitrarily large values except for glucose, formate and o which were set to zero. predicted flux values were then normalized by either the flux of glucose to glucose-6-phosphate or overall flux through enzyme associated reactions. it is worth noting that this protein cost term does not account for the enzyme turnover rates or substrate affinity. we have opted for this cost term since the kinetic parameters of enzymes are extremely difficult to obtain experimentally, and vary greatly based on ph and temperature values. furthermore, the simplification of a constant turnover rate and substrate affinity for all enzymes has been shown to yield significant results. fundamental pathways analysis different methods have been proposed to characterize and decompose the fba solutions space, including elementary modes, extreme pathways and minimal generators. however, these tools have mostly been applied to small networks. the number of pathways obtained for each of these methods, as well as the computational cost of these analyses, increase exponentially with the size of the reconstruction. with large-scale metabolic reconstructions, these decompositions quickly become impractical tools. here we propose a modified version of extreme pathways which yields a significantly lower number of pathways during the decomposition. each of these pathways is associated with a protein cost and an atp production potential, which helps elucidate the model behavior as the value of the objective function varies. we find that this large number of pathways stems partially from the need to balance each pathway to a steady-state, such that no metabolite has a net production or consumption rate. certain metabolites, however, known as currency metabolites, are responsible for basic cell functions and occur in numerous reactions in the model. standard reactions involving these metabolites, such as atp hydrolysis, can be obtained through different combinations of other model reactions. these combinations, or loops, can then be used multiple times with the same purpose, increasing the number of pathways which have essentially the same core reactions. this phenomenon can be illustrated by extreme pathways analysis applied to the pathway and set of reactions presented in figure the pathway depicted in figure 1a is the production of energy through the conversion of glucose into lactate, which produces two mols of atp for every mol of glucose. reaction names are extracted from the e. coli reconstruction by orth et al.. figure 1b illustrates loops and reactions present in the same reconstruction that hydrolyze atp, any of which can be combined with the reactions in figure 1a to produce a pathway satisfying the atp steady-state condition. furthermore, atps4r also introduces a proton imbalance when hydrolyzing atp. this imbalance can then be corrected by any of the proton transport loops also presented in figure 1b. the combination of this pathway with the given loops then give us multiple extreme pathways with essentially the same core reactions; and they represent the same phenomenon - the conversion of glucose to lactate through glycolysis.figure metabolic loops increase pathway multiplicity. pathway converting mol of glucose to mols of lactate and producing mols of atp. reaction fluxes are represented in parentheses. extracellular metabolites are defined by. all other metabolites are in the cytosolic space. network loops and reactions that can hydrolyze atp and transport protons across the cell membrane. these reactions and loops can be used to balance the reactions in, yielding multiple balanced pathways with the same set of backbone reactions in order to extract only the core reactions in each pathway, we then remove a selected set of currency metabolites from the reconstruction. this removal reduces the number of pathways we obtain from extreme pathway analysis. this metabolite removal was performed under certain considerations: metabolites removed should be completely deleted from all reactions and compartments in the model. removing a metabolite from one reaction and not others could block certain reactions, since fewer options to balance that metabolite would be present in the reconstruction. no reaction should become a sink or source of metabolites. allowing reactions to become sinks or sources would introduce exchange reactions for certain metabolites, affecting the model decomposition. reactions removed from the reconstruction should not contain any metabolite which has not been removed from the reconstruction. the removal of these reactions would remove any pathway which uses these reactions to balance the metabolite still in the model. while we should be aware of these considerations when modifying the model, breaking any of these specifications should not invalidate the fundamental pathway analysis. the model decomposition could still be performed, but some limitations might have to be considered. we have illustrated this fundamental pathways method using the e. coli model by orth et al., and have demonstrated that the original pathways can be recovered by adding back the identified loops. this analysis can be found in the supplemental information. while we obtain fewer pathways in this tailored reconstruction, each set of reactions is associated with an imbalance in the original model. for example, instead of multiple extreme pathways, the set of reactions in figure 1a would yield a single fundamental pathway after the removal of currency metabolites from the reconstruction. these reactions alone, however, would yield atp, adp, phosphate, proton and water imbalances in the original model. in the fundamental pathways analysis, each resulting fundamental pathway is then associated with this imbalance, which can be used to further analyze each pathway. in this study, this imbalance is used to analyze the atp production capacity of each pathway obtained.  comparison between different cost functions cost-optimal flux simulations were performed using four different types of cost: molecular weights alone, thermodynamic penalties alone, the previously mentioned combination of molecular weights and thermodynamic penalties and uniform costs. the use of uniform costs gives a simple minimization of the overall flux through the enzyme associated reactions. this minimization strategy is a modified version of a widely regarded two-step optimization method called pfba. pfba has been previously used to predict unique flux distributions at the predicted fba optima. simulations comparing these four costs have been performed for growth rates ranging from to 100% of the predicted optima. flux distributions were normalized by the flux of glucose to glucose-6-phosphate, and results for several central carbon metabolism reactions can be found in the supplemental information. while simulations using molecular weights, thermodynamic costs, or a combination of both present considerably different flux distributions at different sub-optimal values, simulations using uniform costs quickly converge to a unique flux distribution as the objective value is lowered. the same set of simulations was also compared to several data points from two metabolic flux analysis experiments: ishii et al. and yao et al.. correlation and sum of squared error between all simulated and experimental flux distributions for several central carbon metabolism reactions were calculated. these results are presented in the supplemental information. a similar method to that used by holzhütter was also included in additional file 8: figure s while pfba simulations yield the highest correlation values in the suboptimal space when compared to the ishii et al. dataset, these same simulations also yield the highest sum of squared error. furthermore, when the objective function was lowered from 100% down to 65% of optima, corsofba yields the lowest sum of squared error for all data points except one, and the highest correlation values for the higher dilution rates in the yao et al. data. comparison between metabolic flux analysis and simulated fluxes in order to visualize the trend of how different fluxes change in the sub-optimal space, simulations utilizing molecular weight costs, thermodynamic costs and a combination of both, normalized again by the flux of glucose to glucose-6-phosphate, were plotted alongside the experimental values from the two mfa experiments. fluxes were plotted starting at 100% down to 75% of the predicted optima, and the results for several central carbon metabolism reactions are presented in figure the threshold of 75% was chosen in order to match the trend observed in our simulations to the experimental data.figure comparison between simulations and mfa experiments. comparison between selected simulated fluxes using molecular weights, thermodynamic penalties, and a combination of the two costs, and metabolic flux analysis data. fluxes are normalized by glucose to glucose-6-phosphate conversion rates. x-axes values are kept constant for all plots. reaction names are taken from the ecoli ijr model. yellow boxes indicate pentose phosphate pathway, blue indicates glycolysis, green indicates tca cycle and red indicates other reactions as the objective function value decreases, corsofba flux distributions show close agreement with mfa data for higher growth rates, particularly for simulations considering both thermodynamic and molecular weight costs. fluxes through the pentose phosphate pathway and glycolysis remain relatively constant, while the flux through the tca cycle gradually decreases. the decrease in tca cycle usage leads to higher levels of acetate release. quantitatively, the fluxes through glycolysis and the tca cycle are slightly under-predicted by our simulations, while the flux through the ppp is generally over-predicted. similar patterns have been observed in other fba method approaches, and this may arise due to the optimization of the biomass function alone. the optimization of atp production alongside biomass has been shown to yield better predictions. atp production optimization could also lead to higher fluxes through the tca cycle and lower ppp fluxes. it is also worth noting that, although fluxes through the tca cycle are under-predicted according to figure mfa experimental results can be inconsistent, and our predicted tca fluxes show good agreement with two other mfa studies. these simulations mirror a behavior known as overflow metabolism, where e. coli, at high growth rates, moves away from the full oxidation of glucose through the tca cycle and uses less atp efficient pathways, releasing acetate instead of co our simulations support the hypothesis that this behavior stems from a tradeoff between enzyme cost and energy yield. that is, when more glucose is available, e. coli uses a higher flux through pathways that are less enzymatically costly, but which produce fewer atp per mmol of glucose. predicted fluxes using molecular weights only and a combined cost also found excellent agreement with the experimental values in the reactions ppc and icl. both experiments show a decrease in icl flux and an increase in ppc flux, from near-optimal conditions to a growth rate of h- this result is further supported by the experimental results from nanchen et al., that found a lower flux through icl during a growth rate of approximately h- when compared to a growth rate of h- suggesting a transient response through this pathway as glucose concentration increases. the same transient response was found in our simulations. these results suggest that e. coli utilizes these anaplerotic reactions to relieve enzymatic cost, and that considering these costs during fba may increase internal flux predictions when comparing to the minimization of metabolic steps alone. simulation results also suggest that the optimization of the biomass function may yield fluxes that are fundamentally different from those at near-optimal conditions. at optimal growth, there is no predicted flux through the glyoxylate shunt, considerably lower fluxes through the tca cycle and ppp, and a flux through the ppc reaction not present at near-optimal conditions. although comparisons between experimental data and simulated flux distributions show that the highest correlation and lowest error are found at optimal growth, our simulations also indicate that the implementation of these costs yield comparable results in the sub-optimal space. moreover, the suboptimal corsofba approach can better predict fluxes through certain reactions, such as ppc and icl. while previous studies have shown that the fba solution space increases drastically when considering the objective function at near-optimal to optimal values, here we show that the optimization of the enzymatic cost at near-optimal conditions yields results that are more consistent with experimental data for certain reactions. this result strongly suggests that exploring the fba solution space at near-optimal to optimal objective values may increase the predictive accuracy of internal cell fluxes. the near optimal solution space was also analyzed for knockout strains. flux distributions for six e. coli knockout strains reported by ishii et al. were compared to cost-optimal simulations using the combined cost function. results show that cost-optimal flux distribution from 95% to optimal yield similar or better predictions according to both correlation and sum of square error. predictions using uniform costs quickly diverge from experimental data. these simulations further support the “cloud theory” proposed by wintermute et al., and suggest that the experimental data can also be in good agreement with near-optimal flux distributions. comparison between gene expression data and simulated fluxes to further validate this analysis, simulations using the combined cost function were also compared to gene expression values reported in the mfa studies. while the overall correlation between reaction flux and gene expression is moderate at best, increased expression of all genes participating in a particular pathway can be considered a good indication of increased flux. simulations were performed under the same conditions as before, utilizing the combination of molecular weights and thermodynamic penalty, but this time the flux distribution was normalized by the overall flux through the enzyme associated reactions. comparisons between simulated values and gene expression data are shown in figure figure comparison between simulations and gene expression data. comparison between selected simulated fluxes, using a combination of molecular weights and thermodynamic penalties, and associated gene expression. plot axes are the same as the ones defined in the ptsg plot unless otherwise specified. multiple genes associated with the same reaction are included in the same box. reaction names are taken from the e. coli ijr reconstruction. yellow boxes indicate pentose phosphate pathway, blue indicates glycolysis, green indicates tca cycle and red indicates other reactions although gene expression data can also be inconsistent between experiments, and some reactions are associated with multiple genes, this comparison further supports the qualitative results presented in the previous section. when normalizing the predicted fluxes by the overall flux through the enzyme associated reactions, an increased relative uptake of glucose is observed as the objective value decreases. as a result, an increased relative flux through glycolysis is also predicted. furthermore, the relative flux through the tca cycle still decreases as the growth rate increases. the increase in relative glucose uptake and flux through glycolysis, as well as the decreased flux through the tca cycle, are supported by the relative gene expression associated with these pathways. simulated values also find good agreement with fluxes through the entner-doudoroff pathway. genes edd and eda, associated with this pathway, have long been known to exist in e. coli, but the activity of their associated reactions has been observed mainly under growth on gluconate, glucuronate, and methyl-beta-d-glucuronide; phosphate limitation; and carbon starvation. due to this fact, these reactions are generally not included in mfa experiment networks. in contrast with murray et al., our simulations using molecular weight and combined costs predict the use of the ed pathway under excess glucose conditions, not starvation. these predictions are supported both by the genetic data presented by the mfa studies, which show an increase in edd and eda expression at high glucose concentrations, and mfa experimental results by harcombe et al., where fluxes through the ed pathway were measured in bacteria growing at growth rates above h- fundamental pathways analysis to better understand the transition between metabolic states as the value of the objective function is decreased, the energy producing pathways of the e. coli ijr reconstruction were decomposed using the fundamental pathways analysis described in the methods section. details on how these reactions were calculated can be found in the supplemental information. briefly, starting with the atp imbalance associated with each fundamental pathway, the total atp potential of each pathway was estimated based on the imbalance of other energy generating metabolites, such as nadh and ubiquinol- the lowest cost of converting these metabolites into atp was then added to the total cost associated with the pathway. the final atp production by enzymatic cost was then compared to the potential atp production by mmol of glucose.figure fundamental pathway analysis. atp production potential of energy producing fundamental pathways, compared to the associated protein cost of each mmol of atp. optimal and near optimal pathways are highlighted in the scatter plot. ops are plotted in blue and nops in red in the central carbon metabolism diagrams in accordance with general e. coli metabolism knowledge, the fundamental pathway analysis found two optimal pathways, according to atp production by both mmol of glucose and protein cost. the most energy productive pathway is the full oxidation of glucose through the tca cycle. the most cost efficient pathway was found to be the production of acetate through glycolysis. although this pathway potentially produces less than half the amount of atp as op the total enzymatic cost per mmol of atp is lower. also in agreement with general knowledge, this analysis predicted the release of acetate to be more efficient than the release of both lactate and ethanol. in the simulations presented in previous sections, a gradual transition from op to op is observed as the objective function value decreases. this trend also supports the idea that overflow metabolism takes place in order to alleviate enzymatic cost. under low glucose conditions, e. coli would need to extract from glucose the largest possible amount of energy in order to sustain growth. as more glucose becomes available, however, this bacteria can afford to consume more substrate and produce energy through cheaper pathways. one interesting observation from this analysis is the existence of near optimal pathways, which combine the previously described optimal pathways with the ppp or ed pathways. nop and nop both produce more atp per mmol of glucose than op but less than op while remaining less cost efficient than op but slightly more so than op nop and nop reduce the flux through upper glycolysis by partially re-routing the flux through the ppp, and they exhibit no flux through phosphoglucose isomerase. a third near optimal pathway, termed nop skips the upper glycolysis through the ed pathway. this pathway is both less cost and energy effective than op the use of these near optimal pathways could explain the inconsistent ppp fluxes reported by the mfa experiments considered here. a short recovery in the simulated ppp flux near a growth rate of h − demonstrates that this pathway is a viable option at certain dilution rates, due to its coupling with the ppp and biomass production. it has also been demonstrated that Δpfka deficient e. coli strains reduce flux through phosphofructokinase by diverting fluxes through the ppp, hence using nop and nop instead of op or op when the cost through upper glycolysis is increased. this fundamental pathway analysis elucidates how we would expect the model to behave in terms of energy generation as we optimize the protein cost at sub-optimal states. based on this analysis, in order to fulfill its catabolic needs, the model transitions from op to op as we move away from optimal growth. although this is in fact the general trend we observe, anaplerotic needs also need to be considered. these needs are addressed most efficiently by reactions not in op or op such as pep carboxylase and the glyoxylate shunt. interestingly, the model also predicts the use of the ed pathway with anaplerotic purposes at high growth rates. while this pathway has been studied mostly for its catabolic activity in z. mobilis and several pseudomonas species, simulations here predict this pathway to be the cheapest way to get glucose shuttled to the tca cycle for the production of building blocks, while giving up little efficiency in energy pathways through nop  here we present corsofba, a method for analyzing and predicting metabolic flux distributions at sub-optimal states based on protein and thermodynamic cost minimization. corsofba demonstrates that the optimization of protein cost at near-optimal states can produce significantly different results from those at the optimal objective. furthermore, although correlation and error calculations indicate better predictions at the optimal solution, protein-optimal results at sub-optimal objectives show better agreement with mfa experiments and gene expression profiles for anaplerotic reactions. these results suggest it is important to explore the sub-optimal fba space in order to better predict internal cell fluxes. although the method described here is not suited for predicting growth rates, it provides a platform for analyzing internal cell fluxes in the sub-optimal space. we believe corsofba will be particularly useful given recent studies have suggested e. coli can exist freely within a near-optimal space. we also believe this method can be useful in predicting fluxes in healthy, multicellular organisms, which have more complex objectives than the production of biomass. future studies are merited to explore whether the predictive power of methods which optimize the objective function under enzymatic level constraints, such as moment and fbawmc, would also benefit from exploring the sub-optimal solution space. our results suggest that significant conclusions could be drawn by adapting these methods to optimize the enzymatic cost they implement in a near-optimal space, rather than using these simply as a model constraint. additional files additional file figure s sensitivity analysis of parameter α. model simulations ranging from to 75% of optima with varying values of α. simulations are performed using both thermodynamic costs only and combined costs. additional file gibbs free energy change distribution. file showing plot and distribution of the gibbs free energy change of reactions considered. additional file molecular weight and cost calculation. details on how molecular weight values were obtained and final costs calculated. molecular weights and final costs distribution are also plotted. additional file table s thermodynamic, molecular weight and final cost values. table includes all values used during simulations: gibbs energy of formation for all metabolites, calculated gibbs energy change for all reactions based on metabolite energy of formation and molecular weights obtained from databases. additional file fundamental pathway analysis. example demonstrating the recovery of extreme pathways by adding reactions loops back into the fundamental pathways. details on how the fundamental pathway analysis presented in the manuscript was performed are also available. additional file table s model decomposition. extreme pathways and fundamental pathways calculated both for the analysis presented in the manuscript and the pathway recovery analysis in additional file additional file figure s cost values comparison. cost-optimal simulations using molecular weights alone, thermodynamic penalties alone, a combination of the two and uniform costs are presented for major central carbon metabolism reactions in the range of to 50% of maximal growth. additional file figure s error and correlation calculations. correlation and sum of squared error are calculated between all simulated and experimental flux distributions for major central carbon metabolism reactions. additional file figure s error and correlation calculations for knockout strains. correlation and sum of squared error are calculated between mfa experimental data and simulated flux distributions using the combined cost function for six knockout strains. abbreviations tcatricarboxylic acid fbaflux balance analysis ppppentose phosphate pathway mfametabolic flux analysis edentner-doudoroff glcglucose g6pglucose-6-phosphate f6pfructose-6-phosphate fdpfructose- -biphosphate dhapdihydroxyacetone phosphate g3pglyceraldehyde 3-phosphate 13dpg3-phospho-d-glyceroyl phosphate 3pg3-phospho-d-glycerate 2pgd-glycerate 2-phosphate pepphosphoenolpyruvate pyrpyruvate accoaacetyl-coa actpacetyl phosphate acacetate citcitrate icitisocitrate akg2-oxoglutarate succoasuccinyl-coa succsuccinate fumfumarate mal-ll-malate oaaoxaloacetate 6pgl6-phospho-d-glucono- -lactone 6pgc6-phospho-d-gluconate ru5p-dd-ribulose 5-phosphate r5palpha-d-ribose 5-phosphate xu5p-dd-xylulose 5-phosphate s7psedoheptulose 7-phosphate e4pd-erythrose 4-phosphate 2ddg6p2-dehydro-3-deoxy-d-gluconate 6-phosphate competing interests the authors declare that they have no competing interests. authors’ contributions as and aq conceived the methods and revised the manuscript. as drafted the manuscript and generated the results. both authors have read the manuscript and approved the final version. 
maintaining the objective function close to its maximum value predicts metabolic states that closely resemble low dilution rates; while higher dilution rates can be mirrored by lowering the biomass production value. we present an innovative fba method called corsofba based on the optimization of protein cost at sub-optimal objective levels. metabolic flux distributions at the optimal objective can be substantially different from the near-optimal distributions. this objective is generally defined for unicellular organisms by a theoretical reaction which simulates biomass production. studies have suggested, however, that unicellular organisms - like multicellular organisms - do not grow at optimal rates. fba has been extremely successful at predicting in e. coli growth rates under different media and gene essentiality, amongst other things. our method shows good agreement with experimental data of e. coli grown at different dilution rates. importantly, the behavior of e. coli central carbon metabolism can be better predicted by exploring the sub-optimal fba solution space. flux balance analysis is a widely used tool to model metabolic behavior and cellular function. this technique can be used to better predict internal cell fluxes under different conditions, and corsofba will be of great help for the study of cells from multicellular organisms using flux balance analysis. flux balance analysiscobrabigggenome-wide metabolic reconstructionsconstraint based modelsextreme pathwayssub-optimal growthprotein costissue-copyright-statement© the author 2015 to further improve fba predictions, particularly of internal cell fluxes, new techniques to explore the sub-optimal solution space need to be developed. the corsofba method presented here is able to predict the behavior of pep carboxylase, the glyoxylate shunt and the entner-doudoroff pathway at different glucose levels, a behavior not predicted by the minimization of metabolic steps and fba alone. by using a modified version of extreme pathways, we are also able to quantify the energy production and overall protein cost for all possible pathways in the central carbon metabolism. in order to improve predictions, additional constraints are coupled with optimization of the biomass function.
 flux balance analysis is a widely used tool to model metabolic behavior and cellular function. applications of fba span a breadth of research from synthetic engineering of biofuels to understanding evolutionary adaptations. fba predicts metabolic reaction fluxes that optimize a given objective. this objective is generally defined for unicellular organisms by a theoretical reaction which simulates biomass production. fba has been extremely successful at predicting in e. coli growth rates under different media and gene essentiality, amongst other things. in order to improve predictions, additional constraints are coupled with optimization of the biomass function. studies have suggested, however, that unicellular organisms - like multicellular organisms - do not grow at optimal rates. to further improve fba predictions, particularly of internal cell fluxes, new techniques to explore the sub-optimal solution space need to be developed.  we present an innovative fba method called corsofba based on the optimization of protein cost at sub-optimal objective levels. our method shows good agreement with experimental data of e. coli grown at different dilution rates. maintaining the objective function close to its maximum value predicts metabolic states that closely resemble low dilution rates; while higher dilution rates can be mirrored by lowering the biomass production value. by using a modified version of extreme pathways, we are also able to quantify the energy production and overall protein cost for all possible pathways in the central carbon metabolism.  metabolic flux distributions at the optimal objective can be substantially different from the near-optimal distributions. importantly, the behavior of e. coli central carbon metabolism can be better predicted by exploring the sub-optimal fba solution space. the corsofba method presented here is able to predict the behavior of pep carboxylase, the glyoxylate shunt and the entner-doudoroff pathway at different glucose levels, a behavior not predicted by the minimization of metabolic steps and fba alone. this technique can be used to better predict internal cell fluxes under different conditions, and corsofba will be of great help for the study of cells from multicellular organisms using flux balance analysis. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords flux balance analysiscobrabigggenome-wide metabolic reconstructionsconstraint based modelsextreme pathwayssub-optimal growthprotein costissue-copyright-statement© the author 2015 