SCIENTIFIC REPpRTS

OPEN Robust Analysis of Fluxes in

Genome-Scale Metabolic Pathways

Received: 23 August 2016 Accepted: 13 February 2017 Published online: 21 March 2017

Michael MacGillivray1, Amy Ko1, Emily Gruber1, Miranda Sawyer1, Eivind Almaas 2,3 & Allen Holder1

Constraint-based optimization, such as flux balance analysis (FBA), has become a standard systems-biology computational method to study cellular metabolisms that are assumed to be in a steady state of optimal growth. The methods are based on optimization while assuming (i) equilibrium of a linear system of ordinary differential equations, and (ii) deterministic data. However, the steady-state assumption is biologically imperfect, and several key stoichiometric coefficients are experimentally inferred from situations of inherent variation. We propose an approach that explicitly acknowledges heterogeneity and conducts a robust analysis of metabolic pathways (RAMP). The basic assumption of steady state is relaxed, and we model the innate heterogeneity of cells probabilistically. Our mathematical study of the stochastic problem shows that FBA is a limiting case of our RAMP method. Moreover, RAMP has the properties that: A) metabolic states are (Lipschitz) continuous with regards to the probabilistic modeling parameters, B) convergent metabolic states are solutions to the deterministic FBA paradigm as the stochastic elements dissipate, and C) RAMP can identify biologically tolerable diversity of a metabolic network in an optimized culture. We benchmark RAMP against traditional FBA on genome-scale metabolic reconstructed models of E. coli, calculating essential genes and comparing with experimental flux data.

Constraint-based analysis is an approach to study metabolic networks that has become a primary computational tool with a well established literature1-3. The foundation of constraint-based analysis is a representation of a cell's metabolism as a linear system of differential equations in the (unknown) fluxes of the metabolic reactions. The additional assumption that the metabolic network has reached steady state results in a homogeneous system. It is customary to optimize an appropriate objective over the fluxes that achieve steady state to model how a cellular metabolism reacts as it is held in a constant environment4, 5. Different objectives have been proposed and studied6, 7, with the most prevalent being the maximization of the rate at which biomass is created.

The most commonly applied version of constraint-based analysis is known as flux balance analysis (FBA). The original study of FBA investigated small portions of the central metabolism8, but the use of constraint-based methods has expanded to encompass genome-scale metabolic models and a diverse range of applications and computational formulations9. However, the basic premises of optimality and steady state have remained intact, although these assumptions are well known to be inexact from a biochemistry perspective. That said, there are two somewhat evident reasons to leave these premises unchecked. First, the resulting optimization problems are most often linear and less often quadratic (convex), and in both cases the problems are easily solved with standard methods and software. Second, traditional constraint-based models are useful even with questionable premises. For example, with a high level of accuracy, FBA models can identify essential genes10, predict metabolic responses as pathways are interrupted11, determine central metabolic pathways12-14, develop synthetic biology engineering strategies15-18, as well as identify genes showing promise as inhibitors of cancer migration19. Hence, a justified retort by the constraint-based analysis community to those who might question the two defining premises is that the modeling approach yields meaningful and computationally tractable science even though the cornerstone premises are imperfect.

The traditional constraint-based paradigm assumes a model of an individual cell that has reached an ideal, limiting setting of steady state. Models are benchmarked with experimental outcomes under the assumption that experimental measurements are conducted on a population of identical cells that have evolved to optimal

department of Mathematics, Rose-Hulman Institute of Technology, Terre Haute, IN, USA. 2Department of Biotechnology, NTNU - Norwegian University of Science and Technology, Trondheim, Norway. 3K.G. Jebsen Center for Genetic Epidemiology, Department of Public Health and General Practice, NTNU - Norwegian University of Science and Technology, Trondheim, Norway. Correspondence and requests for materials should be addressed to E.A. (email: eivind.almaas@ntnu.no) or A.H. (email: holder@rose-hulman.edu)

performance. However, experimental reality differs from this assumption. First, even well optimized cellular populations of identical cells exhibit prominent heterogeneity in uptake, secretion, and growth rates4. Second, experimental flux measurements constitute averages over such heterogeneous populations. Moreover, these averages provide the only possible experimental prototype of cellular flux since single cell measurements are currently impossible. Third, heterogeneity of metabolic rates is innate within a culture because cells vary within, e.g., their cell cycles and replication states. The traditional paradigm that a constraint-based model represents a single, optimized cell is certainly challenged by such experimental realities upon which models are vetted.

Here, we challenge and relax the foundational assumption of a system being in a steady state, and we instead introduce a robust optimization20 counterpart to FBA, called Robust Analysis of Metabolic Pathways (RAMP). The robust optimization implementation allows us to model a system stochastically where, instead of imposing the rigid and simplistic condition of deterministic coefficients and steady state, we allow these assumptions to be relaxed. In RAMP we control departures from steady state by limiting their likelihood of deviation. Since the stochastic flexibility of RAMP permits us to postulate metabolic models for cultures that deviate from steady state, this opens the possibility of computationally studying functional states of cellular metabolisms as they transition toward a steady state. Such an investigation begs the questions of (1) whether or not RAMP solutions are continuous in their stochastic elements and (2) whether or not convergent trajectories of RAMP solutions yield deterministic FBA solutions as the probabilistic assumptions abate. We mathematically answer both questions in the affirmative, and hence, FBA is justly thought of as the limit of evolutionary adaptations of a cellular metabolic network that achieves steady state.

Additionally, the development and stochastic formulation of RAMP explicitly allows us to systematically address heterogeneity in metabolic phenotypes that exists in isogenic cellular populations21-23. Our mathematical analysis identifies the probabilistic situations that coincide with our modeling framework once the limiting deterministic setting is reached. Consequently, if we assume that RAMP's modeling paradigm aligns with the stochastic nature of a culture, then we can mathematically characterize the ways in which a culture may harbor variation.

Possible concerns with RAMP are that it could lose the predictive successes of FBA, and that it could be less computationally tractable. We have designed our computational experiments to compare RAMP with FBA in such a way that both models can be solved with the default simplex algorithms of several linear solvers, and our results on two E. coli models24, 25 demonstrate that RAMP rivals FBA in its ability to identify essential genes. We further show that RAMP's efficacy in identifying essential genes is predominantly stable as individual coefficients of the biomass equation are assumed to be uncertain. The ranges over which individual coefficients can vary are disparate, having a surprising spread of many orders of magnitude. All coefficients could accommodate uncertainty of at least 0.42%, with most extending beyond 100% uncertainty. We also compare FBA's and RAMP's consistency with experimentally determined fluxes6, 26, 27, finding that RAMP significantly outperforms FBA for both aerobic and anaerobic conditions.

With regards to computational tractability, RAMP is a robust linear program, which is a second-order cone program (SOCP) that is known to be solvable in polynomial time28, 29. Such robust models were developed to overcome problems with over-optimizing designs, meaning that a design could be optimal with respect to the estimated data on which it was built but (much) less so as the data varied over realistic possibilities30. A similar concern about over-optimization has been expressed in the FBA literature7, and RAMP directly addresses this issue by expressing variability within the model itself.

We note that an alternative robust model has been suggested31. This method differs significantly from RAMP since it uses a robust least squares method in a bi-objective model. Such a model can illustrate the trade-off between the optimal growth rate and the deviation from steady state, whereas RAMP uses robust linear programming. Furthermore, the numerical work in ref. 31 is conducted on a small, illustrative example, and there is no extrapolation of the bi-objective model to a realistic genome-scale metabolic reconstruction.

Our conclusion is that RAMP opens a new paradigm for constraint-based optimization and analysis of genome-scale metabolic models that is also directly applicable to a large number of the recent constraint-based extensions to FBA. Our computational implementation and benchmarking of a specific version of RAMP using both linear and non-nonlinear objective functions, show that it is a computationally tractable alternative and complement to traditional FBA.

Results

We begin by re-visiting the foundational assumptions of FBA, where the metabolism of an individual cell c is defined by a system of n reactions, with the j-th reaction being

ajxt] + aCjixC] + ... + acmj[xcm] ^ b¿[xfl + bc2j[xc2] + ... + bcmj[xcm].

Here, [x,c] denotes the concentration of the i-th metabolite in cell c. The cellular superscript c is absent in FBA's standard development due to the assumption that an average cell is being studied. Thus, the ordinary differential equation for the i-th metabolic concentration in cell c is,

d[x,c] E b - afj)(k+ [xtfv [x^i... [xm ]am - k-tfpi [x^i... [xm ]bmj)

dt i=1

i=1 (1)

Input Metabolites

ACCOA(-0.00005) ALA(-0.488) AMP(-0.001) ARG(-0.281)

ASN(-0.229) ASP(-0.229) ATP(-45.73) CL(-0.00645)

COA(-0.000006) CTP(-0.126) CYS(-0.087) DATP(-0.0247)

DCTP(-0.0254) DGTP(-0.0254) DTTP(-0.0247) FAD(-0.00001)

GLN(-0.25) GLU(-0.25) GLY(-0.582) GLYCOGEN(-0.154)

GTP(-0.203) HIS(-0.09) ILE(-0.276) LEU(-0.428)

LPS(-0.0084) LYS(-0.326) MET(-0.146) MTHF(-0.05)

NAD(-0.00215) NADH(-0.00005) NADP(-0.00013) NADPH(-0.0004)

PE(-0.09675) PEPTIDO(-0.0276) PG(-0.02322) PHE(-0.176)

PRO(-0.21) PS(-0.00258) PTRC(-0.035) SER(-0.205)

SPMD(-0.007) SUCCOA(-0.000003) THR(-0.241) TRP(-0.054)

TYR(-0.131) UDPG(-0.003) UTP(-0.136) VAL(-0.402)

ADP(45.560000) Biomass(1.000000) PI(45.560000) PPI(0.730200)

Output Metabolites

Table 1. An example of the non-integer coefficients, in parentheses, of the input and output metabolites of the E. coli metabolic model iJR904 growth reaction36.

where Sc is the stoichiometric coefficient for metabolite i in reaction j, the flux of reaction j is vc, and c denotes the cell. Consequently, genetically identical individuals of a culture share the same collection of metabolic reactions. Hence, the preponderance of stoichiometric coefficients of any particular genome-scale metabolic reconstruction are static across the individuals of a culture. However, a complete metabolic model includes several inferred and less certain reactions to model processes, such as those involved in the creation of biomass. In many of the recent extensions of FBA9, a combination of alternative objective functions, global constraints or upper/lower bounds on sets of reaction fluxes are added to the system. Consequently, many extensions, such as FBAwMC32, GIMME33, MOMENT34, and GX-FBA35, are implemented by adding rows and/or columns to the FBA stoichiometric matrix, and these new stoichiometric coefficients do not generally consist of integer values. Instead, depending on the method of interest, these added stoichiometric coefficients may consist of a combination of data from large-scale' omics sources, biochemistry knowledge, or theoretical modelling assumptions. Thus, they are carriers of inherent uncertainty.

In particular, the experimentally determined coefficients of a biomass equation are noticeably different from their integral counterparts in the shared metabolism, see Table 1 as an example from the E. coli metabolic model iJR90436. These less certain coefficients would be stochastic, time varying parameters per individual cell. Thus, instead of the growth reaction in Table 1 being imposed on all cells in a population, these coefficients should be able to vary among the individual cells.

There exists also other sources of uncertainty for the performance of a genome-scale metabolic reconstruction, such as the value of the phosphate/oxygen (P/O) ratio. Changes in how the network performs the ATP production and drain has a significant potential to affect model-predicted growth performance. In the following discussion, however, we will base our presentation and discussion on uncertainty in the biomass growth coefficients.

Suppose c indexes the single, average cell assumed by FBA. The limiting steady state assumption of FBA constrains the fluxes so that,

= ±4j = 0, vi, dt £j j

which is the equilibrium condition as t^ <. Furthermore, it has been suggested that cellular metabolic performance evolves toward an optimal growth state as a culture's environment is left unchanged and as populations are repeatedly sampled and re-grown4' 10. FBA is predicated on the assumption that cells have been tuned through such an evolutionary processes to an "optimal" use of their resources, where "optimal" is measured by a function of the fluxes, g(v). Hence, FBA studies metabolic processes through optimization problems, and adaptations thereof, of the form

max{g(vc) : Scvc = 0, Lc < vc < Uc}, (2)

where Sc is the stoichiometric matrix whose components are Sc and vc is the associated (decision) vector of fluxes. The objective function g(v) is typically chosen to estimate the rate at which biomass is created, and for convenience, we assume throughout that g(v) is the sole flux of the growth reaction, i.e. g(v) = vGrowth. While several alternative objectives have been proposed and investigated in different biological settings6, 7, the production of biomass is the most prevalent choice and is the default in the COBRA37 toolbox. The vectors of lower bounds, L, and upper bounds, Uc, may contain ±<, or some suitably large value, to indicate that a flux is unbounded. If Lc < 0 and U > 0, then reaction j is reversible.

A probabilistic development of RAMP. We assume that a culture has gone through a long-enough time of growth and re-sampling to be largely optimized to the resources of its (invariant) environment. Heterogeneity within the culture is modeled as the random experiment that is the growth of the culture, and we consider the variational states of the genetically identical cells of the culture as independent and identically distributed outcomes from that growth experiment. Hence, the known heterogeneity within an optimized culture is the result of a large sample of possible outcomes. Consequently, for each metabolite i the collection of d[xf]/dt across the cells c of the culture is a random sample of the stochastic rate of change of the concentration. Let C be the collection of cells in the culture and d[xf]/dt be the average rate of change of the concentration (of metabolite i) over the sample, that is

d[x,c] _ v d[x,c]

dt _ |C| ctC dt '

The central limit theorem mandates that d[x,c]/dt be approximately normal, and we let ¡, and a, be the unknown mean and standard deviation of this normal distribution. This approximation is expected to be trustworthy since cultures regularly contain hundreds of millions of cells.

We replace the traditional static assumption of steady state with the probability constraints

P(d[x,c]/dt > M) < s and P(d[x,c]/dt < -M) < s, (3)

which combine to ensure that P(—M, < d[x,c]/dt < M,) > 1 — 2e. Here, |M,| is a bound on the feasible fluxes that define the allowed deviation from steady state. Unlike the FBA model formulation that collapses onto an ideal cell and ignores cultural variations, these stochastic constraints account for the innate variations that are known to exist. Since (d[x,c]/dt — ¡1)/ai is a standard normal variable, we write the constraints in Eq. (3) as

Mt - f, -Mj - f, -- > S1-£ and-- < SE,

where Se and Sl—e are the e and 1 — e percentiles respectively, i.e. P((d[x,c]/dt — ¡j,)/Oj > ¿1—e) = e. Rearranging these inequalities and using the fact that Sl—e = _Se, we find

[ Mi - f , M, + f, 1

< < min-!-'-,-- !-.

i 6ls 6ls(4)

The stochastic constraints in Eq. (4) are expressed in terms of the fluxes by extending the static differential equation of an individual cell in Eq. (1) to the stochastic differential equation of a culture. The standard assumption has been that the ideal cell of an FBA model is an average representation that matches measurements of a culture. If we allow S, to be the i-th row of the random stoichiometric matrix in which the uncertain elements are random variables, e.g. for the creation of biomass, then the implicit FBA assumption of an ideal cell results in the following stochastic differential equation for each metabolite,

d[xc] d[xf]

-=-= Sv. .r.

dt dt i (5)

Consequently, for any flux vector v we have

= E(d[xf]/dt) = E(d[x f]/dt) = E(Siv) = E(S,)v.

An apt interpretation of FBA's projection of a culture's heterogeneity onto an average is somewhat laid bare by these equalities: FBA assumes that ¡i i = 0. The first and second equalities then constrain the expectation of the sample mean, d[xf] /dt to be zero. The ideal FBA cell c is interpreted as a stochastic representation of the culture; a representation that inherits the expected averages of the metabolic rates of change over a culture. The FBA paradigm thus acknowledges, but ignores within its modeling framework, deviation from steady state since FBA only insists that the expectation be zero without regard to variation. The third equality is a direct consequence of Eq. (5), and the last equality follows from the linearity of the expectation. Hence FBA assumes E(S )v = 0, which shows that the stoichiometric coefficients of an FBA model are reasonably interpreted as averages, i.e. that Sf = E(S'). Hence, growth coefficients, such as those in Table 1, are expectations attributed to an ideal cell in hopes of representing an entire culture. If a stoichiometric coefficient is certain because it is part of the metabolic network that is common to the taxa, then the average is simply the known value that is shared among all cells. If the coefficient is otherwise uncertain and assumed to be random, then FBA assumes the mean value.

The equation E(S) v = i, shows how the vector of fluxes v defines the expected value of the rate of change of the concentration (of metabolite ), and the equation begs for detailed random models of the uncertain stoichiometric coefficients. Such detail is challenged by experimental limitations. For example, it is impossible to study individual cells or even groups of cells that are in identical metabolic or proteomic states. Hence, inferring accurate distributional characteristics on parameters such as the growth coefficients is currently impossible. Beyond averages, our ability to express Eq. (4) in terms of the fluxes further necessitates the calculation of standard deviations.

Scenario analysis in RAMP. A common statistical supposition is to impose additional requirements. For example, we could outright assume that the individual random elements of S are independent and identically

distributed. However, such a requirement would be dubious in our setting due to the certain dependence among growth and metabolic uptake and excretion. We instead promote the use of scenario analysis to query the behavior of the stochastic model. Scenarios have been used in other application domains in which similar optimization problems arise. For instance, scenarios are motivated even with an assumption of independent and identically distributed random variables to better manage computational space38. In addition to overcoming our inability to accurately model the separate random stoichiometric coefficients, scenarios have the practical advantage of being easily interpreted as biological possibilities. Importantly, the normality of d [xf] / dt is independent of the distributions used to model the uncertain stoichiometric coefficients. "tte scenarios in RAMP's development and mathematical analysis are purposely arbitrary so that RAMP remains viable no matter which (discrete) distributions might actually model a biochemical reality.

Assuming that each St has q random scenarios, we let Sik be the (non-random) stoichiometric coefficients for the i-th metabolite in scenario k and let pik be the probability of scenario k, for k = 1, 2, 3, ..., q. Allowing pi to be the (positive) vector of probabilities indexed by k, Pi to be the (positive definite) diagonal matrix formed by p i, and Sj to be the matrix whose k-th row is Sik, the mean and variance of Siv are

J = E(Siv) = E(Si)v = itpikSkv = P?S,v

k=1 (6)

of = Var(Sv) = j2pik (Sikv - E(Sikv))2 = (§,v)T(I - epj)TP(I - epj)§,v,

where e is a vector of ones. Defining R{ = S1—£^[Pj (I — epT)Sj, we have that the standard deviation is

q = vvarcSv = V pv №e = ii-mi/¿w (7)

The expected value and standard deviation are related through constraint (4), from which we have

or equivalently,

Pvll < min{Mi - pTSiv, Ml + pTSiv},

PVI - Mt < pTSv < Mt - HRvll. (8)

Although Eq. (8) is more complicated than a traditional linear constraint, it is the combination of two second order cone constraints. Such constraints share several desirable properties with linear constraints, such as being convex. Replacing the traditional linear constraints Scv = 0 with Eq. (8), we have arrived at an SOCP that is our RAMP model:

max vGrowth

subject to

PVI - Mi < pjSv < Mi - \\RjV\\, i = 1, 2, L < v < U.

"tte RAMP model defaults to a commonality of the upper and lower flux bounds among the cells so that Lc = L and Uc = U for all c. "ttis default is somewhat for notational convenience since stochastic variation in L or U could be remodeled as variation in the coefficient of the associated flux, and the resulting constraint would then be part of the stoichiometric matrix20.

A feasible RAMP solution is a collection of fluxes under which the mean and variance of the normal distribution of d[xic]/dt satisfies the likelihood of being in steady state with regard to the spectrum of stoichiometric possibilities. A mathematical expression of this fact is possible upon recognizing that

PVI = max {uTRv : ||u|| < 1}, which allows the right-hand side of the SOCP constraint to be re-written as

pjsv + Pvll < Mt o pT§iv + max{uTRiv: ||u|| < 1} < Mt O max^TSv + uTRv: ||u|| < 1} < Mi

O Sv < M, V S{ e {pJSi + uTRi: ||u|| < 1}. (10)

"tte left-hand side SOCP constraint can be similarly re-written as,

-pjsvv + pivy < Ml O SiV < Ml V Si e {-pJSi + uTRi: ||u|| < 1}. (11)

The sets from which Si are drawn are called uncertainty sets in the robust optimization literature, and these sets define the collection of stoichiometric possibility. A feasible flux must satisfy being within |Mi| of steady state

for all possible stoichiometric possibilities, meaning that S i can deviate from its average pTSj by as much as uTRi as long as |\u\\ < 1.

Illustrative example ofthe RAMP method. The re-expression of the SOCP constraints in terms of their uncertainty sets provides a geometric description. As an illustrative example, consider the two variable system

where the average values are s1 = 1 and s2 = — 1, corresponding to the standard static values of a stoichiometric matrix. Assume we have q = 3 scenarios, and

which means that

1/4' 0.9 -1.2

p1 = 1/2 and S1 = 1 -1

1/4 1.1 -0.8

P(sj = 0.9 and s2 = -1.2) = P(sj = 1.1 and s2 = -0.8) = 1/4

P(sj = 1 and s2 = — 1) The expected value of the constraint is simply

v = v - v2

3 so that e 0.0015, we have that

pTSv :

which is the value of the original, static constraint. Setting 61—e

-0.15 -0.30 0 0 0.15 0.30

If M1 = 0.2, then Eqs (10) and (11) combine to show that vl and v2 are feasible in the stochastic model if and only if

-0.2 < (1 - 0.15(uj - u3))v + (- 1 - 0.30(uj - u3))v2 < 0.2

for all vectors u such that uj2 + u2 + u3 < 1. Hence, v1 and v2 are feasible in the stochastic case only if they satisfy an infinite number of linear inequalities that are perturbations of a relaxed version of the original, static equality.

This infinite collection of constraints is neither more nor less restrictive than the FBA constraint. One interpretation is that RAMP first relaxes the average equilibrium constraint of pTSjV = 0 by replacing this last expression with —Mj < pTSjV < but then RAMP further restricts the relaxed constraint by replacing it with ||R,v|| — Mj < pTSjV < Mj — ||Riv||, which adds an infinite number of linear constraints (per original constraint).

A depiction of the geometry for the above example is shown in Fig. 1. The original, static (FBA) constraint with the variable bounds is depicted by the bold dashed line segment through the origin. The shaded region is the feasibility set formed by the combined stochastic constraints, and the light dashed lines are samples of the infinite set of linear inequalities added by the stochastic model. The RAMP feasible region does not contain the FBA feasible region. As an example, the point (0.75, 0.75) is feasible to the original static FBA constraint but is infeasible to the stochastic RAMP constraints. Likewise, there are feasible solutions to the stochastic model that are infeasible in the static case.

Mathematical observations of RAMP properties. In the following we develop three important mathematical qualities of the RAMP model by establishing three Theorems, proofs of which are contained in the Methods section. Our first result demonstrates that feasible RAMP flux vectors satisfy a Lipschitz continuity property in their probabilistic elements, provided that bounding constraints can be marginally relaxed. The second result demonstrates that convergent RAMP solutions are indeed FBA solutions as stochastic uncertainty diminishes. The third theorem characterizes the types of uncertainty in the stoichiometric matrix that can be tolerated by an optimal FBA solution.

Continuityproperty ofRAMP solutions. Corollary 1 (see Methods) suggests that RAMP's SOCP constraints might imbue continuity of the feasible flux vectors with regard to the stochastic modeling elements pj and i>j. However, the linear bounds L < v < U prevent an immediate application of Corollary 1 because realistic genome-scale metabolic models commonly enforce some fluxes with implied equalities, i.e. L j = U for some i. An example is the often called ATP maintenance flux, which is the rate at which ATP is lost due to non-growth-associated processes. These fixed fluxes are included to ensure that FBA models account for cellular demands that are not explicitly part of the genome-scale metabolic reconstructed model39. The continuity results depend on scaling, and Lemma 1 and Corollary 1 (see Methods) do not provide the Lipschitz continuity for RAMP since scaled fixed fluxes do not remain feasible under the requirement that v = L = U .

Figure 1. Graphical representation of the RAMP metho d. A depiction of the difference between an FBA (static) equality (the strong dashed line segment through the origin) and its stochastic RAMP counterpart (the bounded, shaded region). The light dashed lines represent the infinite set of linear constraints added by Eqns (10) and (11).

The equalities imposed by the bounding constraints could be re-cast probabilistically. However, unlike the steady state equations of Sv = 0, these constraints cannot be guaranteed to be of the form required by Lemma 1. For instance, if vt is the flux for the reaction that removes ATP for non-growth-associated demands, then the equality vi = Lt would naturally transition to P^j — Mt < vt < Lt + M) > 1 — 2e, which would associate with an SOCP constraint of the form Li - Mt + IRvH < v < Li + Mt - p,v||. Since the signs of Li — M, and Li + M, would agree for small values of M ¡, which would be required to remain realistic, the requirement of a positive right-hand side in Lemma 1 can not be ensured after writing both inequalities in the correct form.

Nevertheless, if we adjust the variable bounds dependent on the perturbed stochastic elements instead of prescribing them independent of perturbation, then we can relax the bounds to ensure feasibility. Probabilistically this means that we can set the values of Mt for the variable bounds so that we satisfy them with probability 1.

Theorem 1. For a collection ofprobability vectors p, and scenario matrices Sj andfor the lower and upper bounds U and L, let T([(p., S) : i = 1, 2, ..., m}, L, U) be the nonempty set of feasible fluxes satisfying the constraints of the RAMP model in Eq. (9). Assuming that each pt + Ap j is a probability vector, we then have for each

Theorem 1 highlights that any questionable discontinuities that FBA might exhibit with regard to parametric update are due to the imposed linear equalities that are the byproduct of enforcing the biologically unrealistic assumption of a uniform steady state. The RAMP approach adjusts FBA's modeling paradigm to include its inherent stochastic nature, and in the process, RAMP ensures the continuity that would be expected of the feasible flux states as long as the imposed linear bounds can also be relaxed commensurate with the magnitude of the perturbation. In conclusion, any discontinuity of FBA that could be caused by minute model adjustments are the outcome of an overly rigid (linear) model.

We note that a recent discussion in the literature, see refs 40-42, suggests that the creation of biomass might be sensitive to small perturbations in the growth coefficients. Indeed, the authors of ref. 40 have argued that in some cases perturbations to the biomass equation are necessary to achieve growth. The continuity of Theorem 1 counters any such concern with regard to RAMP, as flux states for any growth equation would remain close to those with a perturbed equation. Specifically, a change in the creation of biomass is controlled by the amount of coefficient perturbation, and hence, biomass creation can not depend on small adjustments to the growth coefficients in RAMP.

v e T([(p., S) : i = 1, 2, ..., m}, L, U)

that there is a A > 0 such that

mini v — v'\\ < XT

r = max{|| Ap,|| + ||AS(|| + ||Ap(||| lAS^I},

v' e T({(pi + Ap., S{ + as,.) : i = 1,2, ..., m}, L - Are, U + Are) and A is independent of all || A£>^|.

Connection between FBA and RAMP solutions. A reasonable question is if the deterministic FBA model is a limiting case of RAMP as the stochastic elements diminish? Theorems 2 and 3 show that convergent RAMP solutions are indeed FBA solutions, and hence, they answer this question in the affirmative. They also show that interpreting FBA as a limiting RAMP model characterizes the possible random flux states among cells in an optimal growth, steady state culture. The convergence result of Theorem 2 assumes an interiority condition, which is tacit as long as the equalities imposed by the bounding constraints are relaxed to allow arbitrarily small adjustments, i.e. as long as Li = vi = Ui is replaced with Li — n < vi < Ut + n for any arbitrarily small n > 0. The proof of Theorem 2 (see Methods) is straightforward and follows from the fact that if the primal and dual solutions converge as the stochastic elements dissipate, then the resulting necessary and sufficient conditions are that of FBA. However, the proof importantly identifies that not all random elements need to disappear, a point that prompts Theorem 3.

Theorem 2. Let pf, Si, Mp be sequences such that for each t the corresponding RAMP model satisfies Slater's interiority condition. Let v1 be an optimal solution of the RAMP model corresponding to pf, Si, and Mt, and assume V^ v as t ^ <. Assume likewise that a corresponding dual sequence of optimal solutions converges. Further assume that as t ^ < we have (pf )TSj ^ Si, Rj ^ 0, and Mt ^ 0. Then v is a solution to the FBA model

max{vGrowth: Sv = 0,L < v < M},

where S is the matrix whose i-th row is Si.

The assumption that Rf ^ 0 in Theorem 2 can be relaxed since the proof remains valid as long as the limiting matrices Rj satisfy Riv = 0 and RT(y + wi) = 0 (see the proof in Methods). Since Ri can have low rank, e.g. the rank of R1 in Eq. (12) is 1, we see that Ri need not generally vanish. This observation suggests that not all uncertainty needs to be removed to recover an FBA solution with RAMP. It is this interesting observation that motivates our third result.

Allowed probabilistic variation for optimal flux solution. Suppose v is a FBA solution corresponding with an average stoichiometric matrix, i.e. the rows of the matrix are p S¡. We pose the question of whether or not it is possible to identify nontrivial scenarios and probabilities so that v is a solution to both an FBA problem and to a limiting stochastic RAMP counterpart. If such probabilities and scenarios exist, then the optimal flux state v would be optimal for a range of stochastic variation within an optimal growth, steady state culture. We describe scenarios S{ as biologically possible for v if there exists (positive) probability vectors pt such that v maximizes cellular growth under the conditions that v satisfies L < v < U and

lim P(-M, < pTS,v < M) < 1 - 2e, V i.

Theorem 3 characterizes the biologically possible scenarios of any FBA solution. The argument arranges and partitions v into (v', v"), where v' ^ 0 for allj and v" = 0. This way Sv = S'v', where S' is the submatrix ofS whose columns correspond with v'.

Theorem 3. Let V = (V, 0), with V' ^ 0 for all j, be a solution to the FBA problem

max{vGrowth: pTSiV = ° V u L < v < U}. (13)

Then the scenarios of Si are biologically possible for V if and only if for all i we have Si v' = ate for some scalar ai ^ 0.

In biological terms, Theorem 3 identifies the probabilistic variations that are possible in the stoichiometric matrix for cells with the same optimal, steady state fluxes. Thus, if we trust that FBA solutions do indeed identify the limiting behavior of cells as they evolve to optimize to their environmental resources, then Theorem 3 concludes that optimized cells with the same flux state can be described by different probabilistic models.

Additionally, since Eq. (17) in the proof of Theorem 3 only requires that pt be a probability vector with nonzero entries, we can use this fact to construct a test that will check if any particular collection of scenarios is possible for optimal states. For example, suppose we are interested in investigating if some of the stoichiometric coefficients can vary for the cells having a specified flux state V in an optimized culture. Any collection of scenarios for which SjV does not have identical components is impossible, since there is no possible way to assign probabilities that make the flux state optimal once probabilistic variation is included in the model. If our mathematical models accurately assess biological reality, then we can claim that it is biologically impossible to have a collection of optimized cells with a common flux state that vary according to the suggested scenarios.

RAMP implemented as a computational model. We now compare the computational ability of several RAMP implementations with FBA. The computational results herein only initiate the broad comparisons that would benchmark RAMP against the numerous FBA adaptations. The metrics we use in our benchmarking are direct comparisons between FBA and RAMP predictions of experimentally determined fluxes, and RAMP's and FBA's ability to identify essential and viable genes. The latter since this is one of the commonly used tests to assess a metabolic model's potential for producing biologically relevant predictions. All flux tests were conducted on the iJO1366 metabolic model of Escherichia coli2S, whereas knockouts also used the iAF1260 metabolic model24. These high-quality reconstructed models have, respectively, 1366 and 1260 genes, 2583 and 2382 metabolic reactions, and 1805 and 1668 metabolites.

Figure 2. Plot of multipliers and growth coefficients. The largest permissible multiplier a (red) sorted in descending order for each of the 72 growth coefficients in E. coli model iJO136625. The absolute value of the stoichiometric growth coefficient is shown for comparison (black). Note that the multiplier values for 18 indices (see Table 2 for names) have an undetermined maximal value > 1026 and are placed at the vertical axis (104) for completeness.

2fe2s 4fe4s cl cobalt2 fe2 fe3

h h2o k mg2 mn2 mobd

nh4 ni2 pi so4 val-L zn2

Table 2. Tabulation of model-specific names of the (cytosolic) coefficients in the biomass equation of iJO136625 associated with numerically unbounded multipliers a, thus allowing RAMP to predict gene essentiality identical to that of standard FBA.

Gene essentiality is decided by simulating gene knockouts with a corresponding removal of the associated set of metabolic reactions2' 3. If the modified metabolic model affects the optimal growth rate, typically by at least a 50% reduction, then the gene is considered to be essential to the organism. When a predictive model correctly identifies a gene as essential, we label this as a true positive. Similarly, a true negative indicates a correctly identified non-essential gene by the predictive model. A false positive occurs if the analysis incorrectly labels a gene as essential, and a false negative if the analysis incorrectly labels a gene as non-essential. The predictive power of the model is the ratio of the sum of true positives and negatives to the number of genes.

The forthcoming computational experiments introduces various stochastic assumptions by querying either different percentiles 61-s or different scenarios for the growth equation. At this point, however, it is necessary to emphasize that the values of M¡ bounds are also part of the stochastic interpretation since they define the permissible variation from steady state. Specifically, the set of feasible flux states either grows or stays the same with an increase in any M, from which we know that the growth rate is non-decreasing as a function of the individual M¡ values. So if the M¡ values are too small, then we will likely under-approximate the growth rate, but if the M¡ values are too large, then we will likely over-approximate the growth rate. We solve an optimization problem (see Methods) to appropriately set the values of M¡ so that RAMP's growth rate exactly matches that of FBA's. The problem ensures that the sum of permissible variations from steady state is as small as possible.

RAMP assessment ofuncertainty in growth reaction. The goal of our first computational experiment is to gauge the relationship between RAMP and FBA for individual uncertainties in the growth reaction. This is motivated by the observation that stoichiometric coefficients of the growth reaction are empirically determined non-integers, with values spanning many orders of magnitude (see e.g. Table 1). We address this goal by posing the question, What is the maximum level of uncertainty tolerated in a single stoichiometric coefficient in the growth reaction before RAMP's predictions of gene-knockout essentiality depart from those of FBA in a given environment?

We tackle this challenge by individually increasing the level of randomness in a single growth coefficient, multiplying it with a factor containing ±a based on five scenarios, and we increase the multiplier a from zero in an iterative fashion (see Methods for details). As expected, RAMP and FBA coincide for the case with no randomness in the growth stoichiometric coefficients, i.e. with a = 0. However, as a increases it is possible that RAMP will find a different set of essential genes than FBA. Consequently, we determine the maximal value of a for which RAMP and FBA return identical sets of essential genes.

Assuming that the genome-scale metabolic reconstruction is a high-quality rendering of the biochemical processes associated with the individual biomass components, we make the following observations: If RAMP returns identical essential gene sets as FBA (i.e. stable predictive ability) for large percent variations of a single coefficient, then the cells of the culture could possibly be disparate in how they incorporate the associated metabolite in biomass. If the predictive ability instead degrades with only slight deviations, then the growth coefficient is more likely exhibiting little variation across the culture.

Experimental

Essential Nonessential

FBA 171(160) FBA 44 (35)

RAMPj 171(160) RAMP1 44 (35)

ramp2 162 (160) ramp2 43 (35)

Essential True Positive RAMP3o 166(155) RAMP3o 44 (35) False Positive

RAMP3b 166(152) RAMP3b 44 (30)

RAMP3c 163(152) RAMP3c 39 (30)

Computational ramp4 171(160) ramp4 44 (35)

FBA 77 (78) FBA 1074 (987)

RAMP1 77 (78) RAMP1 1074 (987)

ramp2 86 (78) ramp2 1075 (987)

Nonessential False Negative RAMP3o 82 (83) RAMP3o 1074 (987) True Negative

RAMP3b 82 (86) RAMP3b 1074 (992)

RAMP3c 85 (86) RAMP3c 1079 (992)

ramp4 77 (78) ramp4 1074 (987)

Table 3. Predictability of gene knockouts using stochastic models 1, 2, 3, and 4 for RAMP compared with FBA. RAMP models 3a, 3b, and 3c have increasing levels of stochastic variation, respectively. Results for the iAF1260 metabolic model are in parentheses, all other results are for the iJO1366 model.

Figure 2 displays, on a logarithmic scale, the calculated maximal values of the multiplier a before RAMP departs from FBA's gene-knockout predictions, for each of the 72 growth coefficients in iJO136625. We find that for 18 of the coefficients (see Table 2 for their names), we were only able to determine a lower bound on the value of the multiplier of a > 1026. Upon inspection of the numerical results for these 18 coefficients, we noticed the following two properties: within numerical accuracy, (i) the value of the mean (see Eq. (6)) is ¡i = 0, and (ii) the standard deviation (see Eq. (7)) scales with Mi. A direct consequence of these two observations is that the RAMP-specific constraints of Eq. (8) reduces to the standard FBA constraint Siv = 0. Hence, for these 18 coefficients of the biomass equation, the RAMP problem reduced to that of FBA.

The biochemical explanation for why RAMP collapses to the standard FBA problem for all values of a is immediate for the coefficients in Table 2 with the exception of 2fe2s, 4fe4s, fe3, and mobd: Upon inspection of iJO1366, we find that the biomass components are (i) freely available in the environment, and (ii) may be imported into the cytosol at no energetic or metabolic burden. Hence, any change of their stoichiometric coefficients in the biomass equation leaves the remainder of the biochemical network flux-patterns unchanged. However, the introduction or generation of 2fe2s, 4fe4s, fe3, and mobd into the cytosolic compartment is associated with metabolic burden through more complex pathways in the biochemical network reconstruction. For these four cases we were unable to determine a clear pattern, making a detailed biochemical explanation difficult.

We find for the remaining 54 biomass coefficients that the maximally allowed multiplier spans a surprising 4 orders of magnitude. Upon investigating a possible relationship between these maximal multiplier values and the corresponding stoichiometric coefficients, we find that they are uncorrelated (p = —0.004). Note that, for all of the 54 coefficients, the numerically determined RAMP solutions are more permissive than the FBA solutions, meaning that at the determined maximal multiplier values, RAMP identified as essential a subset of the FBA-determined essential genes. In total, we found that 23 of the 54 coefficients supported maximal permissible multipliers a > 1. This is quite surprising, as it suggests that both the sign and the magnitude of the corresponding stoichiometric coefficient may vary while the growth rate and RAMP's predictive power are maintained.

RAMP computational predictions of gene essentiality. The essential gene predictions of RAMP and FBA for our second computational experiment, in which all coefficients in the growth reaction are simultaneously uncertain (see Methods), are tabulated in Table 3. We assumed in all experiments an oxygen-unlimited environment with glucose as the sole, and limiting, carbon source as the default environment. We additionally conducted the same tests with glycerol as the only carbon source, finding results that mirror those of glucose (not shown).

The tallies in Table 3 show some interesting differences between the two genome-scale metabolic reconstructions iAF126024 and iJO136625 in response to the different RAMP implementations. While the majority of biochemical reactions are the same in these two models, the RAMP analysis emphasizes that small details of the reconstructed model affect the identification of essential genes. Here, RAMP may be an additional approach in vetting the quality of model reconstructions.

Table 3 shows that the default RAMP model (RAMPj), which only permits uncertainty past the stated accuracy of FBA's growth equation, precisely replicates the essential gene predictions of FBA for both genome-scale metabolic model reconstructions. This fact computationally substantiates Theorem 2 and numerically demonstrates further that minor probabilistic variations mimic FBA.

Model 2 (RAMP2) is designed to magnify the perturbations of Model 1 by a multiplicative factor p and 'spread' the scenarios in the last few digits. From Table 3, we see that this approach impacts RAMP's predictive ability. For the iJO1366 model, all values of p decreased the true positives by 9 and increased the true negatives by 1, which means that RAMP2 had 8 fewer correct predictions than FBA with increased scenario variability as magnified by

20 15 10 5 0 -5

20 15 10

_ I I | I I | I I | I I | I I | I I | I I | I I | I r_

E- (a) *

"i i I i i I i i I i i I i i I i i I i i I i i I i i"

M-11111111111111111111 'Al '-H (c)

□ —

¡** *

b'111111 'SA' 1111111111111 =-<b)

□ RAMP

• EXP

r". I ■ ■ I ■ ■ I ■ ■ I ■ ■ I ■ ■ I ■ ■ I-I-. I .73

LL1 | 1 1 | 1 1 |1 1 | 11 | 11 | 1 1 | 1 1 |"JJ

1 4 7 10 13 16 19 22 25 28 1 4 7 10 13 16 19 22 25 28 Reaction index Reaction index

Figure 3. Comparison of experimental fluxes with RAMP and FBA predictions. Experimentally determined fluxes in E. coli (black circles with whiskers) for 28 reactions in the central carbon metabolism are compared with the predictions from FBA (blue stars) and from RAMP3 with a = 0.2 (red squares). All calculations are made with model iJO136625. The panels correspond to experiments conducted under conditions of (a) aerobic batch growth26, (b) anaerobic batch growth26, (c) carbon-limited chemostat at dilution rate 0.1/h27, and (d) carbon-limited chemostat at dilution rate 0.4/h27. The correspondence between reaction index and biochemical reaction is given in Supplementary Table S1.

p. In contrast, the iAF1260 model agreed with FBA for all values of p. RAMP's predictive power compared with FBA for RAMP2 decreased from 91.14% to 90.56% with the iJO1366 model and matched FBA's 91.32% with the iAF1260 model. The similarity between FBA and RAMP for RAMPj and RAMP2 shows that probabilistic deviations in the growth coefficients beyond what is reported by the genome-scale metabolic reconstructions only lead to minor adjustments in predictive power.

Scaling the scenarios proportionately (RAMP3a, RAMP36, and RAMP3c representing increasing levels of variation) has a slightly more varied effect than multiplying the default variations by p. Results with a > 0.4 had little value since these models relaxed the FBA constraints to the point at which few, if any, essential genes were identified, and for this reason we only report the results for a < 0.3. The trend as a increases is that RAMP3 reduces the number of predicted essential genes, which lowers the number of true and false positives but increases the number of true and false negatives. Five true positives move to false negatives in the iJO1366 model for RAMP3a and RAMP36, which reduces RAMP's predictive ability to 90.78%. However, for larger a, RAMP3c increases its number of true negatives and reduces the number of false positives, which increased predictive ability to 90.92%. The trend for the iAF1260 model is similar.

Altering the parameter e (RAMP4) had no effect on RAMP's ability to predict essential genes, and for all tested values of e, RAMP agreed with FBA. This clearly suggests that it is more important to tune RAMP with regard to the scenarios than it is with the probabilistic guarantee of satisfying the near equilibrium constraints in Eq. (3).

RAMP comparison with experimentally determined fluxes. To investigate RAMP's consistency with experimentally determined fluxes, we formulated a non-linear (quadratic) objective function for RAMP that minimizes the distance between a set of experimentally determined flux data and the corresponding reaction fluxes in a genome-scale metabolic model. We directly compare the results from RAMP3 with a = 0.2 with the corresponding results from FBA. The details of the problem formulation are given in Eqs (22) and (23). We used flux data6 measured for 28 reactions from the central carbon metabolism of E. coli growing on glucose in batch conditions26 both aerobically (Fig. 3(a)) and anaerobically (Fig. 3(b)), as well as in aerobic chemostat conditions27 at the two dilution rates of 0.1/h (Fig. 3(c)) and 0.4/h (Fig. 3d)).

A striking visual feature of Fig. 3 is that the RAMP-predicted fluxes consistently provide a better fit to the experimental fluxes than the FBA-predicted ones. In Table 4, we quantify the performance of RAMP and FBA in predicting fluxes by calculating the mean square error MSE = \\vpredicted — vexp\\2 /N, with N = 28. The results

RAMP MSE Relative MSE

Fig. 3(a) 0.171 5.8%

Fig. 3(b) 1.105 24.4%

Fig. 3(c) 3.116 31.7%

Fig. 3(d) 0.106 15.1%

Table 4. Quantification of RAMP and FBA ability to predict experimental fluxes. Mean square error (MSE, see text for definition) for RAMP, and relative MSE: the RAMP MSE divided by the FBA MSE.

reported in Fig. 3 and Table 4 demonstrate that RAMP significantly outperforms FBA in predicting experimental fluxes in a wide range of conditions.

Discussion

The expanding suite of systems-biology computational methods described as constraint-based analyses has significantly improved our ability to probe the properties of genome-scale metabolic reconstructions. However, the basic premises of these methods, the assumptions of steady state and no heterogeneity in a cellular population, has stayed unchanged. In this article, we have developed the conceptual framework for a new paradigm of constraint-based methods by explicitly relaxing the basic assumption of steady-state while allowing for heterogeneity. We have presented the mathematical foundation proving that, as levels of uncertainty are reduced, the RAMP method will become identical to standard FBA.

Building on the RAMP foundation, it will be possible to develop robust versions of many of the contemporary constraint-based methods. For some of the new extensions to FBA, the inclusion of RAMP formalism is straightforward, e.g. FBAwMC32, GIMME33, MOMENT34, and GX-FBA35. However, also methods that add gene-regulatory effects, such as rFBA43, and thermodynamic considerations44 may be directly implemented with RAMP, since many of these methods generate modified metabolic networks, or temporal sequences of such. Note that, since it is known that uncertainty in variable bounds, e.g. restrictions such as upper bounds on reaction fluxes, mathematically is not different from uncertainty in the stoichiometric coefficients20, the RAMP formulation is also capable of modelling uncertainty in flux constraints.

Theoretically, we are able to formulate a wide host of linear and nonlinear RAMP models as SOCPs. However, the nonlinear ones have not been uniformly translatable into computational reality. While we have been able to solve the L2 RAMP flux-minimization problem (see Eq. (22)) and show that RAMP outperforms FBA for flux predictions, we have achieved interspersed success with native SOCP solvers for general RAMP models. Specifically, we have attempted to

• solve L2 versions of Eq. (21),

• predict essential genes with a RAMP version of the minimization of metabolic adjustment (MOMA) method45, which is a quadratic model, and

• increase variation beyond the growth coefficients to replicate stochastic environmental bounds.

While we are able to solve these problems occasionally, we were unable to identify a choice of solver settings and convergence criteria to achieve consistent and stable success for the numerical implementation with current SOCP software for any of these extensions. Establishing a stable computational platform for the general RAMP paradigm is an important goal for future work.

The SOCP constraints of RAMP appropriately inject statistical variation into the FBA paradigm, but the conclusion of whether or not RAMP's stochastic perspective is more or less restrictive than FBA's deterministic setting remains unclear. From a mathematical point of view, RAMP is neither a direct relaxation of FBA nor a direct restriction of FBA. However, our computational results demonstrate that RAMP is, in practice, a relaxation of FBA. When investigating uncertainty in individual growth coefficients, we uniformly show that discrepancies in gene essentiality are due to RAMP being less restrictive, as evidenced by RAMP having predicted at least one essential gene less than FBA, i.e. a subset of the FBA results. Similarly with all growth coefficients assumed uncertain, RAMP's stochastic adaptation results in a model whose trend is to predict fewer essential genes, especially as uncertainty increases.

These computational outcomes favorably agrees with the biological sentiment that cellular metabolic networks are more robust in their wild-type setting and are increasingly fragile as they evolve toward an optimal, and specialized, steady-state. Hence, the number of essential genes should reduce as statistical variation increases, which is predicted by RAMP. Moreover, our mathematical analysis and computational outcomes show that RAMP mimics FBA if the stochastic elements are sufficiently small. In particular, if variation is not allowed beyond the assumed accuracy of an FBA model, then RAMP behaves like FBA for predictions of essential genes.

Theorem 3 provides a new approach to discriminate FBA's optimal flux states: For example, some FBA solutions may have no biologically possible scenarios while others may have many. This suggests that the former is less stable with regards to random variation, whereas the latter is optimal under a wider range of stochastic possibilities. One may intuit that actual flux states should be robust against random variation46, and hence, learning to calculate FBA solutions that can witness as much variability as possible is a promising direction of future research.

Finally, our choice of letting the objective function coincide with the growth reaction was purely a choice based on making RAMP easily comparable to traditional implementations of FBA, and RAMP is valid beyond this choice.

Methods

In this section we provide the mathematical foundation and proofs of the theorems presented in the main text. We also discuss aspects of the computational stability and implementation of RAMP.

Lemma 1

The following lemma establishes a general property from which our result follows.

Lemma 1. Let A be an n-element row vector, R be a q x n matrix, and b be a positive scalar. Let T(A, R) = {v : Av + ||Rv|| < b}. Then, for any v e F(A, R), there are scalars z and X satisfying 0 < z < 1 and X > 0 such that

min{||v - v'|| : v' e T(A + AA, R + AR)} < ||v - zv\\ < A(||AA || + ||AR||), where zv e T(A + AA, R + AR) and X is independent of\|AA|| and ||AR||.

Proof. Should v satisfy (A + AA )v + ||(R + AR)v|| < b, then we immediately have

min {||v - v'\\ : v' e T(A + AA, R + AR)} = ||v - v|| = 0,

and we are done with z = 1 and X = 0.

Suppose instead that (A + AA)v + ||(R + AR)v|| > b. From the Intermediate Value Theorem there is a t such that 0 < t < 1 and

0 < (A + (1 - t)AA)v + ||(R + (1 - t)AR)v|| = b. To ease notation, let A' = A + (1 - t)AA and R' = R + (1 - t)AR. Also select z such that

z =-=-.

b + r(||AA II + ||AR||)||v||

(A + AA )(zv) + ||(R + AR)(zv)|| - b = (A' + tAA )(zv) + ||(R' + rAR)(zv)|| - b

< z(A'v + ||R'v|| + t(||AA\\ + ||AR|) - b = z(b + r^lAA|| + ||AR|)|) - b = 0.

We conclude that (A + AA)(zv) + ||(R + AR)(zv^ < b, and hence,

zv e T(A + AA, R + AR).

Let A = || v |2 /b, from which we have that

IIv - zvll

b + rai AA || + ||AR|v||

ai AA\\ + ||AR|)

b + rai AA || + ||AR|| v|| l|v|f

AA || + ||AR|)

Hence,

A(|| AA\\ + ||AR|).

min^|v - v'|| : v' e T(A + AA, R + AR)} < ||v - zv\\ < A^| A^l + ||AR|),

and X is independent of 11 AA || and ||AR|. □

Lemma 1 shows that vectors satisfying SOCP constraints of the form Av + ||Rv|| < b, with b > 0, can be scaled to remain feasible, and moreover, that the magnitude of the adjustment to remain feasible is uniformly bounded by the magnitude of the perturbations. Importantly, this result immediately extends to finite collections of SOCP constraints of the same form, a fact we formalize in Corollary 1.

Corollary 1. For i = 1,..., m, let At be an n-element row vector, Rj be a q x n matrix, and bj be a positive scalar. Let

F({(Aj, Rj) : i = 1, 2, ..., m}) = {v : A, v + \\R^\\ < b, for i = 1, 2, ... m,}.

Then, for any v e F({(At, R,) : i = 1,2, ..., m}), there are scalars z and X satisfying0 < z < 1 and X > 0 such that

min{||v - v'\\ : v' e T({Ai + AA i, Ri + AR,, i = 1, 2, ..., m})} < II v - zv\\ < A max ,{(|| AA,|| + ||ARJ|)},

where zv e T({(Ai + AA,., Rj + AR,) : i = 1, 2, ..., m}) and A is independent of all HAA.H and | |AR,| |.

Proof. The proof follows directly from Lemma 1 and its proof upon letting z, and AI be the scalars for each I and setting z = min,{z,} and A = max,{ A,}. □

Unfortunately, systems of linear equalities like those of FBA do not satisfy similar continuity properties, and changing the stoichiometric matrix in FBA can lead to discontinuities. As a simple example, the system of homogeneous equations

v - v2 = 0 vl — (1 + a)v2 = 0

has substantially different solution sets around a = 0. Notice that vl = v2 = t solves the system for any t if a = 0, but that the only solution for a ^ 0 is vl = v2 = 0. Hence the solution vl = v2 = 1 at a = 0 is not arbitrarily close to a solution with a ^ 0.

Proof for Theorem 1

Proof. Let A = pT§i and Ai + AA,. = (p. + Ap.)T(S,. + AS,.). Then,

l|AA,l| = ||(pT + ApT)(S,. + AS> - pJSvW < \\pjAS,. + ApJSt + Ap.AS^||v||

< AS J + ||Ap,|| + ||AS^||ApJ|), (14)

where k,a = \\v\\max{||S,||, ||p,||, 1}. Further note that

1^.11 = ||(Ri. + AR() - R,||

= S—yPi + Ap.(I - e(p. + Ap.)T)(Si. + AS) - 4P(I - epT)Si = S1-£|| Q p + Ap - .JP, )(I - epj )Si+%/ p + Ap ((I - epT )AS, - eAp. S, - eAp. AS,.)]

< Sm^pTApt — -v/p^y||I - ep,T 111 Si-| |

+ ||I - ep. 11 x|| ASj| + Jq ||ApJ||S| + Jq\Ap\\| ASj|). (15)

From the Mean Value Theorem there is a jjiik between pk and pik + Ap,k such that

Jpik + Apk - -p

i^ik Hence,

Iivp+^p — V^n = ji( + Apik )2

where j, = min^^}. We now have from inequality Eq. (15) that

||AR,|| < «iRa|ASM + l|Ap,l| + IIAS^I|ApJ|), (16)

k,R = &1_s maxl-^ai I - epj\\ + 4q)\\Si||, ||I - epj\\, Vqi.

Each of RAMP's SOCP constraints may be re-written as

(A,. + AA,.)v + ||(R,. + AR,)v|| < M, i = 1, 2, ..., m and -(A,. + AA,.)v + ||(R,. + AR>|| < M, i = 1, 2, ..., m.

From Corollary 1 and inequalities Eqs (14) and (16) we know that there is a v satisfying these perturbed SOCP constraints and a A > 0 defined independent of|| AA,|| and 1| AR,||, and subsequently independent of 1| AS,||, such that

||v - v'll < Amax{||AA¡|| + ||AR¡||}

< W{(«¡A + k¡r)(||AS¡|| + ||ApJ + ||AS¡||||Ap¡||)}

< Amax{K¡A + K¡R}max{(||ASJ + ||Ap.|| + ||AS¡||||Ap.||)}

ce A, k¡a, and kr are

satisfy

Since A, and kr are all independent of 11AS,||, the proof is complete upon noticing that v' must also

A re < v' < U + Are,

r = max{|| Ap l| + ||AS¡|| + ||Ap||| |AS¡|

Proof for Theorem 2

Proof. For each t let (y*, y*, w*, w*, p*, a*) be the assumed dual solution that converges to (y, y, w, w, p, a). Further let f be a vector so that fTv = vGrowth. The dual can be assumed to satisfy Slater's interiority condition since the necessary and sufficient primal-dual conditions in this case would be satisfied by the convergent sequences:

-M* + \\R jvt\\ < (pp*)TS¡vt < Mt - ||R,V||, V ¡

L < v* <

^((S!)Tp¡ (y* - wt) - (Rt)T (y + w¡)) + p - o * = f

\\y!\\ < y*, v ¡ \\w¡W < w¡, V ¡

p*, o * > 0

YmMí($¡ + W) + Utp - LTo* - fTv* = 0.

Dual feasibility is defined by the 3rd through 6th constraints, from which we can see that the dual is always strictly feasible, and hence, the dual satisfies Slater's interiority condition. Select any y and w variables such that \\y \\ < y* and 11 w¡|| < wj. Since p* — a* can attain any vector, these two variables can be chosen so that p* > 0, a* > 0 and that

E((s",t)Tp,t y - */) - (R/)T y + w *)) + p - o * = f.

Hence the system is indeed necessary and sufficient for optimality.

Allowing * ^ <, we have by assumption that (p*) S¡ ^ S, R¡ ^ 0, and M * —> 0, from which we have

Sv = 0 L < v < U

S (y - w) + p - o = f

y, w, p, o > 0

UTp - Lto - fTv = 0.

Since these are the necessary and sufficient conditions for FBA, v is an optimal solution to the FBA model. Proof for Theorem 3

Proof. Let v = (v', 0) be as stated. Then,

R¡v = R' v' = JP¡ (I - epT )§' v' = 0 if and only if epj (§' v') = (S' v'). (17)

The last equality shows that R¡v = 0 if and only if S'v' is an eigenvector of epT for the eigenvalue 1. Since p¡ is a probability vector, ep T has only two eigenspaces, one of dimension q — 1 for the eigenvalue 0, and one of dimension 1 for the eigenvalue 1. All eigenvectors for the eigenvalue of 1 are scalar multiples of the all ones vector e. Hence, R¡v = 0 if and only if the scenarios for each I satisfy S¡ v' = a¡e for some aI ^ 0. RAMP's SOCP constraints as M¡ j 0 are

||R¡v|| < prS,tv < - ||R¡v|| o J p¡ S¡v = 0 ' Rv = 0.

where A = Amax.K

Hence, the limiting RAMP model as Mi [ 0 is the linear program

max{vGrowth : p^S^ = 0, Riv = 0, V i, L < v < U}. The necessary and sufficient conditions are

p.Sjv = 0, V i

Riv = 0, V i L < v < U

S' PiX

- P + ^ = f

P, <y >

UTa - LTp - fTv :

where f is the vector so that fTv = vGrowth.

Suppose for each I that S'.v' = ae for some aI ^ 0. Then v satisfies R,v = 0 for all i. Moreover, since v solves Eq. (13), the strong duality theorem of linear programming guarantees a solution to Eq. (19) with v = v and Wj = 0. Alternatively, if for some I we have S'¡v' ^ a e for all nonzero aI, then Rvv ^ 0 and v is infeasible in Eq. (18). Hence, in this case S{ is not biologically possible for v. □

Computational stability of the RAMP framework. A straightforward implementation of RAMP (see Eq. (9)) as a computational model proved elusive. While we were able to determine optimal solutions in many situations for the tested genome-scale metabolic reconstructions, we were unable to find consistent achievement of optimality for the general RAMP problem with a wide variety of state-of-the-art algorithms to solve SOCPs. The computational challenge is likely caused by two characteristics specific to genome-scale metabolic reconstructed networks. First, a substantial number of variables are unsigned and, essentially, unbounded. Second, the FBA problems arising from genome-scale reconstructed metabolic networks are highly degenerate (a detailed discussion of this point is provided in ref. 47) and often have high-dimensional solution sets. These characteristics combine in a way that seems to hamper the underlying interior-point algorithms employed in native SOCP solvers. We also tried to use standard nonlinear reformulations that allowed us to experiment with nonlinear solvers, but again success was tepid at best.

Motivated by biological considerations, we therefore implemented a limited version of RAMP as a computational model. If the probabilistic variability of each constraint is restricted to a single coefficient, then the associated SOCP is linear. The linear constraint is constructed by, for example, assuming that the last coefficient of the i-th constraint is the sole random variable. Then for some q-vector s we have

IIRivll

epT )S v||

,(m-1)'

e, s]v||

epj e), (I

sfii (I

^¡P,(I - epT)[Si ,e, ^¡P, [s,. i(e - epTe), dP[0,..., 0, (I - epT)s]v|| E4Pt(I - epT)s\\vn.

This calculation shows that if we restrict probabilistic variability to, for example, the coefficients of the growth reaction, then each of the SOCP constraints of the form

can be re-written linearly as

pJSvv -pTSv -

IRH - Mt < plStv < Mt - ||R,v||,

Wl-dP(I - epT)S,.Growth«vGrowth > - M, I

^ Pl-^ (I - epT )Si, Growth I vGrowth < M I

where §j Growth is the column of S{ containing the scenarios for the growth coefficient.

We choose to use the growth coefficients for two reasons. First, we have already noted that this is an empirically derived reaction that is based on the aggregate properties of a cell culture, and thus, the stoichiometric values are inherently associated with uncertainty. Furthermore, it is well known that biomass composition (and thus the values of the stoichiometric coefficients of the growth reaction) of a cellular culture is affected by nutrient conditions and the culture's growth phase48, 49.

The linear re-formulation in Eq. (20) allows RAMP to be solved with standard linear simplex solvers, which proved to be computationally stable. All numerical work was conducted with the freeware GLPK or with the commercial solver Gurobi©, both of which worked well with the COBRA toolbox37.

RAMP prediction of gene essentiality. We present two computational experiments to assist in assessing how stochastic variation affects RAMP's ability to predict gene essentiality. The first experiment iteratively induces individual randomness in each of the growth coefficients, holding the others at their nominal value as

FBA Growth Scenario

Metabolite Coefficient 1 2 3 4 5

2ohph[c] -0.000223 -0.0002232 -0.0002231 -0.000223 -0.0002229 -0.0002228

adp[c] 53.95 53.948 53.949 53.95 53.951 53.952

probability 0.0351 0.2389 0.4520 0.2389 0.0351

Table 5. The first column is the metabolite name from the iJO1366 model, and the second column contains the associated growth coefficient from the FBA model. The sign indicates whether the metabolite is an input (negative) or output of the growth reaction. The remaining columns are the scenarios for the default case along with their probabilities.

stated in the model. The second experiment assumes simultaneous randomness in all growth coefficients. In all of the computational experiments, we used the E. coli genome-scale metabolic models iAF126024 and iJO136625 in a minimal nutrient environment consisting of (i) the models' presets to emulate M9 medium, (ii) a maximal glucose uptake rate of 10 mmol/gDw/h, and (iii) unlimited oxygen uptake.

RAMP variation in a single growth coefficient. Here, we increase the level of uncertainty in a single growth coefficient in an iterative fashion, and we systematically go through all of the coefficients. The goal of this experiment is to gauge predictive ability against individual uncertainties in the growth equation, and we measure the predictive ability simply by comparing the number of essential genes with that predicted by FBA. If the predictive ability is stable for large percent variations of a single coefficient, then the cells of the culture could possibly be disparate in how they use the associated metabolite. If instead the predictive ability degrades with slight deviations, then the growth coefficient is more likely conserved across the culture.

The experiment uses q = 5 scenarios for each coefficient. The two most extreme scenarios multiply the coefficient by ±a with probability P(3/2 < z) = 0.0351, where z is standard normal variable. Two other scenarios multiply the coefficient by ±a/2 with probability P(1/2 < z < 3/2) = 0.2389, and the other scenario leaves the coefficient unchanged with probability P(-1/2 < z < 1/2) = 0.4520. The value of 6l-s was 3, which meant that the probabilistic constraints were guaranteed with probability 0.9997.

The largest value of a was calculated within an accuracy of0.0001 and with a final value that allowed RAMP's predictive accuracy to match that of FBA's. Thus, starting from an initial value of a = 0.0001, we (i) increased a by a factor 2 until RAMP's gene knockout predictions no longer matched that of FBA, and (ii) at this point, we initiated a binary search to identify the maximal value of the multiplier a, which resulted in RAMP and FBA gene knockout predictions matching.

RAMP simultaneous variation in allgrowth coefficients. In this experiment we assumed simultaneous randomness in all growth coefficients. Four probabilistic models were used to assess RAMP's overall sensitivity to stochastic variation, and each was compared against FBA's gene knockout predictions as calculated by the COBRA toolbox (see Table 3 for results).

Model l(default) Our default RAMP model assumes that the means of the growth coefficients are the values stated in the FBA model and that probabilistic variation is restricted to the first unspecified significant digit. Assuming 10-di identifies this digit for the i-th growth coefficient, we further assume that each growth coefficient has the 5 scenarios in which it is perturbed by ±r\ ■ 10-di, where n is one of 0, 1, or 2. The probability of the scenario with no perturbation is P(-1/2 < z < 1/2) = 0.4520, with perturbation ±10-di is P(1/2 < z < 3/2) = 0.2389, and with perturbation ±2 • 10-di is P(3/2 < z) = 0.0351, where z is a standard normal variable. Two examples of these scenarios are illustrated in Table 5. We choose to set 6l-£ = 3, so that the probabilistic constraints have a 99.97% guarantee of satisfaction.

Model 2 The second probabilistic model multiplies n in model 1 by the scalar p to assess how RAMP solutions adjust as the scenarios deviate from those of FBA. Eight tests with p = 2, 3,..., 9 were considered. Integers beyond

9 resulted in sign changes and were not considered. The probability scenarios from Model 1 were used, as well as Sl-s similarly being set to 3.

Model 3 Since Model 2's scaling by p disproportionately effects small growth coefficients, we compensate by replacing the scenarios of Model 1 with percentages of the growth coefficient itself. We choose different ranges for a, using the scenarios (1 ± a ■ n) ■ §j Growth, where n is one of 0, ±1/2, or ±1. The probabilities are unchanged from those in Model 1, and 6l-s is 3.

Model 4 To assess how RAMP reacts to changes in the certainty of satisfying the probabilistic constraints, model 4 changes 6l-s. All other model parameters are inherited from the default model.

In our testing of the predictive ability of these four computational RAMP models (for computational results, see Table 3), we used the short-hand notation given in Table 6 for our choice of combination of model and parameters.

Determination of bounds M.. The selection of Mt (see Eq. (9)) is an important parameter to decide. However, instead of imposing these bounds arbitrarily, these parameters are determined so that each RAMP model accurately returns the targeted, optimal growth rate. Let 7* be the optimal growth rate as calculated by the FBA model. The Mj values are calculated by solving the following optimization problem.

Model 1 RAMP1

Model 2 ramp2

Model 3 RAMP3„ 0.001 < a < 0.01

RAMP3i a = 0.02

RAMP3c 0.03 < a < 0.3

Model 4 ramp4 0.01 < 2s < 0.9

Table 6. The choice of RAMP parameters for models 1-4 in our computational simulations.

min ||M||j subject to

vGrowth > Y *

P,v|| - M¡ < pTS¡v < M¡ - p¡v||, L < v < U M > 0,

where M is the vector whose i-th component is M, and 11 • ||¡ is the L1 norm. Setting M to be the calculated optimal solution of Eq. (21) tightens the SOCP constraints while ensuring that the optimal growth rate is held at its desired value. We experimented with L2 and counterparts; however, the L2 norm suffered from inconsistent solves, and the norm overly relaxed constraints, which was not surprising.

Comparison of RAMP and FBA with experimentally determined fluxes. We tested the ability of FBA and RAMP to agree with experimentally determined fluxes. For RAMP we solved,

min \\v - vEXp\I2 subject to

vGrowth > 9 Y *

\\RjV\\ - M¡ < pTS¡v < M¡ - p,v||, i = 1, 2, ..., m

L < v < U. (22)

where vEXP is the set of experimentally determined fluxes. The M¡ values were calculated to achieve the FBA-optimal growth rate as noted in the previous section. The scenarios were those of Model 3 (RAMP3) with a 20% scaling of the growth equation, i.e. a = 0.2. The (FBA) optimal growth rate 7* was scaled down by 6 = 0.9 because the experimental fluxes were not guaranteed to coincide with optimized cultures. Note that Eq. (22) is an example of a non-linear (quadratic) objective within the RAMP formalism. Consequently, we determine a set of RAMP flux values that are as close to the experimental flux set as possible.

Likewise, the ability of FBA to realize the experimental fluxes was tested by solving the following adaptation of Eq. (2),

min llv - vexp\\ subject to

vGrowth > 9 Y *

Scv = 0 Lc < v < Uc.

As with RAMP, d was set to 0.9, and the solution was again as close as possible to vEXP..

References

1. Bordbar, A., Monk, J. M., King, Z. A. & Palsson, B. O. Constraint-based models predict metabolic and associated cellular functions. Nat Rev Genet 15, 107-120 (2014).

2. Palsson, B. O. Systems Biology: Constraint-based Reconstruction and Analysis (Cambridge Univ Press, 2015).

3. Maranas, C. D. & Zomorrodi, A. R. Optimization Methods in Metabolic Networks (Wiley, 2016).

4. Ibarra, R. U., Edwards, J. S. & Palsson, B. O. Escherichia coli k-12 undergoes adaptive evolution to achieve in silico predicted optimal growth. Nature 420, 186-189 (2002).

5. Dekel, E. & Alon, U. Optimality and evolutionary tuning of the expression level of a protein. Nature 436, 588-592 (2005).

6. Schuetz, R., Kuepfer, L. & Sauer, U. Systematic evaluation of objective functions for predicting intracellular fluxes in escherichia coli. Mol Syst Biol 3, 119 (2007).

7. Schuetz, R., Zamboni, N., Zampieri, M., Heinemann, M. & Sauer, U. Multidimensional optimality of microbial metabolism. Science 336 (2012).

8. Varma, A., Boesch, B. W. & Palsson, B. O. Biochemical production capabilities of escherichia coli. Biotechnol Bioeng42, 59-73 (1993).

9. Lewis, N. E., Nagarajan, H. & Palsson, B. O. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nat RevMicrob 10, 291-305 (2012).

10. Fong, S. S. & Palsson, B. O. Metabolic gene-deletion strains of escherichia coli evolve to computationally predicted growth phenotypes. Nat Genet 36, 1056-1058 (2004).

11. Burgard, A. P., Pharkya, P. & Maranas, C. D. Optknock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng 84, 647-657 (2003).

12. Burgard, A. P., Vaidyaraman, S. & Maranas, C. D. Minimal reaction sets for escherichia coli metabolism under different growth requirements and uptake environments. Biotechnology Progress 17, 791-797 (2001).

13. Almaas, E., Kovacs, B., Vicsek, T., Oltvai, Z. N. & Barabasi, A. L. Global organization of metabolic fluxes in the bacterium escherichia coli. Nature 427, 839-843 (2004).

14. Almaas, E., Oltvai, Z. N. & Barabasi, A.-L. The activity reaction core and plasticity of metabolic networks. PLoS Comput Biol 1, e68 (2005).

15. Lee, K. H., Park, J. H., Kim, T. Y., Kim, H. U. & Lee, S. Y. Systems metabolic engineering of escherichia coli for l-threonine production. Mol Syst Biol 3, 149 (2007).

16. Park, J. H., Lee, K. H., Kim, T. Y. & Lee, S. Y. Metabolic engineering of escherichia coli for the production of l-valine based on transcriptome analysis and in silico gene knockout simulation. Proc Natl Acad Sci USA 104, 7797-7802 (2007).

17. Yim, H. et al. Metabolic engineering of escherichia coli for direct production of 1,4-butanediol. Nat Chem Biol 7, 445-452 (2011).

18. Khodayari, A., Chowdhury, A. & Maranas, C. D. Succinate overproduction: A case study of computational strain design using a comprehensive escherichia coli kinetic model. Front Bioeng Biotechnol 2, 76 (2014).

19. Yizhak, K. et al. A computational study of the warburg effect identifies metabolic targets inhibiting cancer migration. Mol Syst Biol 10, 744 (2014).

20. Ben-Tal, A., El Ghaoui, L. & Nemirovski, A. Robust Optimization (Princeton Univ. Press, 2009).

21. Ko, E. P., Yomo, T. & Urabe, I. Dynamic clustering of bacterial population. Physica D 75, 81-88 (1994).

22. Avery, S. V. Cell individuality: the bistability of competence development. Trends in Microbiology 13, 459-462 (2005).

23. Smits, W. K., Kuipers, O. P. & Veening, J.-W. Phenotypic variation in bacteria: the role of feedback regulation. Nat Rev Micro 4, 259-271 (2006).

24. Feist, A. M. et al. A genome-scale metabolic reconstruction for escherichia coli k-12 mg1655 that accounts for 1260 orfs and thermodynamic information. Mol Syst Biol 3, 121 (2007).

25. Orth, J. D. et al. A comprehensive genome-scale reconstruction of escherichia coli metabolism-2011. Mol Syst Biol 7, 535 (2011).

26. Perrenoud, A. & Sauer, U. Impact of global transcriptional regulation by arca, arcb, cra, crp, cya, fnr, and mlc on glucose catabolism in escherichia coli. J Bacterial 187, 3171-3179 (2005).

27. Nanchen, A., Schicker, A. & Sauer, U. Nonlinear dependency of intracellular fluxes on growth rate in miniaturized continuous cultures of escherichia coli. Appl Environ Microbiol 72, 1164-1172 (2006).

28. Ben-Tal, A. & Nemirovski, A. Robust convex optimization. Mathematics of Operations Research 23, 769-805 (1998).

29. Ben-Tal, A. & Nemirovski, A. Robust solutions of uncertain linear programs. Operations Research Letters 25, 1-13 (1999).

30. Ben-Tal, A. & Nemirovski, A. Robust truss topology design via semidefinite programming. SIAM journal on optimization 7, 991-1016 (1997).

31. Zavlanos, M. & Julius, A. Robust flux balance analysis of metabolic networks. In American Control Conference (ACC), 2011, 2915-2920 (2011).

32. Beg, Q. K. et al. Intracellular crowding defines the mode and sequence of substrate uptake by escherichia coli and constrains its metabolic activity. Proc Natl Acad Sci USA 104, 12663-12668 (2007).

33. Becker, S. A. & Palsson, B. O. Context-specific metabolic networks are consistent with experiments. PLoS Comput Biol 4, e1000082 (2008).

34. Adadi, R., Volkmer, B., Milo, R., Heinemann, M. & Shlomi, T. Prediction of microbial growth rate versus biomass yield by a metabolic network with kinetic parameters. PLoS Comput Biol 8, e1002575 (2012).

35. Navid, A. & Almaas, E. Genome-level transcription data of yersinia pestis analyzed with a new metabolic constraint-based approach. BMC Systems Biology 6, 150 (2012).

36. Reed, J. L., Vo, T. D., Schilling, C. H. & Palsson, B. O. An expanded genome-scale model of escherichia coli k-12 (ijr904 gsm/gpr). Genome Biol4, R54 (2003).

37. Schellenberger, J. et al. Quantitative prediction of cellular metabolism with constraint-based models: the cobra toolbox v2.0. Nat.

Protocols 6, 1290-1307 (2011).

38. Chu, M., Zinchenko, Y., Henderson, S. & Sharpe, M. Robust optimization for intensity modulated radiation therapy treatment planning under uncertainty. Physics in Medicine and Biology 50, 5463-5477 (2005).

39. Thiele, I. & Palsson, B. O. A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat. Protocols 5, 93-121 (2010).

40. Chindelevitch, L., Trigg, J., Regev, A. & Berger, B. An exact arithmetic toolbox for a consistent and reproducible structural analysis of metabolic network models. Nat Commun 5 (2014).

41. Ebrahim, A. et al. Do genome-scale models need exact solvers or clearer standards? Molecular Systems Biology 11 (2015).

42. Chindelevitch, L., Trigg, J., Regev, A. & Berger, B. Reply to "do genome-scale models need exact solvers or clearer standards?". Molecular Systems Biology 11 (2015).

43. Covert, M. W., Schilling, C. H. & Palsson, B. Regulation of gene expression in flux balance models of metabolism. J Theor Biol 213, 73-88 (2001).

44. Ataman, M. & Hatzimanikatis, V. Heading in the right direction: thermodynamics-based network analysis and pathway engineering. Curr Opin Biotechnol 36, 176-182 (2015).

45. Segre, D., Vitkup, D. & Church, G. M. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA 99, 15112-15117 (2002).

46. Stelling, J., Klamt, S., Bettenbrock, K., Schuster, S. & Gilles, E. D. Metabolic network structure determines key aspects of functionality and regulation. Nature 420, 190-193 (2002).

47. Livingstone, K., Holder, A. & Almaas, E. Introduction to systems biology for mathematical programmers. In Lim, G. & Lee, E. (eds.) Optimization in medicine and biology, chap. 11, 311-354 (Auerbach Publications, 2008).

48. Mahadevan, R. et al. Characterization of metabolism in the fe(iii)-reducing organism geobacter sulfurreducens by constraint-based modeling. Appl Environ Microbiol 72, 1558-1568 (2006).

49. Neidhardt, F., Ingraham, J. & Schaechter, M. Physiology of the bacterial cell: a molecular approach (Sinauer Associates, 1990).

Acknowledgements

The authors thank Dr. Eric Reyes for his careful reading of an earlier draft, as well as Jarle Pahr for insights on flux comparisons. E.A. would like to thank The Research Council of Norway (RCN) grant 245160 (ERASysAPP: WineSys) for funding.

Author Contributions

E.A. and A.H. initiated the project. Mathematical analysis by M.M., A.K., E.G., M.S., and A.H. A.H. conducted the numerical experiments. E.A. and A.H. analyzed the results and wrote the paper.

Additional Information

Supplementary information accompanies this paper at doi:10.1038/s41598-017-00170-3 Competing Interests: The authors declare no competing financial interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

l/JJv 0 I This work is licensed under a Creative Commons Attribution 4.0 International License. The images ^■ss^M or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

© The Author(s) 2017