Convex Representation of Metabolic Networks with Michaelis–Menten Kinetics

Polyhedral models of metabolic networks are computationally tractable and can predict some cellular functions. A longstanding challenge is incorporating metabolites without losing tractability. In this paper, we do so using a new second-order cone representation of the Michaelis–Menten kinetics. The resulting model consists of linear stoichiometric constraints alongside second-order cone constraints that couple the reaction fluxes to metabolite concentrations. We formulate several new problems around this model: conic flux balance analysis, which augments flux balance analysis with metabolite concentrations; dynamic conic flux balance analysis; and finding minimal cut sets of networks with both reactions and metabolites. Solving these problems yields information about both fluxes and metabolite concentrations. They are second-order cone or mixed-integer second-order cone programs, which, while not as tractable as their linear counterparts, can nonetheless be solved at practical scales using existing software.


Introduction
The structure of a metabolic network contains useful information about its cellular functions.Two techniques for analyzing this structure are • flux balance analysis (FBA), in which optimization is used to predict reaction fluxes (Orth et al. 2010), and • minimal cut set (MCS) analysis, which attempts to find critical subsets of reactions that, when removed, disable certain functions (Klamt and Gilles 2004).
Standard formulations of FBA and MCS analysis are based on a linear approximation in which only the reaction fluxes are variables.The benefit of this simplification is that FBA is a linear program (LP), which can be reliably solved at large scales, and powerful analytical tools like Farkas' Lemma are available for finding the MCSs of polyhedral systems.Quoting Orth et al. (2010), "FBA has limitations, however.Because it does not use kinetic parameters, it cannot predict metabolite concentrations."For the same reason, MCS analysis cannot explicitly identify critical metabolites.
In this paper, we augment FBA and MCS analysis with Michaelis-Menten kinetics and metabolite concentrations.Using the results of Taylor and Rapaport (2021), we represent the Michaelis-Menten kinetics as a second-order cone (SOC) constraint.This leads to several original problem formulations.
• Conic FBA (CFBA).CFBA predicts both reaction fluxes and metabolite concentrations in steady state.It is a single SOC program (SOCP), which, while not as tractable as LP, can be solved at practical scales (Boyd and Vandenberghe 2004).
We use the dual to derive sensitivities to maximum reaction rates and Michaelis constants.We also formulate dynamic CFBA, in which the SOC representations of the reaction kinetics are used in dynamic FBA (Mahadevan et al. 2002).

• Conic MCS (CMCS).
A CMCS is a cut set through a network of paths from reactions to metabolites, as specified by the stoichiometric matrix, and from metabolites to reactions, as specified by the Michaelis-Menten kinetics.To solve for CMCSs, we follow the strategy of Ballerstein et al. (2012), which uses Farkas' Lemma (Boyd and Vandenberghe 2004) and results from Gleeson and Ryan (1990) and Parker and Ryan (1996) on irreducible infeasible subsystems (IIS) to identify MCSs.We generalize this to CMCSs using the recent results of Kellner et al. (2019) on the IISs of semidefinite systems.We also formulate a linear approximation that, due to the discrete nature of cut sets, produces similar results.
We remark that CFBA and CMCS analysis might not produce better predictions of reaction fluxes than existing methods.The main benefit of CFBA and CMCS analysis is that they incorporate metabolite concentrations and reaction kinetics while retaining much of the tractability of standard linear formulations.We now describe how our contributions relate to the existing literature.Using LP to analyze metabolic networks was first suggested by Watson (1984).Since then FBA has gained wide acceptance (Varma and Palsson 1994;Orth et al. 2010), and is available in open source implementations (Schellenberger et al. 2011;Ebrahim et al. 2013).Dynamic FBA (DFBA) is an extension that incorporates reaction kinetics, which are typically nonlinear, and metabolite concentrations, potentially as well as other transient information such as reprogramming and light intensity.Reference (Mahadevan et al. 2002) first formulated the two main types of DFBA.The 'dynamic optimization approach' is a large nonlinear program, in which the reaction kinetics constrain the fluxes through time.In the 'static optimization approach', one solves a sequence of LPs and integrates the solution between time periods.There have been several refinements such as putting the optimization directly into the ODE simulation (Hanly and Henson 2011), and using lexicographic optimization to improve robustness when the LP has multiple solutions (Gomez et al. 2014).
CFBA is similar to DFBA in that it also captures metabolite concentrations and reaction kinetics, but only when modeled as Michaelis-Menten.They differ in that CFBA is in steady state, and hence does not capture transients.If one approximates the biomass concentration as a constant, or approximates Michaelis-Menten with the Contois function (Contois 1959), one can formulate the dynamic optimization approach of Mahadevan et al. (2002) as an SOCP.We refer to this as Dynamic CFBA.Dynamic CFBA can capture transients and accommodate time-varying parameters.CFBA and dynamic CFBA are not necessarily more accurate than dynamic FBA, but are more tractable in that they are single SOCPs.
CFBA is also related to resource balance analysis (RBA) (Goelzer et al. 2011;Goelzer and Fromion 2011), a more general problem that can predict fluxes, metabolites, macromolecular cellular processes, and proteins.RBA differs from CFBA in that it does not contain the Michaelis-Menten function or any other nonlinearities, and as a result is an LP.In principle, our SOC representation of the Michaelis-Menten function could be incorporated into RBA, leading to an SOCP.
Another relevant literature stream focuses on the analysis of pathways through metabolic networks (see, e.g., Clarke 1980;Schilling et al. 2000).In the linear case, the nonnegativity and stoichiometric constraints form a polyhedral cone, the extreme rays of which correspond to the elementary flux modes of the metabolic network (Schuster and Hilgetag 1994).It is not clear that this perspective extends to our setup because SOC constraints are nonpolyhedral, and the SOC representation of the Michaelis-Menten kinetics is not in fact a cone.The Minkowski-Weyl Theorem states that a cone is finitely generated if and only if it is polyhedral (Schrijver 1998).Therefore, such a system could have an infinite number of elementary modes if it has any at all.If elementary modes do exist, at present we are not aware of any reliable techniques for obtaining them.
For these reasons, we focus on the adjacent problem of identifying MCSs.Reference (Ballerstein et al. 2012) showed that the MCSs of a metabolic network are the elementary modes of a dual network specified by Farkas' Lemma (Boyd and Vandenberghe 2004).We make use of the results of Kellner et al. (2019) to generalize this strategy to networks with metabolite concentrations coupled to the reaction fluxes through Michaelis-Menten kinetics.By then using a linear approximation of the Michaelis-Menten kinetics, we recover the use of tools for polyhedral systems, which we find produce more reliable results.
The paper is organized as follows.Section 2 reviews metabolic network modeling and the SOC representation of the Michaelis-Menten kinetics.Section 3 presents FBA, CFBA, and dynamic CFBA.Section 4 presents the MCS analysis of Ballerstein et al. (2012) and our extensions to systems with Michaelis-Menten kinetics.In Sect.5, we apply CFBA, dynamic CFBA, and CMCS analysis to a model of Escherichia coli.

Metabolic Networks
The reactions are indexed by the set N , where n = |N |.Let R ⊆ N and I ⊆ N be the sets of reversible and irreversible reactions, where R ∪ I = N and R ∩ I = ∅.There are a total of m metabolites.
Let z ∈ R m + be a vector of metabolite concentrations and v ∈ R n a vector of fluxes due to the reactions.S ∈ R m ×n is the stoichiometric matrix.The dynamics of the metabolites are given by where z b is the element of z corresponding to biomass (Mahadevan et al. 2002;Hanly and Henson 2011;Höffner et al. 2013).While the biomass is not a metabolite, its evolution can be written as z b times a linear combination of the reaction fluxes, and hence can be represented as a row of (1).
In quasi-steady state, ż = 0. We then have z b Sv = 0, which, because z b is a nonzero scalar, implies Sv = 0.The fluxes are subject to the bounds v ≤ v ≤ v. Typically, v i = 0 for i ∈ I and −∞ otherwise.When v = 0 and v = ∞, the system is a polyhedral cone (Schuster and Hilgetag 1994).Quasi-steady state models generally do not explicitly model the metabolites.Here we will include a subset of them, M, where m = |M|.We will henceforth let z ∈ R m + .Some of the fluxes are also bounded by a nonlinear function, usually the Michaelis-Menten kinetics (Michaelis and Menten 1913).We denote this set by Q ⊆ N , where q = |Q|.We expect that q ≥ m, because otherwise there are metabolites that are not substrates in any reaction, and thus have no coupling to the rest of the model.To simplify notation, we order N so that its first q elements are Q.Let V max ∈ R q + and K m ∈ R q + be vectors of maximum reaction rates and Michaelis constants.We can write this bound for each reaction i ∈ Q as . (2) Here σ (i) identifies the index of the metabolite concentration appearing in reaction i.In general σ is not invertible because the same metabolite can appear in multiple reactions.Together with z ≥ 0, the inequality (2) defines a convex set because the right-hand side is concave for z σ (i) ≥ −K m i , as can be shown, e.g., by taking the second derivative.Note that if i is a reversible reaction, (2) may be applied in the reverse direction by putting a minus sign in front of v i .We note that there are other ways to model reversible reactions, e.g., by making the product metabolite the substrate of a different kinetics (Ndiaye and Gouzé 2013).
When a reaction has multiple reactants, its flux may be limited by more complicated functions such as the product of several Michaelis-Menten kinetics (cf. Table IV in Chassagnole et al. (2002)).We can straightforwardly generalize our notation to this case.Suppose that reaction i has p i reactants.For k = 1, ..., p i , let σ k (i) be the index of the k th metabolite in the reaction.Then we can write the upper bound on reaction i ∈ Q as Note that the product, is the maximum reaction rate of reaction i, V max i .The intermediary variables, v ik , k = 1, ..., p i , will be convenient for representing (3) in SOC form.We remark that this is not a general model, and that there are reactions with multiple reactants that are not well-described by (3).

Convex Representation of the Michaelis-Menten Kinetics
As shown by Taylor and Rapaport (2021), we can represent the inequality (2) as the SOC constraint ⎡ ⎣ for i ∈ Q,1 One can confirm equivalence by squaring both sides and simplifying.We refer the reader to Appendix A for a brief introduction to SOCP.As will be seen later, the advantages of the SOC formulation include amenability to powerful, specialized solvers, and more refined analytical tools such as conic duality.Suppose now there are multiple reactants.As with (2), we can write (3a) as an SOC constraint in the form of (4).Unfortunately, (3b) is nonconvex and has no SOC representation, and so does not fit the problems we later formulate.As such, one potential approximation is which has an SOC representation if (Alizadeh and Goldfarb 2003).When π ik = 1/ p i for k = 1, ..., p i , it is the geometric mean.Equation ( 5) is ad hoc in that it cannot in general be derived from first principles.However, we believe it could be a useful approximation because there is considerable room to tune the fit via the exponents and Michaelis-Menten parameters; and, similar ad hoc expressions have been used in the past when Michaelis-Menten alone was not sufficiently descriptive (cf. Table IV in Chassagnole et al. (2002)).The exponents could be fit using nonlinear least squares.
In the simple case where p i = 2 and π i1 = π i2 = 1/2, (5) is a hyperbolic constraint with SOC representation We emphasize that as ( 5) and ( 6) have not been used in prior studies, they require validation via simulation or experiments; this is a topic of future work.
Sometimes LP is more practical to work with than SOCP, e.g., when there are numerous discrete constraints or commercial solvers are too expensive.Fortunately an SOCP can always be approximated to arbitrary accuracy with an LP, albeit a potentially large one.Ben-Tal and Nemirovski (2001) provide a constructive procedure for approximating a generic SOC constraint with a family of linear constraints.Alternatively, the right hand side of (2) can be straightforwardly approximated by a family of line segments.

Flux Balance Analysis
In this section, we first review the standard formulation of FBA, and then formulate CFBA and Dynamic CFBA.We also use the duals to derive several sensitivities.

Linear FBA
The below LP is a basic FBA routine.
Here c ∈ R n + selects and/or weights the fluxes for maximization.Let F denote the optimal objective value.
Let λ ∈ R m be the vector of dual multipliers of (7b), and let δ L ∈ R n + and δ U ∈ R n + be the vectors of dual multipliers associated with the upper and lower bounds in (7c).
The dual of ( 7) is min The dual variables, or shadow prices, can be interpreted as the sensitivities of the objective in (7) to changes in the constraints.More precisely, for each i ∈ N , λ is similarly interpretable as the sensitivity of the objective to perturbations to (7b).This, for example, can be used to determine which reactions are most influential (Schilling et al. 2000;Reznik et al. 2013).
We an also obtain insight by observing that the objectives of ( 7) and ( 8) must be equal due to strong duality: This breaks down the optimal objective into contributions from each constraint.In the common case where we maximize a single reaction, v s , and v = 0, this simplifies to v s = δ U v, a convex combination of the upper flux limits.

Conic FBA
We now generalize the FBA (7) to include metabolites, which are coupled to the reactions by the Michaelis-Menten kinetics.Let z ∈ R m + and z ∈ R m + be vectors of upper and lower bounds on the metabolite concentrations.Let d ∈ R m − be a vector weighting the metabolites in the objective.Consider the below optimization.
We make the following comments.
• This is an SOCP if we write (9e) in the form of (4).
• Within the subset of reactions limited by Michaelis-Menten kinetics, Q, it may be that not all are important, and therefore that not all metabolites are limiting.We can identify which metabolites are limiting through sensitivity analysis, as described later in this section.It is these influential metabolite concentrations that we expect CFBA to predict well.• We could include reaction limits in the form of ( 5) that depend on multiple Michaelis-Menten kinetics, and retain the SOC structure.We have not done so as to retain simplicity in (9) and its dual, and because the numerical examples in Sect. 5 only have upper bounds with single Michaelis-Menten reactions.
• (9e) is only meaningful if the corresponding metabolite concentrations are minimized in the objective, i.e., d < 0. Otherwise, it is trivially optimal to fix the metabolites at their maximum concentrations, z, in which case (9e) can be represented as an upper bound on v in (9c).• The objective is interpretable as maximizing the fluxes specified by c while minimizing the metabolite concentration.We discuss an alternative formulation below in Sect.3.2.1.
We now derive the dual.For matrices A and B of the same size, let A • B denote the element-wise product.We let A i: and A : j denote the i th row and j th column of A.
Let λ ∈ R m be the vector of dual multipliers of (9b), and zero otherwise.The interpretation of M is complementary to that of S. S encodes a network representing how the reactions influence the evolution of the metabolite concentrations.Similarly, M encodes a network representing which metabolites appear in each reaction.
The dual of ( 9), which is also an SOCP, is below.
Strong duality holds if a constraint qualification is satisfied, e.g., Slater's condition (Boyd and Vandenberghe 2004).In this case λ, δ L , λ, δ U , γ L , and γ U all have the usual LP sensitivity interpretations and complementary slackness with their respective constraints.
We can similarly interpret and φ; see, e.g., Section 5.9.3 of Boyd and Vandenberghe (2004).φ i is the sensitivity of the optimal objective value, F, to perturbations to the right hand side of (4).:i is a vector of sensitivities of F to perturbations to each element in the left hand side.We can use the chain rule to derive sensitivities for the parameters of the Michaelis-Menten function.For i ∈ Q we have If the sensitivity to K m i is high, then the corresponding metabolite concentration, z σ (i) , is limiting in the sense that a small change will substantially change the reaction flux, v i , and the objective value.On the other hand, if the sensitivity to V max i is high and K m i low, then the optimal solution is on the flatter, rightward part of the Michaelis-Menten function, and so the corresponding metabolite concentration has little influence on the solution.In this case, the Michaelis-Menten constraint could be replaced with a simple upper limit on v i .
As in the previous section, the objectives of ( 9) and ( 13) match if strong duality holds.If we assume that v = 0, z = 0, and z = ∞, the equality simplifies to

Alternative Objectives
A shortcoming of ( 9) is that the objective-flux rates minus metabolite concentrations-mixes units.This is problematic because the corresponding value does not have a clear interpretation, e.g., the concentration of a metabolite of interest; and, the two terms might differ numerically by orders of magnitude.We need not directly compare fluxes and metabolite concentrations if we instead maximize d z alone subject to a constraint on v, e.g., v s ≥ 1, where s is the index of a reaction of interest.This is interpretable as the minimum metabolite concentration necessary to carry out a certain function.Such a constraint could be incorporated into (9c).
An objective consisting only of metabolite concentrations might be more physically interpretable.For instance, minimizing internal metabolite concentrations was used as a model of cellular function in Schuster and Heinrich (1991).A second advantage is that the dual sensitivities in (11) also have physical meanings.For these reasons, we maximize −1 z in the example in Sect.5.1, where 1 is the appropriately sized vector of all ones.
We could also maximize c v alone and add constraints on z.This would be nontrivial, e.g., with coupled polyhedral constraints.However, box constraints like (9d) would simply lead to z binding at its upper bounds.
Note that while the dual is slightly different for these alternative formulations, the expressions for the sensitivities to V max and K m in (11) are unchanged.

Dynamic FBA
Dynamic FBA extends conventional FBA in two ways: the incorporation of metabolite concentrations, and changes over time not captured in a steady state model.Metabolite concentrations are present in CFBA, but not transients.We thus formulate dynamic CFBA for the purpose of capturing transient phenomena.
There are two main approaches to dynamic FBA, as described by Mahadevan et al. (2002).
• In the dynamic optimization approach, the full trajectory is optimized.The concentrations across time periods are coupled by a finite difference approximation of the derivative; note that we could use a more accurate approximation, e.g., a higher order Runge-Kutta scheme (Betts 1998).The optimization is nonlinear due to the Michaelis-Menten kinetics and bilinearities between the fluxes and biomass variables.• In the static optimization approach, a conventional FBA is solved in each time period.The resulting flux vector is used to propagate the concentrations to the next time period.
The first is harder to solve because it is larger and nonlinear.Here, we optimize the full trajectory as in the dynamic optimization approach.We render the problem more computationally tractable by representing the Michaelis-Menten kinetics in SOC form, (4).
A difficulty not present in steady state FBA is the biomass, z b , which we recall is an element of the vector z.Because we are not setting the derivatives in (1) to zero, we cannot divide it out, and the product of the Michaelis-Menten function and z b does not have an SOC representation.To maintain consistent notation, we introduce the variable ν ∈ R n , which takes the place of the product z b v.We can recover convexity in two ways.
1.If the biomass evolution is predictable, we can approximate it as a fixed parameter in each time period, zb (t).Then the kinetics can be represented as an SOC constraint in the form of (4).2. Instead of Michaelis-Menten, we can model the reaction limits with the Contois function (Contois 1959): The Contois function is commonly used to model biochemical processes (Bastin and Dochain 1990), and differs from Michaelis-Menten only in that the denominator depends on z b .Given its similarity to the Michaelis-Menten kinetics, (12) could also be an acceptable approximation.As shown by Taylor and Rapaport (2021), we can represent the inequality (12) as the SOC constraint Which approximation is more appropriate depends on the context.For example, if the biomass does not vary much, holding it constant in Michaelis-Menten is a clear choice.
If the evolution of the biomass does depend strongly on the other metabolites, then the Contois function might be a better choice.Alternatively, we could fit a generic SOC constraint to reaction data, which would retain compatibility with SOCP and potentially be more accurate than either of the above approximations (Tan et al. 2022).
There are time periods t = 0, ..., τ , each of length .Let ν(t) and z(t) denote the fluxes and metabolite concentrations in period t.We thus obtain the below SOCP for dynamic CFBA.max (13b) is the initial condition.(13c) is the Euler approximation of (1), and couples the variables across time periods.If we approximate the biomass as a fixed parameter, we also have for t = 1, ..., τ : If we model the reaction limits with the Contois function, (12), then instead of (13f) and (13g), we have: (13d) and (13f) (or (13h)) are bounds on the reactions and metabolites.(13e) limits the rate of change of the reactions.In (13g) (or (13i)), ν i (t) is less than either the Michaelis-Menten kinetics with fixed biomass or the Contois function, depending on which of the above approximations is used.
In dynamic FBA, the evolution of the metabolites generally depends on more than just the stoichiometry, as in (13c).The same should also be the case for dynamic CFBA.For example, there could be inflows and outflows with endogenous or exogenous metabolite concentrations, or biomass death (Gomez et al. 2014).

Minimal Cut Sets
In (11) in Sect.3.2, we used dual variables to compute sensitivities to kinetic parameters.This was one of several potential applications of the CFBA dual system.We now explore another in which we use Farkas' Lemma to identify CMCSs.
An MCS is the smallest set of reactions which, if constrained to be zero, disables some function of interest.In this regard, the reactions in an MCS are the lynchpins of the system.Ballerstein et al. (2012) showed that the MCSs of a metabolic network correspond to the elementary modes of a dual network, which is specified by Farkas' Lemma (Schrijver 1998).This was an application of a result from Gleeson and Ryan (1990) and Parker and Ryan (1996), which states that there is a one-to-one mapping between the IISs (irreducible infeasible subsystems) of a linear system and the vertices of its dual polyhedron.One can therefore identify the MCSs of a metabolic network by solving for the vertices of the dual polyhedron, which can be done via mixed-integer LP (MILP).
In Kellner et al. (2019), a weaker version of the result of Gleeson and Ryan (1990) and Parker and Ryan (1996) was extended to semidefinite systems, which generalize SOC systems.We make use of this in Sect.4.1 to extend the results of Ballerstein et al. (2012) to networks with metabolites linked by Michaelis-Menten kinetics.Whereas an MCS contains only reactions, a CMCS contains the smallest set of reactions and/or metabolites that the system cannot function without.We then formulate a linear approximation that produces similar results in Sect.4.2.

Conic Minimal Cut Sets
We now describe the setup, starting with the linear part as given by Ballerstein et al. (2012).Let T v ⊆ N denote the polyhedral set of target reactions, parametrized by the matrix T and vector v * .The target set, which does not contain the origin, encodes some function of interest, which the removal of a cut set disables.The constraint forces the reactions in T v to be active.The below two constraints encode the steady state operation of the metabolic network: 14b) and (14c).It is an MCS if it contains no smaller cut sets, i.e., cut sets with fewer elements, for T v .
We now incorporate metabolite concentrations and reaction kinetics.We assume there is also a target set of metabolites, T z , each of which is constrained to be in some concentration range by where W and z * are an appropriately dimensioned matrix and vector.T z similarly does not contain the origin.The below two constraints prohibit negative concentrations and relate the concentrations to the reactions: Similar to (14d), we make this portion of the system infeasible by adding the constraint which, together with (14f), implies z = 0. We note that infeasibility is not always precise enough when moving from linear to nonlinear systems because in the latter case, the set of solutions might not be closed.As in Kellner et al. (2019), we say that a semidefinite system is weakly feasible if any positive perturbation to its eigenvalues makes it feasible, and weakly infeasible if it is not weakly feasible.Feasibility implies weak feasibility, and weak infeasibility implies infeasibility.( 14) is weakly infeasible, and hence also infeasible, because there are positive perturbations that do not make it feasible.We henceforth take the definition of an IIS to refer to weak infeasibility.
We now extend the definition of an MCS.14c), (14f), and (14g) if the additional constraints v i = 0 for i ∈ C v and z i = 0 for i ∈ C z imply that v i = 0 for i ∈ T v and z i = 0 for i ∈ T z .It is a CMCS if it contains no smaller cut sets.

Lemma 1 Each CMCS
The proof is similar to that of Lemma 1 of Ballerstein et al. (2012).
The definition specifies an infeasible system consisting of (14a)-( 14c), ( 14e)-(14g), and (15) Denote this system .is a subsystem of ( 14) because ( 15) is a subsystem of v = 0 and z = 0.If is not irreducible, it must contain an IIS.Because the cut set is minimal, the removal of any element from C v or C z makes feasible.Therefore, any IIS of must contain (15).Because C = {C v , C z } is distinct to the CMCS, any IIS of is also distinct to the CMCS.
As discussed for MCSs by Ballerstein et al. (2012), Lemma 1 has the following two implications: while each IIS corresponds to at most one CMCS, a CMCS can correspond to multiple IISs; and there may be IISs that do not correspond to a CMCS.We now seek to relate the IISs of ( 14) to some dual system, which we recall was specified by Farkas' Lemma in the linear case by Ballerstein et al. (2012).Farkas' Lemma does not apply to ( 14) due to the SOC constraint, (14g).There are several extensions to semidefinite systems, e.g., (Klep and Schweighofer 2013), which (Kellner et al. 2019) employs to generalize the results of Gleeson and Ryan (1990) and Parker and Ryan (1996) to semidefinite systems.The dual system in this case is also referred to as the alternative spectrahedron.Because any SOC constraint can be written as a semidefinite constraint, the results of Kellner et al. (2019) specialize to SOC systems like ( 14) without modification.
The dual system of ( 14) is λ is the dual variable of (14b); δ of v = 0; γ of z = 0; ρ v of (14a); ρ z of (14e); and ( i , φ i ), i ∈ Q, of the SOC form of (14g).Theorem 3.2 of Kellner et al. (2019) states that if ( 14) is weakly infeasible, then each of its IISs corresponds to an extremal point of ( 16).More precisely, each IIS determines the nonzero entries of some extremal point of ( 16).We can thus identify CMCSs by finding the extremal points of ( 16).
Unfortunately, ( 16) may have extremal points that do not correspond to an IIS of ( 14).Theorem 4.1 of Kellner et al. (2019) provides conditions under which the alternative spectrahedron has a single solution; however, it does not appear to apply in general to (16).
Lemma 3.1 of Kellner et al. (2019) states that the indices of a minimal cardinality solution of ( 16) correspond to an IIS of ( 14).In our case, similar to Ballerstein et al. (2012), we seek solutions in which the vectors δ and γ have minimal cardinality.
In the linear case, there are several different MILPs for solving for IISs.Likewise, one could formulate multiple mixed-integer SOCPs (MISOCP) for finding the IISs of ( 14); e.g., (9) of Kellner et al. (2019) can be expressed as an MISOCP when specialized to our problem.The MISOCP below differs from Kellner et al. (2019) in that we only seek minimal cardinality in δ and γ , and not the other variables.The portion of the problem corresponding to reaction fluxes is based on the MILP ( 14) of Klamt et al. (2020).
β δ and β γ are vectors of binary variables.( 17b)-( 17f) are disjunctive constraints that, for large enough ϒ, either force their left hand sides to be zero or have no effect.The left side of ( 17f) is the 2-norm of the SOC variable ( :i , φ i ) (Alizadeh and Goldfarb 2003).We include this constraint because if β γ i = 0, the dual variables of ( 14g) and (14h) should be zero.Given a solution, the IIS is specified by the entries of β δ and β γ that are equal to one.
There are a number of further refinements one can make to (17).For example, by weighting the terms in the objective, one can promote cuts with reactions or with metabolites.To exclude either the reactions or metabolites in the target set, or a cut that has already been found, say βδ , βγ , we can add the constraint

Linear Approximation
There are severals reasons why solving (17) might not be the best way to identify CMCSs.First, as discussed, the underlying theoretical results are weaker than in the linear case-an extremal point of the alternative spectrahedron might not correspond to an IIS (Kellner et al. 2019).Second, MISOCP is less tractable than MILP.And third, given the discrete nature of cut sets, it is not clear that the nonlinearity could not be replaced with something simpler.
The following is a standard linear approximation for when the metabolite concentration is smaller than the Michaelis constant: where is in a cut set, i.e., we set it to zero, then v i ≤ 0 under either (14g) or (18).Definition 1 and Lemma 1 also hold in both cases.Note that while this definition of η has precedent, any positive value would serve similarly.
Here (20d) serves the same role as (17f), but is linear because ( 18) is a linear inequality.Note that the system Sv = 0, v ≥ 0, z ≥ 0, and ( 18) is a polyhedral cone, of which we could therefore analyze the extreme rays.This is a topic of future work.

Escherichia coli Example
We now apply the tools we have developed to an example based on the model e_coli_core in the BiGG database (King et al. 2016), which corresponds to Escherichia coli str.K-12 substr.MG1655.The model has 95 reactions and 72 metabolites, of which we explicitly model the twelve that appear in Table 1.
We take all parameters for the FBA routine (7) from the COBRA Toolbox (Schellenberger et al. 2011;Ebrahim et al. 2013).We augment the model with Michaelis-Menten kinetics for the reactions listed in Table 1.The values of K m and V max for each reaction were taken from Table 1 from Meadows et al. (2010).
All optimizations were carried out in Python using CVXPy (Diamond and Boyd 2016) and the solver Gurobi (Gurobi 2021).All figures were made with Matplotlib (Hunter 2007).

CFBA
We first apply CFBA.We solve the alternative formulation described in Sect.3.2.1, which is identical to (9) except in two respects.
• The objective is −1 z, which corresponds to minimizing the sum of the concentrations of the metabolites in Table 1.• We add the constraint v b ≥ .v b is the flux of the biomass reaction, which is BIOMASS_Ecoli_core_w_GAM in the BiGG database.

Table 1
Michaelis-Menten parameters.The of K m and V max are mM and mmol/g dw/h, and the values are from Meadows et al. (2010).The metabolites and reactions are listed along with their ID in the BiGG database (King et al. 2016) # We are thus finding the minimum metabolite concentrations necessary to keep the biomass flux above .
The SOCP for CFBA took 0.08 s to solve.For comparison, the LP for FBA took 0.007 s, roughly an order of magnitude less.
Figure 1 shows the sensitivities of the optimal objective to V max and K m for three of the reactions in Table 1 for = 0.1.Given primal and dual CFBA solutions, the sensitivities are computed via (11).Each represents the change in total metabolite usage resulting from a small change in a Michaelis-Menten parameter, given that the biomass reaction flux cannot go below .Of the reactions that are not shown, 1 and 2 have sensitivities on the order of 10 −4 , and the rest 10 −8 .
We can see that the objective is most sensitive to the kinetics of Reaction 3 (ACONTa), which depends on the concentration of citrate.This reaction produces H 2 0 for a number of other reactions.Via inspection of the network structure and optimal reaction fluxes, we can see that it also directly enables the reactions aconitase (half-reaction B, Isocitrate hydrolyase) and then isocitrate dehydrogenase (NADP).
The objective is also sensitive to K m for Reactions 12 (GLCpts) and 17 (O2t), which depend on glucose and oxygen, indicating that an increase in K m for either reaction will significantly increase the amount of metabolite needed to keep the biomass reaction flux at .
We may interpret this as follows.High sensitivity to K m indicates that at the optimal solution, we are near zero, on the steeply increasing part of the Michaelis-Menten function-where a slight increase in concentration significantly increases the reaction rate.On the other hand, high sensitivity to V max indicates that we are on the flatter, right side of the Michaelis-Menten function-where changing the metabolite concentration does not significantly affect the reaction rate.Note that conventional FBA does not provide this information because it identifies which reactions are important, but not how they depend on the metabolite concentrations.
In this example, CFBA reveals that Reactions 3, 12, and 17 are successively further to the left on the steeply increasing part of the Michaelis-Menten function.This means that a slight increase in oxygen would significantly increase the maximum rate of Reaction 17, whereas a slight increase in citrate and glucose would moderately increase the maximum rates of Reactions 3 and 12, respectively.
Figure 2 shows the concentrations of citrate, glucose, and oxygen as , the required biomass flux, increases from 0.05 to 0.5.Because these are limiting metabolites, their concentrations increase with , and do so at a greater than linear rate.The concentrations of citrate and glucose increase rapidly with because, at the optimal solution, they are further to the right on flatter part of the Michaelis-Menten function.The concentration of oxygen increases least rapidly because, as described above, a slight increase dramatically increases the reaction rate.
These concentrations are consistent with reported ranges.When = 0.5, the predicted concentration of citrate is 1.17 mM, within the range of 1.1 to 3.5 mM reported in Supplementary Table 3 of Bennett et al. (2009).Escherichia coli's metabolism can function over a wide range of oxygen concentrations.When = 0.5, the predicted oxygen concentration is 0.01 mM, which, for example, falls well within the range depicted in Figure 4A of Henkel et al. (2014).

Dynamic CFBA
We now test dynamic CFBA by solving (13).Our secondary goal in this section is to understand the scalability of CFBA and its dynamic extension, which we do by varying the time horizon in (13), τ .
The objective, (13a), is to maximize the biomass concentration in the last period, z b (τ ), as in equation (6b) in Case 2 of Mahadevan et al. (2002).The initial biomass concentration is z b0 = 0.001, and the remaining elements of z 0 are ones.The biomass concentration evolves as where v b (t), as described in Sect.5.1, is the biomass reaction flux.
We solved (13) for time horizons ranging from τ = 50 to τ = 500.In each instance the time step was = 0.001 hours.Figure 3 shows the time taken by the solver as a function of τ .When τ = 50, there are 8,473 variables, and when τ = 500, there are 84,073 variables.The trend is increasing with problem size, though some smaller instances take longer than larger ones, presumably due to specific problem structure and solver behavior.The longest time taken by the solver was roughly six minutes.This confirms that, by virtue of being SOCPs, CFBA and dynamic CFBA are highly tractable problems.
The upper plot in Fig. 4 shows the biomass, citrate, glucose, and glutamine concentrations through time for the case when τ = 500.Recall that in Sect.5.1, citrate The biomass concentration increases from near zero, first exponentially, and then more gradually-this is a standard behavior.The glutaminase reaction first increases as the biomass concentration increases, and then decreases as the concentration of glutamine drops.The aconitase reaction also increases as biomass concentration increases and then levels out.
These observations point to several potential refinements of dynamic CFBA.First, to induce more complicated behaviors like loss or death of biomass, one must include more complicated exogenous conditions; e.g., in Meadows et al. (2010) there are five different metabolic phases, some of which have distinct exogenous inputs.
A second potential shortcoming is that dynamic CFBA assumes too much foresight, in that the entire trajectory of 500 periods is optimized at once, and with full knowledge of future exogenous inputs.Full foresight is consistent with the dynamic optimization approach, as described in Sect.3.3.However, it is the static optimization approach that finds more usage today, wherein optimization only occurs within individual time periods.A potential remedy is to implement dynamic CFBA with a shorter horizong in receding horizon fashion (Mattingley et al. 2011).This would both limit foresight and enhance scalability.

CMCS Analysis
We compute CMCSs of our example by solving (17), an MISOCP, and ( 20), an MILP.In addition to the setup at the start of this section, we must also specify target sets, T v and T z , in the form of (14a) and (14e).Table 2 lists reaction 'knockouts' and the corresponding CMCSs found by ( 17) and ( 20), specified by their IDs in the BiGG database (King et al. 2016).The target set in each case consists of constraining the flux to be greater than one.Note that we only list the first CMCSs found by the solver, that one target set can have many CMCSs, and that a given CMCS can be a solution for both (17) and (20).
The CMCSs in Table 2 consist of reactions, metabolites, or combinations of both.For example, to disable the glutamine synthetase reaction (GLNS), we can eliminate the substrate L-Glutamine (gln__L_c) and disable the glutamate dehydrogenase (NADP) (GLUDy) and ammonia reversible transport (NH4t) reactions.On the other hand, conventional MCS analysis can only identify critical reactions, not metabolites.
The prevalence of metabolites in the CMCSs depends on which reactions are limited by Michaelis-Menten kinetics.For instance, fum_c and succ_c are both CMCSs because they limit FUM and SUCDi, among other reactions.Note that adding more metabolites and Michaelis-Menten constraints will increase the number of CMCSs, but does not invalidate those that exist for a smaller number of metabolites.
The computation time in all but the final case was under one second.In the last case, (20) took six seconds to solve, and (17) was still not solved after an hour.This highlights the fact that MILP and MISOCP are NP-hard, and similar instances of a problem can take very different times to solve.
Problem (20) often has the same or similar solutions to (17) and, due to the higher tractability and maturity of MILP, is easier to solve in some cases.For these reasons, (20) appears to be more practical than (17).

Conclusion
By representing the Michaelis-Menten kinetics as second-order cone constraint, we can add metabolite concentrations to standard models of metabolic networks without losing much tractability.This has enabled us to formulate several new tools: conic flux balance analysis, dynamic conic flux balance analysis, and conic minimal cut set analysis.In our numerical examples, we demonstrated that each of these new problems is tractable and can provide insight into both reaction fluxes and metabolite concentrations.
There are several directions for future work.We that there are numerous potential applications to the many different organisms there are.Such studies could both provide new insights into metabolic networks and further clarify when these tools are appropriate.A starting point for this is applying them to larger metabolic network models.This entails augmenting more existing models with Michaelis-Menten parameter data.There is also room for methodological advancements.For example, a receding horizon implementation (Mattingley et al. 2011) could make dynamic conic flux balance analysis both more realistic and tractable.The SOC representation of Michaelis-Menten could be incorporated into flux variability analysis (Gudmundsson and Thiele 2010) so as to find near optimal ranges of both fluxes and metabolites.There is certainly more to understand about the basic geometry of our setup, which, though convex, is not amenable to many of the techniques used to analyze polyhedral models.

A Second-order Cone Programming
We here provide cursory background on SOCP, and refer the reader to Alizadeh and Goldfarb (2003) and Boyd and Vandenberghe (2004) for in-depth coverage.Let A ∈ R m×n , b ∈ R m , c ∈ R n , and d ∈ R. A standard form SOC constraint on the variable x ∈ R n is written where the left-hand side is the two-norm.If A = 0, this reduces to a linear constraint.
A commonly occurring constraint with SOC form is the hyperbolic constraint x 2 1 ≤ x 2 x 3 .It can be written 2x 1 x 2 − x 3 ≤ x 2 + x 3 .
We must also require x 2 ≥ 0 and x 3 ≥ 0.
An optimization problem with a linear objective and p SOC constraints is an SOCP, and can in general be written

Fig. 1
Fig. 1 Sensitivities of the CFBA objective, the sum of the metabolite concentrations, to V max and K m for ACONTa, GLCpts, and O2t.The lower bound on the biomass reaction flux is = 0.1

Fig. 2
Fig. 2 Optimal concentrations of citrate, and oxygen as a function of the lower bound on the biomass reaction flux, , in CFBA

Fig. 3
Fig. 3 Computation time of CFBA as a function of the number of periods, τ