Next Article in Journal
Nonlinear Dynamics and Performance Analysis of a Buck Converter with Hysteresis Control
Next Article in Special Issue
Exact Boolean Abstraction of Linear Equation Systems
Previous Article in Journal
Forecasting Multivariate Chaotic Processes with Precedent Analysis
Previous Article in Special Issue
A Language for Modeling and Optimizing Experimental Biological Protocols
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Metabolic Pathway Analysis in the Presence of Biological Constraints

Université Paris-Saclay, CNRS, ENS Paris-Saclay, Inria, Laboratoire Méthodes Formelles, 91190 Gif-sur-Yvette, France
Computation 2021, 9(10), 111; https://doi.org/10.3390/computation9100111
Submission received: 29 August 2021 / Revised: 13 October 2021 / Accepted: 13 October 2021 / Published: 19 October 2021
(This article belongs to the Special Issue Formal Method for Biological Systems Modelling)

Abstract

:
Metabolic pathway analysis is a key method to study a metabolism in its steady state, and the concept of elementary fluxes (EFs) plays a major role in the analysis of a network in terms of non-decomposable pathways. The supports of the EFs contain in particular those of the elementary flux modes (EFMs), which are the support-minimal pathways, and EFs coincide with EFMs when the only flux constraints are given by the irreversibility of certain reactions. Practical use of both EFMs and EFs has been hampered by the combinatorial explosion of their number in large, genome-scale systems. The EFs give the possible pathways in a steady state but the real pathways are limited by biological constraints, such as thermodynamic or, more generally, kinetic constraints and regulatory constraints from the genetic network. We provide results on the mathematical structure and geometrical characterization of the solution space in the presence of such biological constraints (which is no longer a convex polyhedral cone or a convex polyhedron) and revisit the concept of EFMs and EFs in this framework. We show that most of the results depend only on very general properties of compatibility of constraints with vector signs: either sign-invariance, satisfied by regulatory constraints, or sign-monotonicity (a stronger property), satisfied by thermodynamic and kinetic constraints. We show in particular that the solution space for sign-monotone constraints is a union of particular faces of the original polyhedral cone or polyhedron and that EFs still coincide with EFMs and are just those of the original EFs that satisfy the constraint, and we show how to integrate their computation efficiently in the double description method, the most widely used method in the tools dedicated to EFs computation. We show that, for sign-invariant constraints, the situation is more complex: the solution space is a disjoint union of particular semi-open faces (i.e., without some of their own faces of lesser dimension) of the original polyhedral cone or polyhedron and, if EFs are still those of the original EFs that satisfy the constraint, their computation cannot be incrementally integrated into the double description method, and the result is not true for EFMs, that are in general strictly more numerous than those of the original EFMs that satisfy the constraint.

1. Introduction

1.1. Metabolic Networks

In order to ensure this paper is self-contained and has no prerequisite to be read, we summarize in this introduction the state-of-the-art related to the subject and fix the notations adopted throughout the paper. The results quoted being known, are thus given without proof and the reader is invited to refer to the works in [1,2,3,4,5,6,7], in addition to the references in the text, for the proofs or surveys.
A metabolic network is made up of a set of r biochemical enzymatic reactions. Each reaction consumes certain metabolites (called substrates of the reaction) and produces other metabolites (called products of the reaction). Each metabolite is assigned a coefficient in the reaction, its stoichiometric coefficient (counted negatively for substrates and positively for products). We distinguish internal from external metabolites w.r.t. the system under study (e.g., a bacterium, an eukaryotic cell), the reactions involving both internal and external metabolites being transfer reactions. If m is the number of internal metabolites ( r > m ), the network is thus given by its stoichiometric matrix S R m × r , where coefficient S j i is the stoichiometric coefficient of internal metabolite j in reaction i (positive if j is a product of reaction i and negative if it is a substrate). A state of the network at a given time t is given by the net rates (or fluxes) in each of its reactions at t, i.e., by a flux vector (or rate vector, or flux distribution) v ( t ) R r . Denoting by M ̲ ( t ) R + * m the vector of the concentrations of internal metabolites at t, the time evolution of the network is thus given by
d M ̲ ( t ) d t = S v .

1.2. Steady-State Behavior and Flux Subspace

We are interested in the steady-state behavior of the network. The steady-state assumption means that the concentrations of internal metabolites remain constant along time (no accumulation or reduction of internal metabolites inside the system, an approximation which is valid for short time periods, e.g., a few minutes) and leads thus to the fundamental equation
S v = 0 .
With only this assumption, the solution space S o l S , i.e., the space of all admissible flux vectors v , is thus the linear subspace F S of R r given by the kernel, or nullspace [8], of S :
S o l S = F S = { v R r S v = 0 } .
The dimension of the flux subspace is given by d i m ( F S ) = r r a n k ( S ) r m . Often, possibly after a preprocessing to eliminate its linearly dependent rows, S is assumed to be of full rank and then d i m ( F S ) = r m .
The support of a vector x R r is defined by
s u p p ( x ) = { i x i 0 } .
The support of a flux vector v has thus an important biological signification as it represents the reactions involved in the subnetwork (that we shall call pathway) defined by v (i.e., those through which the flux given by v is not null).

1.3. Irreversible Reactions, Flux Cones and Polyhedral Cones

In addition to the homogeneous equality constraints provided by the steady-state assumption, the flux vectors, in general, must also satisfy homogeneous inequality constraints corresponding to reactions known as irreversible, whose set will be noted Irr :
v i 0 for i Irr .
This means that fluxes in irreversible reactions are constrained to be non-negative (in reversible reactions, fluxes may be either positive, or negative or null and the direction fixed as positive is arbitrary, the role between substrates and products being able to switch). If r I = | Irr | , with 0 r I r , is the number of irreversible reactions, the solution space S o l S , Irr is the intersection of the linear subspace F S with r I non-negative half-spaces, it is thus a particular case of a convex polyhedral cone, called s-cone (subspace cone or special cone) or flux cone, noted F C :
S o l S , Irr = F C = { v R r S v = 0 , v i 0 for i Irr } .
A (general) convex polyhedral cone is defined implicitly (or by intension) by finitely many homogeneous linear inequalities:
C = { x R r A x 0 }
with a suitable matrix A R n × r (called a representation matrix of C) and is thus the intersection of n half-spaces whose frontiers contain the origin. The dimension of C, noted d i m ( C ) , is defined as the dimension of its affine span. A flux cone corresponds thus to a particular matrix A R ( 2 m + r I ) × r given by
A = S S I Irr
where I Irr R r I × r is the extension of the ( r I × r I ) identity matrix by columns of zeros corresponding to reversible reactions. This means that, for a flux cone, the homogeneous linear inequalities are of a special type: part of these (m) are actually equalities defining a lower-dimensional subspace given by the nullspace of the stoichiometric matrix S , the others ( r I ) being non-negativity constraints regarding some single coordinate variables (given by the irreversible reactions) corresponding thus to particular half-spaces defined by such positive coordinate axes.
Conversely, to any convex polyhedral cone C in R r defined by a representation matrix A R n × r , we can associate a flux cone F C C of the same dimension in R r + n defined by
F C C = { x A x R r + n x C } = { v R r + n A I v = 0 , v i 0 for r + 1 i r + n } .
Using this correspondence (which defines a bijection of C onto F C C ), several properties, proven for flux cones, can actually be lifted to general convex polyhedral cones.
From the definition of a cone, for every nonzero element x of C, the whole half-line { α x α 0 } is contained in C. This is called a ray of C. Thus a flux vector is defined up to a positive scalar multiplication. The lineality space of C is the union of all lines of C, i.e., { x C x C } . If C is defined by a representation matrix A , its lineality space thus equals the nullspace of A , i.e., { x R r A x = 0 } . For a flux cone F C given by Equation (6), its lineality space is thus constituted by the flux vectors v such that v i = 0 for i Irr , i.e., flux vectors involving only reversible reactions (and thus the global flux can go in either one direction or the other). The cone C is called pointed if it does not contain a line, i.e., if its lineality space is reduced to { 0 } . For example, if C is contained in a closed orthant, it is pointed (where the 3 r closed orthants are defined by { x R r x i op i 0 for i = 1 , , r } for an operator vector op { , = , } r ). In particular, a flux cone F C with only irreversible reactions ( r I = r ) is necessarily pointed as it is included in the positive r-orthant (i.e., of dimension r). Actually, reversible flux vectors very rarely occur in metabolic networks, which therefore often give rise to pointed flux cones.

1.4. Extreme Vectors and Generating Sets

We are interested in finding an explicit (by extension) representation of C in the form of a (minimal) set of generators. A nonzero vector x C is called extreme (or extreme pathway [9,10] if x is a flux vector), if
x = x 1 + x 2 , with nonzero x 1 , x 2 C , implies x 1 = λ x 2 with λ > 0 .
If x C is extreme, then { α x α 0 } is called an extreme ray of C as all its nonzero elements are extreme (and thus for simplifying notations, we will not distinguish extreme vectors and extreme rays when it does not create confusion). In fact, C has an extreme ray if and only if C is pointed and, in this case, the extreme rays are the edges (faces of dimension one) of C and, according to Minkowski’s theorem, constitute the unique minimal (finite) set of generators of C for conical (i.e., non-negative linear) combination:
C = { k K β k y k β k 0 } = c o n e ( { y k } )
where the y k ’s, k K (finite index set), are representatives of the extreme rays (unique up to positive scalar multiplication) and c o n e is the conical hull. More precisely, we get an upper bound for the number of extreme vectors that are sufficient to decompose any given nonzero vector x C : x = k K x β k y k with | K x | m i n ( d i m ( C ) , | s u p p ( x ) | + | s u p p ( A x ) | ) . This result can be demonstrated first for a flux vector v of a pointed flux cone F C with an upper bound given by | K v | m i n ( d i m ( C ) , | s u p p ( v ) | ) and then for a vector x of a general pointed polyhedral cone C by using the correspondence (9) between C and F C C which maps the extreme vectors of C onto the extreme vectors of F C C . Extreme vectors of a pointed cone C will be noted ExVs.
The Double Description (DD) method [11,12], known as Fourier–Motzkin, is an incremental algorithm (which processes one by one each linear inequality ( A x ) j 0 ) to build an explicit description of a pointed cone C, as a minimal generating matrix (whose columns are in 1-to-1 correspondence with the extreme rays), from an implicit description of C by a representation matrix, i.e., to enumerate its extreme rays.
If C is not pointed, it is still finitely generated:
C = { k K β k y k + l L γ l z l β k 0 , γ l R } = c o n e ( { y k , z l , z l } )
with (not unique this time) minimal set of generators consisting of basis vectors z l of the lineality space and suitable vectors y k not in the lineality space (e.g., the extreme vectors of the pointed cone obtained by intersecting C with the orthogonal complement of its lineality space). Actually, Minkowski–Weyl theorem for cones states that it is equivalent for a set C to being a polyhedral cone (7) or to being a finitely generated cone, i.e., the conical hull of a finite set of vectors (as the y k ’s, z l ’s and z l ’s).

1.5. Elementary Vectors and Conformal Generating Sets

Now, if the existence and uniqueness of a minimal set of generators for conical decomposition in a pointed polyhedral cone C is satisfactory for an explicit geometric description of C, it is not in general meaningful for a flux cone F C representing the steady-state flux vectors of a metabolic network. In fact, for a metabolic pathway, only a conical decomposition without any cancellations is biochemically meaningful, as a reversible reaction cannot have a net rate in opposite directions in the contributing pathways. Indeed, the second law of thermodynamics states that a reaction can only carry flux in the direction of negative Gibbs free energy of the reaction, which is imposed by the values of the concentrations of the metabolites. This means that, when decomposing a flux vector, only so-called conformal sums, i.e., sums without cancellations, are biochemically admissible. A sum v = v 1 + v 2 of vectors is called conformal if, for all i { 1 , , r } :
v i = 0 implies v i 1 = v i 2 = 0 , v i > 0 implies v i 1 , v i 2 0 , v i < 0 implies v i 1 , v i 2 0 .
An equivalent definition is to define a sum v = v 1 + v 2 as conformal if
s i g n ( v 1 ) , s i g n ( v 2 ) s i g n ( v )
where the sign vector s i g n ( v ) { , 0 , + } r is defined by applying the sign function component-wise, i.e., s i g n ( v ) i = s i g n ( v i ) , and the partial order ≤ on { , 0 , + } r is defined by applying component-wise the partial order on { , 0 , + } induced by 0 < and 0 < + . For example, there is a one-to-one mapping between closed orthants O and sign vectors η , defined by O = { x R r s i g n ( x ) η } (O will be called defined by η and noted O η ), which induces a one-to-one mapping between closed r-orthants and full support (i.e., with only nonzero components) sign vectors. For ξ , η sign vectors and v a vector, we say that ξ conforms to η if ξ η and that v conforms to η if s i g n ( v ) η . We call two vectors v 1 , v 2 conformal if s i g n ( v 1 ) , s i g n ( v 2 ) η for a certain sign vector η or, equivalently, if v i 1 v i 2 0 for all i (so, v 1 + v 2 is a conformal sum if and only if v 1 and v 2 are conformal). A conformal sum v = v 1 + v 2 , i.e., verifying Equation (14), will be noted v = v 1 v 2 . It is therefore natural to look for generators as conformally non-decomposable vectors, where a nonzero vector x of a convex polyhedral cone C is called conformally non-decomposable, if
x = x 1 x 2 , with nonzero x 1 , x 2 C , implies x 1 = λ x 2 with λ > 0 .
A vector x (resp., flux vector v ) of a convex polyhedral cone C (resp., a flux cone F C ) is called elementary [7] if it is conformally non-decomposable. All nonzero elements of the ray defined by an elementary vector are elementary, i.e., elementary vectors are unique up to positive scalar multiplication (and thus we will often not distinguish elementary vectors and elementary rays). Elementary vectors (resp., flux vectors) will be noted EVs (resp., EFVs).
The elementary rays constitute the unique minimal (finite) set of conformal generators of C, i.e., generators for conformal conical sum:
C = { g V G β g e g β g 0 } = c o n e ( { e g } )
where the e g ’s, g V G (finite index set), are representatives of the elementary rays (unique up to positive scalar multiplication) and c o n e is the conical conformal hull. More precisely, we get an upper bound for the number of elementary vectors that are sufficient to decompose any given nonzero vector x C : x = g V G x β g e g with | V G x | m i n ( d i m ( C ) , | s u p p ( x ) | + | s u p p ( A x ) | ) . This result can be demonstrated first for a flux vector v of a flux cone F C with an upper bound given by | V G v | m i n ( d i m ( C ) , | s u p p ( v ) | ) and then for a vector x of a general polyhedral cone C by using the correspondence (9) between C and F C C which maps the elementary vectors of C onto the elementary flux vectors of F C C .
It follows from (10) and (15) that an extreme vector of C is elementary, i.e., ExVs ⊆ EVs, but the converse is generally false as a conformally non-decomposable vector may be conically decomposable. Nevertheless, if C is contained in a closed orthant, there is identity between extreme vectors and elementary vectors, i.e., ExVs = EVs. More precisely, for nonzero x C and O a closed orthant with x O , then x is elementary in C if and only if it is extreme in C O . It results that the elementary vectors of C are the extreme vectors of intersections of C with any closed orthant:
E V s ( C ) = O orthant E x V s ( C O ) .
Elementary vectors of C can thus be obtained by using algorithms, such as DD, to compute extreme vectors ExVs of the pointed polyhedral cones C O . By doing this, it is convenient to select in Equation (17) only a minimal subset of closed orthants O in order to avoid equality or inclusion between the C O ’s (nonempty intersection can obviously not be avoided as orthants are closed). It is clearly enough to consider only the 2 r closed r-orthants of maximal dimension r, but this does not avoid equality or inclusion. Let { η i } be the maximal (for the partial order defined above on { , 0 , + } r ) sign vectors of s i g n ( C ) = { s i g n ( x ) x C } . It is then enough in Equation (17) to consider the closed orthants O i = { x R r s i g n ( x ) η i } , and there is no equality or inclusion between the C O i ’s, where C O i = C η i = { x C s i g n ( x ) η i } . The C η i ’s are called topes, noted Ts (flux topes [13] noted FTs for a flux cone F C ). C is thus decomposed into topes and Equation (17) can be rewritten as
E V s ( C ) = η i maximal in s i g n ( C ) E x V s ( C η i ) .
Note that, for a flux cone F C Equation (6), a FT is defined by specifying a maximal subset of reactions with fixed directions (thus fixing the directions of reversible reactions), the others having necessarily a zero flux. This simplifies if F C is consistent, i.e., without unused reaction, which means that every reaction, in every possible direction for reversible reactions, is supported by a flux vector: i { 1 , , r } v F C v i > 0 and i { 1 , , r } \ Irr v F C v i < 0 . We can always assume F C consistent after a preprocessing step (practically, this can be achieved by using flux variability analysis [14]) that removes all reactions that cannot carry nonzero steady-state flux and changes all reversible reactions that cannot carry flux in both directions into irreversible ones. In this case, every remaining reaction in every possible direction is supported by a flux vector with full support (i.e., with nonzero flux in any reaction) and all FTs F C η i have full support, i.e., the η i ’s have full support or, equivalently, the O i ’s are r-orthants. An obvious upper bound for the number of FTs is thus 2 r r I .
Now, for a flux cone F C , another commonly used method is, at the extreme opposite, to have it included into a single (positive) orthant in a higher dimension by splitting each reversible reaction i into a forward i + and a backward i irreversible reaction. This means decomposing a flux in i as v i = v i + v i with v i + = v i + 0 and v i = v i 0 . Columns of the stoichiometric matrix S corresponding to reversible reactions i are negated (which means exchanging the roles of substrates and products in i) and appended to S as new columns to form the new stoichiometric matrix S ˜ R m × r ˜ , where r ˜ = 2 r r I is the new number of reactions and all r ˜ reactions are now irreversible, Irr ˜ = { 1 , , r ˜ } . F C is in one-to-one correspondence with vectors v of F C ˜ such that v i + . v i = 0 for i reversible. In particular, the fluxes of the form v i + = v i > 0 with all other components being null are obtained as extreme vectors of F C ˜ but represent futile cycles (involving reactions i + and i ) without biological reality and must be eliminated. Finally, EFVs ( F C ) are in one-to-one correspondence with ExVs ( F C ˜ ) ∖{futile cycles} (called at the origin extreme currents in stoichiometric network analysis [15]). F C ˜ is included in the positive r ˜ -orthant and has thus only one FT.
We therefore have two opposite ways of dealing with reversible reactions for computing EFVs of a flux cone F C : either splitting each reversible reaction into two irreversible ones, such that F C is reduced to a single FT at the price of an increase in the space dimension by r r I (which can cause serious efficiency problems to algorithms such as DD) or keeping the reversible reactions unchanged and decomposing F C into FTs, in each of which the directions of reversible reactions are fixed, at the price of a potentially exponential (in terms of r r I ) number of FTs to consider. All intermediate cases, where only a subset of reversible reactions are split into irreversible ones and the others are processed by decomposition into FTs, are obviously possible. Independently of the solution adopted, we will work most of the time in a given FT for F C , defined by a (full support if F C is consistent) sign vector η , and the EFVs of F C which conform to η are thus given by the ExVs of this FT F C η .

1.6. Elementary Modes

The null value 0 plays a component-wise crucial role in definitions of the support of a vector (4), of a flux cone (6) and of a conformal sum Equations (13) and (14). A close relationship results between support-minimal vectors and elementary vectors in a flux cone. A nonzero vector x C is called support-minimal, if
s u p p ( x ) s u p p ( x ) , with nonzero x C , implies s u p p ( x ) = s u p p ( x ) .
A nonzero vector x (resp., flux vector v ) of a convex polyhedral cone C (resp., a flux cone F C ) is called elementary mode (resp., elementary flux mode) if it is support-minimal (the concept of elementary flux mode was first introduced, under the name of elementary vector [16], for a subspace of R r , i.e., for a flux linear subspace F S Equation (3) and then [17] for a flux cone F C with a definition actually closer to that of a support-wise non-decomposable vector). All nonzero elements of the ray defined by an elementary mode are elementary modes having the same support, i.e., elementary modes are unique up to positive scalar multiplication (and we will therefore in general identify two positively proportional elementary modes). Elementary modes (resp., flux modes) will be noted EMs (resp., EFMs). The concept of EFM in a flux cone F C is biologically significant as it represents a minimal pathway operating in a steady state, i.e., with all reactions involved necessarily active (with a nonzero net rate), which means that no proper sub-pathway can operate in a steady state.
Note that if x is an EFM, then s i g n ( x ) is a minimal (for the partial order defined above on { , 0 , + } r ) nonzero element of s i g n ( F C ) and, conversely, it is shown that a minimal nonzero sign vector σ s i g n ( F C ) determines an EFM x with s i g n ( x ) = σ by F C σ = { v F C s i g n ( v ) σ } = { v F C s i g n ( v ) = σ } { 0 } = { λ x λ 0 } . There is thus a one-to-one mapping between EFMs and minimal nonzero sign vectors of s i g n ( F C ) . Comparing with the one-to-one mapping between FTs and maximal sign vectors of s i g n ( F C ) , we see that EFMs and FTs are dual concepts.
Now, for a flux cone F C , support-minimality and conformal non-decomposability are equivalent properties, i.e., there is identity between elementary flux modes and elementary flux vectors: EFMs = EFVs. For metabolic pathways in a flux cone F C , there is therefore identity between minimal (for support inclusion) pathways and non-decomposable (for conformal sum) pathways. From Equation (16), the EFMs constitute a conformal generating set (i.e., generating set for conformal sum) for F C , and in fact the unique minimal such set (for that matter one way of proving Equation (16) for flux cones F C is to prove it with EFMs as a conformal generating set and to prove that EFMs = EFVs). From Equation (18), for any maximal sign vector η of s i g n ( F C ) , the EFMs of F C which conform to η are the EFMs of the FT F C η and coincide with the ExVs of the said FT, this result being the basis of methods for computing EFMs [18,19]. EFMs are thus decomposed into subsets according to the decomposition of F C into flux topes [13]: the EFMs of the FT F C η correspond to the F C σ ’s, where the σ ’s are the minimal nonzero sign vectors of s i g n ( F C ) such that σ η .
For a general polyhedral cone C, there is no direct relationship between elementary modes and elementary vectors. Nevertheless, from Equation (16), it follows that, for any EM in C, there is an EV with the same support. Thus, all minimal support patterns of vectors appear in the set of supports of elementary vectors and are actually the minimal elements in this set for subset inclusion:
s u p p ( E M s ) = Min { s u p p ( E V s ) } .
In addition, for a general polyhedral cone C, the correspondence (9) between C and the higher dimensional flux cone F C C maps the elementary vectors of C onto the elementary flux vectors of F C C , i.e., the elementary flux modes of F C C :
E V s ( C ) = { x R r x y E F M s ( F C C ) }
with F C C given by Equation (9).
Moreover, it remains true that any vector which is support-minimal in a given C η is actually support minimal in C, as it depends only on the convexity of C: if x and x are vectors in a convex set with s u p p ( x ) s u p p ( x ) , then a vector x exists in this convex set with s i g n ( x ) s i g n ( x ) and s u p p ( x ) s u p p ( x ) (we take x = λ x + ( 1 λ ) x with λ minimal in ( 0 , 1 ] such that x i = 0 for a certain i with x i 0 ). Thus, EMs of C can be computed tope by tope:
E M s ( C ) = η i maximal in s i g n ( C ) E M s ( C η i ) .
See Figure 1 for an illustration of the concepts of extreme vectors, elementary vectors and elementary modes in a convex polyhedral cone.

1.7. Inhomogeneous Linear Constraints and Polyhedra

Additionally, in this standard setting, fluxes may be constrained by other constraints, typically lower and upper bounds regarding reaction rates:
v i v i v i +
or, more generally, any set of inhomogeneous linear constraints, noted ILC , that can be written in the general form
G v h
where G R l × r is a matrix and h R l a vector with nonzero components, defining a general inhomogeneous convex polyhedron P ILC = { v R r G v h } . The solution space S o l S , Irr , ILC is thus now a s-polyhedron or flux polyhedron noted F P and defined by
S o l S , Irr , ILC = F P = { v R r S v = 0 , v i 0 for i Irr , G v h } = F C P ILC .
This is a particular case of (general) convex polyhedron that is defined implicitly by finitely many linear inequalities:
P = { x R r A x b }
with a suitable matrix A R n × r and vector b R n and is thus the intersection of n (affine) half-spaces ( R r is equipped with both its structure of affine space, with origin 0, and its underlying structure of vector space and we will identify a point in the affine space with the corresponding vector in the vector space). Its dimension is defined as the dimension of its affine span. In this way, a flux polyhedron F P corresponds to a particular matrix A R ( 2 m + r I + l ) × r and vector b R 2 m + r I + l given by
A = S S I Irr G b = 0 0 0 h
meaning that the inequalities that are homogeneous actually divide into m equalities defining a lower-dimensional subspace given by the nullspace of the stoichiometric matrix S and into r I non-negativity constraints regarding certain single coordinate variables (given by the irreversible reactions). A polyhedral cone (resp., flux cone) is a special case of polyhedron (resp., flux polyhedron) where b = 0 (resp., ILC = ).
To any nonempty polyhedron P given by Equation (26) is associated its so-called recession cone C P = { x R r A x 0 } , which is the polyhedral cone containing all unbounded directions (rays) of P (if P is a polyhedral cone, then P = C P ). A bounded polyhedron is called a polytope and thus P is a polytope if and only if its recession cone is trivial: C P = { 0 } . P is called pointed if its recession cone C P is pointed, i.e., if its lineality space { x R r A x = 0 } , also called the lineality space of P (as it contains all unbounded lines of P), is trivial. For a flux polyhedron F P , we have: C F P = F C C P ILC (thus F P can be a polytope without P ILC being so and be pointed without either F C or P ILC being so). Note that C F P is not in general a flux cone.

1.7.1. Extreme Points and Vectors and Generating Sets

A vector x P is called an extreme point, if it cannot be written as a convex combination of two distinct vectors of P:
x = λ x 1 + ( 1 λ ) x 2 , with x 1 , x 2 P and 0 < λ < 1 , implies x 1 = x 2 .
Note that, as P is convex, it is enough to consider the midpoint of x 1 and x 2 , i.e., to take λ = 0.5 . Extreme points coincide with vertices of P, where a vertex of P is defined as a face of dimension 0. P is pointed if and only if it has a vertex and in this case, according to Minkowski’s theorem, the vertices of P and the extreme rays of C P constitute the unique minimal (finite) set of “bounded” and “unbounded” generators of P for convex and conical combination, respectively:
P = { j J α j p j + k K β k y k α j , β k 0 , j J α j = 1 } = c o n v ( { p j } ) + c o n e ( { y k } )
where the p j ’s, j J (finite index set), are the extreme points (vertices) of P, noted ExPs, and the y k ’s, k K (finite index set), are the extreme vectors of C P (unique up to positive scalar multiplication), noted ExVs, and c o n v is the convex hull (if P is a pointed polyhedral cone, then it has only one vertex, which is the zero vector, and the Equation (29) reduces to Equation (11)). More precisely, we get an upper bound for the number of extreme points and vectors that are sufficient to decompose any given vector x P : x = j J x α j p j + k K x β k y k with | J x | + | K x | m i n ( d i m ( P ) + 1 , | s u p p ( x ) | + | s u p p ( A x b ) | + 1 ) . This result can be deduced from result (11) for a pointed flux cone F C by using the following correspondence between P and such a flux cone.
In fact, to any convex polyhedron P in R r defined by a matrix A R n × r and vector b R n Equation (26), we can associate a flux cone F C P in a higher dimension R r + 1 + n defined by
F C P = { x ξ A x ξ b R r + 1 + n ξ 0 , A x ξ b 0 } = { v R r + 1 + n A b I v = 0 , v i 0 for r + 1 i r + 1 + n } .
This introduces a correspondence between vectors x of P and vectors x 1 A x b of F C P , which maps vertices of P onto extreme vectors of F C P with component ξ = 1 , and between vectors x of C P and vectors x 0 A x of F C P , which maps extreme vectors of C P onto extreme vectors of F C P with component ξ = 0 . Thanks to this correspondence, several properties, proven for flux cones, can be lifted to general convex polyhedra.
For the particular case where P is a flux polyhedron F P given by Equation (25), the correspondence (30) simplifies by associating to F P the flux cone F C F P in dimension R r + 1 + l defined by
F C F P = { ( v ξ G v ξ h ) R r + 1 + l S v = 0 , v i 0 for i Irr , ξ 0 , G v ξ h 0 } = { v R r + 1 + l S 0 0 G h I v = 0 , v i 0 for i Irr and for r + 1 i r + 1 + l } .
The DD method builds an explicit description of a pointed polyhedron P, in the form of two generating matrices whose columns are respectively the p j ’s and the y k ’s, from an implicit description of P as in Equation (26), i.e., enumerates its vertices ExPs and extreme vectors ExVs.
If P is not pointed, it is still finitely generated:
P = { j J α j p j + k K β k y k + l L γ l z l α j , β k 0 , γ l R , j J α j = 1 } = c o n v ( { p j } ) + c o n e ( { y k , z l , z l } )
with a (not unique this time) minimal set of generators consisting of basis vectors z l of the lineality space and suitable vectors p j and y k (e.g., the vertices and extreme vectors of the pointed polyhedron obtained by intersecting P with the orthogonal complement of its lineality space; if P is a non-pointed polyhedral cone, there is no nonzero p j and Equation (32) reduces to Equation (12)). In fact, Minkowski–Weyl theorem for polyhedra states that it is equivalent for a set P to be a polyhedron (26) or to be finitely generated, i.e., to be the Minkowski sum of the convex hull of a finite set of vectors (as the p j ’s) and of the conical hull of a finite set of vectors (as the y k ’s, z l ’s and z l ’s).

1.7.2. Elementary Points and Vectors and Conformal Generating Sets

However, such a decomposition into a finite set of generators is not in general satisfactory for a flux polyhedron as only a decomposition without any cancellations is biochemically meaningful, as was stipulated for a flux cone in Section 1.5. In the same way as we replaced, as generators for a polyhedral cone, extreme vectors by conformally non-decomposable vectors, we will replace, as generators for a polyhedron P, extreme points (vertices) by convex-conformally non-decomposable vectors (and, for its recession cone C P , extreme vectors by conformally non-decomposable vectors). A vector x of a polyhedron P is called convex-conformally non-decomposable, if
x = λ x 1 ( 1 λ ) x 2 , with x 1 , x 2 P and 0 < λ < 1 , implies x 1 = x 2 .
Given a polyhedron P (resp., flux polyhedron F P ), a vector (resp., flux vector) x is called an elementary point (resp., elementary flux point)—also called “bounded” elementary vector—of P if x P is convex-conformally non-decomposable and is called an elementary vector (resp., elementary flux vector)—also called “unbounded” elementary vector—of P if x C P is conformally non-decomposable (it is unique only up to positive scalar multiplication) [7]. Elementary points (resp., flux points) will be noted EPs (resp., EFPs) and elementary vectors (resp., flux vectors) will be noted EVs (resp., EFVs), which is consistent with the same notation for polyhedral cones and flux cones. We will note Es = EPs ∪ EVs (resp., EFs = EFPs ∪ EFVs) the elementary elements (resp., elementary fluxes) of P (resp., F P ).
The elementary points and the elementary rays constitute the unique minimal (finite) set of conformal generators of P, i.e., generators for convex-conformal (for elementary points) and conformal (for elementary vectors) sum:
P = { g P G α g e g g V G β g e g α g , β g 0 , g P G α g = 1 } = c o n v ( { e g g P G } ) c o n e ( { e g g V G } )
where the e g ’s, g P G (finite index set), are the elementary points and the e g ’s, g V G (finite index set), are the elementary vectors (unique up to positive scalar multiplication), and c o n v is the convex conformal hull. More precisely, we get an upper bound for the number of elementary points and vectors that are sufficient to decompose any given vector x P : x = g P G x α g e g g V G x β g e g with | P G x | + | V G x | m i n ( d i m ( P ) + 1 , | s u p p ( x ) | + | s u p p ( A x b ) | + 1 ) . This result can be demonstrated from result (16) for a flux cone by using the correspondence (30) between P and F C P which maps the elementary points of P onto the elementary flux vectors of F C P with component ξ = 1 and the elementary vectors of P onto the elementary flux vectors of F C P with component ξ = 0 .
We already know that an extreme vector of C P is elementary, and is therefore an elementary vector of P, i.e., ExVs ⊆ EVs. It follows from Equations (28) and (33) that an extreme point (vertex) of P is an elementary point of P, i.e., ExPs ⊆ EPs, but the converse is generally false as a convex-conformally non-decomposable vector may be convex decomposable. Nevertheless, if P is contained in a closed orthant (and thus C P too), any sum of vectors of P (resp., C P ) is conformal and thus there is identity between extreme points (vertices) and elementary points (resp., between extreme vectors and elementary vectors), i.e., ExPs = EPs and ExVs = EVs. More precisely, for x P (resp., x C P and nonzero) and O a closed orthant with x O , then x is an elementary point (resp., elementary vector) of P if and only if it is a vertex in P O (resp., an extreme vector in C P O = C P O ). It follows that the elementary points (resp., elementary vectors) of P are the vertices (resp., extreme vectors) of intersections of P (resp., C P ) with any closed orthant, which are pointed subpolyhedra (resp., pointed subcones):
E P s ( P ) = O orthant E x P s ( P O ) E V s ( P ) = O orthant E x V s ( C P O ) .
Note in particular that, if 0 P , then 0 EPs. Elementary points and vectors can therefore be obtained by using algorithms, such as DD, to compute vertices ExPs and extreme vectors ExVs of the pointed polyhedra P O . It is obvious that considering only the 2 r closed r-orthants is enough. Now, as for polyhedral cones, decomposing the polyhedron into topes is better:
E P s ( P ) = η i maximal in s i g n ( P ) E x P s ( P η i ) E V s ( P ) = η j maximal in s i g n ( C P ) E x V s ( C P η j ) .
If O i is the closed orthant defined by η i , then the corresponding tope for P is P O i = P η i = { x P s i g n ( x ) η i } and, as seen for polyhedral cones, C P O j = C P η j is a tope for the recession cone C P . Examine how the equality C P η = C P η , for η an arbitrary sign vector, can be expressed in terms of topes. Note first that any η j is dominated by at least one η i for the partial order ≤ on { , 0 , + } r : η j η i η j η i , which means that O j is a sub-orthant of O i , and we have C P η j = C P η i , expressing the relation between the topes for the recession cone of the polyhedron and the recession cones of certain of the polyhedron topes (precisely those topes P η i for which η i dominates an η j , necessarily unique). More generally, the recession cone of any tope P η i for P can be expressed as a subcone of a tope for the recession cone of P by C P η i = C P c ( η i ) , where c ( η i ) = m a x { η s i g n ( C P ) η η i } is the greatest (it is necessarily unique) sign vector in s i g n ( C P ) dominated by η i (thus, if η i does not dominate any η j , C P c ( η i ) is not a tope for C P ; this is the case for example if P is not a polytope but P η i is, implying that c ( η i ) = 0 and that C P c ( η i ) = { 0 } is not a tope for C P { 0 } ).
For the particular case of a flux polyhedron F P (25), we have F P η = F C η P ILC and C F P η = C F P η = F C η C P ILC , for any sign vector η . Any flux tope F P η i for F P can thus be expressed as F P η i = F C η k P ILC , for a certain flux tope F C η k for FC, i.e., by applying the constraints ILC to a flux tope for F C . The same holds for the recession cone: C F P η i = C F P c ( η i ) = F C η k C P ILC , i.e., by applying the homogeneous counterparts of constraints ILC to a flux tope for F C . If F P is assumed to be consistent (same definition as for a flux cone, i.e., without unused reaction, always with a zero flux), all FTs F P η i then have full support, i.e., the η i ’s have full support or, equivalently, the O i ’s are r-orthants. An obvious upper bound for the number of FTs is thus 2 r r I . Note that C F P to be consistent is a sufficient (but not necessary) condition for F P to be consistent, and that, if F P is consistent, so is F C (but F P can be inconsistent even if both F C and P ILC are consistent).
The method used to include a flux cone into a single (positive) orthant in higher dimension, by splitting each reversible reaction i into a forward i + and a backward i irreversible reaction, applies as well to a flux polyhedron F P . Matrix G is extended into matrix G ˜ R l × r ˜ in the same way that S is extended into S ˜ R m × r ˜ , where r ˜ = 2 r r I is the new number of reactions and all r ˜ reactions are now irreversible, Irr ˜ = { 1 , , r ˜ } . F P is in one-to-one correspondence with vectors v of F P ˜ such that v i + . v i = 0 for i reversible. In particular, the fluxes with v i + , v i > 0 for a certain i, which are obtained as vertices or extreme vectors of F P ˜ , are futile (involving a net rate in opposite directions in reactions i + and i ) without biological reality and must be eliminated (they are not generally limited to futile cycles as in flux cones). Finally, the set of elementary fluxes EFs ( F P ) is in one-to-one correspondence with (ExPs ( F P ˜ ) ∪ ExVs ( F P ˜ )) ∖{futile fluxes}. As F P ˜ is included in the positive r ˜ -orthant, it has only one FT.
As for flux cones, there are two opposite ways of dealing with reversible reactions for computing EFs = EFPs ∪ EFVs of a flux polyhedron F P : either splitting each reversible reaction into two irreversible ones, reducing F P to a single FT in higher dimension, or keeping the reversible reactions unchanged and decomposing F P into FTs without increasing the dimension; all intermediate cases are possible. Whatever the solution adopted, EFs are obtained as the union of ExPs and ExVs for each flux tope for F P and thus we will generally work in a given FT for F P , defined by a (full support if F P is consistent) sign vector η , and the EFs of F P which conform to η are thus given by the ExPs and ExVs of this FT F P η .

1.7.3. Elementary Modes

For a general polyhedron P (and even for a flux polyhedron F P ), there is no direct relationship between support-minimal Equation (19) vectors and elementary elements as it was the case for flux cones F C . Nevertheless, from Equation (34), it follows that, for any support-minimal nonzero vector in P (still called elementary mode EM), there is an elementary element with the same support. Thus, Equation (20) generalizes to the case of polyhedra and all minimal support patterns of vectors appear in the set of supports of elementary elements and are actually the minimal nonempty elements in this set for subset inclusion:
s u p p ( E M s ) = Min ( { s u p p ( E s ) } \ { } ) .
In particular, for a metabolic network, this means that any set of reactions involved in a minimal pathway (EFM) appears as the set of reactions involved in a certain elementary flux (EF).
In addition, for a general polyhedron P, the correspondence (30) between P and the higher dimensional flux cone F C P maps the elementary points (resp., elementary vectors) of P onto the elementary flux vectors, i.e., the elementary flux modes, of F C P with component ξ = 1 (resp., with component ξ = 0 ):
E P s ( P ) = { x R r x 1 y E F M s ( F C P ) } E V s ( P ) = { x R r x 0 y E F M s ( F C P ) }
with F C P given by Equation (30). The same formula holds for the particular case of a flux polyhedron F P (25) with F C F P given by Equation (31).
Moreover, as seen in the proof of Equation (22), any vector which is support-minimal in a certain P η is actually support minimal in P, thus EMs of P can be computed tope by tope:
E M s ( P ) = η i maximal in s i g n ( P ) E M s ( P η i ) .
See Figure 2 for an illustration of the concepts of extreme points and vectors, elementary points and vectors and elementary modes in a convex polyhedron.

1.8. Complexity Results

Enumerating the vertices of an unbounded polyhedron has been proven to be an NP-hard enumeration problem [20]. However, the hardness of vertex generation for bounded polyhedra, and also the complexity of enumerating together vertices and extreme rays of polyhedra, are open problems. Consequently, the complexity of enumerating all EFs or EFMs of a metabolic network remains an open problem [21]. Nevertheless, enumerating all extreme rays of a flux cone that contain a given reaction in their support has been proven not to be in polynomial total time unless P = NP [21], which means that the output cannot be generated in time polynomial in the combined size of the input and the output. The problem of enumerating the extreme vectors of a polyhedral cone or the vertices of a polytope with a polynomial delay (i.e., the time between one output item and the next is bounded by a polynomial function in terms of the input size), or the even weaker question whether this enumeration can be done in polynomial total time, is an open problem too. Therefore, the possibility of enumerating EFs or EFMs of a metabolic network with a polynomial delay is an open problem [22] (except for a flux cone with all reactions being reversible, i.e., a flux linear subspace (3), as the EFMs are then the circuits of a linear matroid [16]). Given a flux cone and two reactions, deciding if there exists an extreme ray of the cone that has both reactions in its support is NP-complete [21]. Thus, given a metabolic network, deciding if an EF or EFM exists with a given support of size at least two is an NP-complete problem. The same result holds for deciding if an EF or EFM exists whose support size is bounded above by a given positive integer [22].
Regarding now the number of EFs or EFMs, counting the extreme rays of a polyhedral cone or the vertices of a polytope has been proven to be a #P-complete problem [23]. Thus, counting the number of EFs or EFMs of a metabolic network is also a #P-complete problem [22]. The McMullen’s upper bound theorem [24] states that, for any fixed positive integers d and n, the maximum number of j-faces of a d-polytope with n facets (i.e., faces of dimension d 1 ) is attained by the dual cyclic polytope c ( d , n ) for all j = 0 , 1 , , d 2 . A consequence is that the maximum number of vertices of a d-polytope with n facets is given by n d / 2 n d + n d / 2 1 n d 2 n d / 2 n d . We thus obtain that the number of EFs or EFMs in a metabolic network (after having split each reversible reaction into two irreversible ones) is bounded above by a quantity approximately equal to 2 ( r ˜ + m ) / 2 m with r ˜ = 2 r r I (so r ˜ varies between r and 2 r ). If the number m of internal metabolites is small compared to the total number r ˜ of reactions (after having split reversible reactions), the number of EFs or EFMs is then bounded above by a quantity of order Θ ( ( r ˜ / 2 ) m ) . If m is close to r ˜ , this number is bounded above by a quantity of order Θ ( r ˜ ( r ˜ m ) / 2 ) . The worst case occurs when m is close to r ˜ / 3 with an upper bound approximately equal to 2 2 m m . We obtain the same results with r instead of r ˜ if we do not split reversible reactions but fix their signs arbitrarily, i.e., if we consider EFs or EFMs in an arbitrary closed r-orthant O, i.e., in an arbitrary flux tope for F C (or F P ). However, note that in practice the actual number of EFs or EFMs of a metabolic network is likely to be much smaller than these upper bounds.

2. Metabolic Pathways in the Presence of Biological Constraints

Although the current improved implementations of the DD method [25,26] allow the computation of millions, even billions, of EFMs or EFs, tackling genome-scale metabolic models (GSMMs) is still beyond our reach. Moreover, most of the computed EFMs or EFs are not biologically valid, because only stoichiometry and certain flux constraints, such as irreversibility of reactions or bounds on reaction rates, are taken into account. For both scaling up to large systems and limiting the number of biologically invalid solutions, it is necessary to consider additional biological constraints, such as thermodynamic, kinetic or regulatory constraints.

2.1. Biological Constraints

2.1.1. Thermodynamic Constraints

Assuming constant pressure and a closed system, according to the second law of thermodynamics, a reaction i proceeds spontaneously only in the direction of its negative Gibbs free energy Δ r G i [27], given by
Δ r G i = Δ r G i 0 + R T ln ( j M j S ¯ j i )
where Δ r G i 0 is the standard Gibbs free energy of reaction i, R the molar gas constant, T the absolute temperature, M j the (positive) concentration of metabolite j and S ¯ j i is the stoichiometric coefficient of metabolite j in reaction i (i.e., S ¯ R ( m + m ¯ ) × r is the extension of the stoichiometric matrix S to the m ¯ external metabolites). This means that Δ r G i < 0 , (resp., Δ r G i > 0 ) is a necessary condition to get v i > 0 (resp., v i < 0 ), which can be expressed by the constraint s i g n ( v i ) s i g n ( Δ r G i ) , and that a flux vector v is thermodynamically feasible [28] if and only if all its components v i satisfy such a constraint (it is enough to consider those i s u p p ( v ) as the constraint is trivially satisfied when v i = 0 ). The thermodynamic constraint for v , that depends on the vector M of metabolite concentrations, can thus be defined as
TC M ( v ) = i s u p p ( v ) , s i g n ( v i ) = s i g n ( Δ r G i 0 + R T ln ( j M j S ¯ j i ) ) .
As at equilibrium Δ r G i is null, we obtain Δ r G i 0 = R T ln ( K e q i ) , where K e q i = j M j e q S ¯ j i is the equilibrium constant of reaction i. Thus, Δ r G i can be rewritten as Δ r G i = R T ln ( j M j S ¯ j i / K e q i ) and the thermodynamic constraint TC M ( v ) as
TC M ( v ) = i s u p p ( v ) , s i g n ( v i ) = s i g n ( j M j S ¯ j i K e q i ) .
Often, the concentrations of external metabolites can be measured and included in the constraint as known parameters, keeping only a dependency of the constraint on the concentrations of internal metabolites. The formula Δ r G i = R T ln ( j M j S ¯ j i / K e q i ) can be rewritten, by dividing the numerator and denominator of the fraction by the terms dealing with external metabolites, as Δ r G i = R T ln ( j M ̲ j S j i / K ^ e q i ) , where K ^ e q i = K e q i / j M ¯ j S ¯ j i is the apparent equilibrium constant of the reaction i and M ̲ (resp., M ¯ ) the vector of internal (resp., external) metabolite concentrations. The thermodynamic constraint (42) can thus be rewritten as
TC ̲ M ̲ ( v ) = i s u p p ( v ) , s i g n ( v i ) = s i g n ( j M ̲ j S j i K ^ e q i ) .
For given metabolite concentrations vector M (resp., internal metabolite concentrations vector M ̲ ), let ts M { , 0 , + } r (resp., ts M ̲ ) be the fixed thermodynamic sign vector defined by:
( ts M ) i = s i g n ( Δ r G i 0 + R T ln ( j M j S ¯ j i ) ) = s i g n ( j M j S ¯ j i K e q i ) ( ts M ̲ ) i = s i g n ( j M ̲ j S j i K ^ e q i ) for 1 i r .
Then, the thermodynamic constraint TC M (resp., TC ̲ M ̲ ) can be rewritten as
TC M ( v ) = s i g n ( v ) ts M TC ̲ M ̲ ( v ) = s i g n ( v ) ts M ̲ .
Thus, the set S o l TC M (resp., S o l TC ̲ M ̲ ) of vectors v satisfying the constraint TC M ( v ) (resp., TC ̲ M ̲ ( v ) ), given by { v R r s i g n ( v ) ts M } (resp., { v R r s i g n ( v ) ts M ̲ } ), is the closed orthant O ts M (resp., O ts M ̲ ) defined by ts M (resp., ts M ̲ ), of dimension r if ts M (resp., ts M ̲ ) has full support, i.e., ( ts M ) i 0 (resp., ( ts M ̲ ) i 0 ) for all i, of lesser dimension otherwise.
Lemma 1.
Given metabolite concentrations M (resp., internal metabolite concentrations M ̲ ), the set S o l TC M (resp., S o l TC ̲ M ̲ ) of vectors in R r satisfying the thermodynamic constraint TC M (resp., TC ̲ M ̲ ) is the set of vectors that conform to the thermodynamic sign vector ts M (resp., ts M ̲ ) Equation (44), i.e., the closed orthant O ts M (resp., O ts M ̲ ) defined by ts M (resp., ts M ̲ ).

2.1.2. Kinetic Constraints

Metabolic reactions are catalyzed by enzymes. The catalytic mechanisms of key enzymes have been investigated in great detail and described by mathematical formulas. However, many kinetic equations are still unknown and have to be substituted by standard rate laws such as mass-action kinetics, power laws, reversible Hill kinetics, lin-log kinetics, convenience kinetics, generic rate equations or TKM rate laws [29]. What is important for our study is that these modular rate laws share the general form
v i = E i κ i ( M ) with κ i = f i r e g T i D i + D i r e g
where E i is the (non-negative) level of the enzyme catalyzing the reaction i (given either as an amount or as a concentration, in which case the rate law is pre-multiplied by the compartment volume) and κ i depends on the concentrations of the metabolites occurring in i (reactants of i) and on reaction i stoichiometry, rate law considered, allosteric regulation and parameters. In the general form of κ i , T i is the thermodynamic numerator (which can be written in the compact form k i + θ i + k i θ i with turnover rate parameters k i ± ) that gives its sign to κ i (and thus to v i ) and reflects the relationship between chemical potentials and reaction directions and ensures that the rate vanishes at chemical equilibrium, f i r e g and D i r e g , both positive, implement enzyme regulation (partial or complete for the first one, specific for the second) and D i is the (positive) kinetic denominator, a polynomial of scaled reactant concentrations whose terms correspond to different binding states of the enzyme (reducing the enzyme amount available for catalysis), which depends on the rate law considered.
The kinetic constraint for v depends both on the vector E of enzyme concentrations and on the vector M of metabolite concentrations and can thus be defined as
KC E , M ( v ) = v = E κ ( M )
where ∘ is the component-wise product of vectors: ( E κ ) i = E i κ i . This means that the flux vector is a component-wise linear function of the vector of enzyme concentrations (and a nonlinear function of the metabolite concentrations vector). For nonzero E i , the sign of T i gives the direction in which reaction i proceeds, so in this sense the kinetic constraint includes the thermodynamic constraint.
This can be highlighted on a widely-used rate law, the reversible Michaelis–Menten kinetics [30]. In the simple case of a reaction i where the enzyme can only exist in one of three distinct states: free, all substrates bound, or all products bound, it can be written as
κ i = k i + j S ¯ j i < 0 ( M j / K i j M ) S ¯ j i k i j S ¯ j i > 0 ( M j / K i j M ) S ¯ j i 1 + j S ¯ j i < 0 ( M j / K i j M ) S ¯ j i + j S ¯ j i > 0 ( M j / K i j M ) S ¯ j i
where k i + and k i are the maximal forward and backward rates in reaction i per unit of enzyme and the K i j M ’s are the Michaelis constants. Equating the numerator to zero at equilibrium, we obtain k i + k i j ( K i j M ) S ¯ j i = K e q i . This gives
κ i = k i + j S ¯ j i < 0 ( M j / K i j M ) S ¯ j i 1 + j S ¯ j i < 0 ( M j / K i j M ) S ¯ j i + j S ¯ j i > 0 ( M j / K i j M ) S ¯ j i 1 j M j S ¯ j i K e q i
that is, the product of three terms: the positive capacity term per unit of enzyme, the positive (smaller than one) fractional saturation term depending on M and the thermodynamic term, which can be rewritten as 1 e Δ r G i / R T and gives the sign of κ i (and thus the sign of v i ): κ i > 0 Δ r G i < 0 , i.e., the thermodynamic constraint. Thus, s i g n ( κ ( M ) ) = ts M and we will assume that this equality holds for all kinetic laws we consider.
Note that, for given metabolite and enzyme concentrations M and E , the kinetic constraint KC E , M defines completely and uniquely the only vector that satisfies it: S o l KC E , M = { E κ ( M ) } .

2.1.3. Regulatory Constraints

Coupling metabolic networks with Boolean transcriptional regulatory networks allows us to express the additional constraints imposed by gene regulatory information on a metabolic network and to take them into account when computing EFMs [31,32,33]. In all generality, such a constraint may be given by an arbitrary Boolean formula in terms of the reactions i, viewed as propositional symbols, i.e., the positive literal i meaning that reaction i is active (nonzero flux) and the negative literal ¬ i meaning it is inactive (zero flux). Thus, when applying this Boolean constraint, noted B c , to a flux vector v , the positive literal i is interpreted as v i 0 , i.e., i s u p p ( v ) (4), and the negative literal ¬ i as v i = 0 , i.e., i s u p p ( v ) .
The regulatory constraint RC B c for v induced by B c can thus be defined as
RC B c ( v ) = B c ( v i ) .
Note that coupled reactions as used by Flux Coupling Analysis (FCA) [34,35,36] can be easily represented by such constraints. For example, i directionally coupled to j, meaning that zero flux through i implies zero flux through j, is expressed by B c = i ¬ j , and i partially coupled to j, meaning that zero flux through i is equivalent to zero flux through j, is expressed by B c = ( i j ) ( ¬ i ¬ j ) .
By rewriting the Boolean constraint B c in DNF (Disjunctive Normal Form), B c = k D k , the set S o l RC B c = { v R r RC B c ( v ) } of vectors v satisfying the constraint RC B c is a union of the solution spaces for each disjunct D k (and this union can be assumed to be disjoint by taking the disjuncts D k two by two inconsistent). Now a disjunct is a conjunction of literals, where a negative literal ¬ i corresponds to the constraint v i = 0 and a positive literal i to the constraint v i 0 , which can be rewritten as the disjunctive constraint ( v i < 0 ) ( v i > 0 ) . Finally, a propositional symbol i that does not appear in the disjunct corresponds to an absence of constraint on v i , which can be rewritten as the disjunction ( v i < 0 ) ( v i = 0 ) ( v i > 0 ) . Thus, the solution space for a disjunct is itself the disjoint union of subspaces each one defined by constraints of type v i o p i 0 for all i, with o p i { < , = , > } , that is, defined by { v R r s i g n ( v ) = rs } for a given sign vector rs { , 0 , + } r , i.e., an open orthant O ˚ rs (which, for rs 0 , is topologically open in the vector subspace it spans and is the interior in this subspace of the closed orthant O rs = { v R r s i g n ( v ) rs } , which is the closure of O ˚ rs in R r ). In summary, S o l RC B c is thus the disjoint union of such open orthants. Now, instead of keeping this partition of S o l RC B c in open orthants, it can be more practical to generalize this concept and deal with what we will call semi-open orthants O , i.e., orthants O without some of their faces of lesser dimension (open orthant is thus a particular case of semi-open orthant, without any facet thus without any face of lesser dimension). To do this, we can group together with any given open orthant O ˚ rs j in S o l RC B c , in a same cluster C , all other open orthants O ˚ rs k in S o l RC B c such that O rs k O rs j , i.e., O rs k is a face of O rs j , and, for any arbitrary orthant O rs l with O rs k O rs l O rs j (which is equivalent to rs k < rs l < rs j ), then O ˚ rs l C . In this case O ˚ C O ˚ is a semi-open orthant O . We can iteratively process like this by considering each time the open orthants in S o l RC B c that have not yet been selected in any cluster, which guarantees that the semi-open orthants built in this way are disjoint. In addition, by choosing as open orthant, O ˚ rs j , to start with at each iteration, a maximal one among those that remain, we are sure that no two of the semi-open orthants thus built can be grouped together to constitute a bigger semi-open orthant, i.e., the collection obtained is minimal in this sense. We finally obtain that S o l RC B c can be written as a disjoint union of semi-open orthants, with no merging possible between any two of them. However, note that such a decomposition is not unique and that an inclusion can still exist between the closures of two such semi-open orthants. If we consider the particular case of a Boolean constraint which is an arbitrary disjunct D, then for any closed r-orthant O η , i.e., for any full support sign vector η , S o l RC D η is a semi-open orthant, which is the face of O η defined by v i = 0 for all i such that ¬ i is a literal of D, without its facets of equation v i = 0 for all i such that i is a literal of D.
Lemma 2.
The set S o l RC B c of vectors in R r satisfying the regulatory constraint RC B c is a disjoint union of open orthants O ˚ . These can be grouped together such that S o l RC B c is a disjoint union of semi-open orthants O , i.e., orthants O without some of their faces, such that no two of them can be grouped together to constitute a bigger semi-open orthant. If B c is a conjunction of literals and η an arbitrary full support sign vector, then S o l RC B c η is itself a semi-open orthant.

2.2. Characterizing the Solution Space

We start from a flux cone (6) (or more generally a flux polyhedron (25)) representing the flux vectors of a metabolic network in a steady state, satisfying the stoichiometric equations, the inequalities regarding irreversibility of reactions (and possibly some inhomogeneous linear constraints). Then, we consider additional biological constraints, such as those described in Section 2.1. In all generality, these constraints will be noted C x ( v ) , where v represents a flux vector in R r and x a vector of biochemical quantities involved in the constraint (typically, metabolite concentrations and enzyme concentrations). Some parameters (such as stoichiometric coefficients), that are not made explicit in the notation for sake of simplicity, are also present. In the following, the solution space in R r of an arbitrary constraint C will be noted S o l C , defined as S o l C = { v R r C ( v ) } . For given values of quantities x , the solution space is thus the constrained flux cone subset (not necessarily a cone, depending on C x ):
S o l S , Irr , C ( x ) = C F C C ( x ) = { v R r S v = 0 , v i 0 for i Irr , C x ( v ) } = S o l S , Irr S o l C x = F C S o l C x
or, more generally, the constrained flux polyhedron subset (not necessarily a polyhedron, depending on C x ):
S o l S , Irr , ILC , C ( x ) = C F P C ( x ) = { v R r S v = 0 , v i 0 for i Irr , G v h , C x ( v ) } = S o l S , Irr , ILC S o l C x = F P S o l C x = C F C C ( x ) P ILC .
Now, very often, most of the values of quantities x , if not all, are unknown. We then only require consistency of the set of constraints C x , i.e., the existence of values for variables x such that the constraints are satisfied. This means that we replace the constraint C x by the constraint x C x whose solution space is S o l x C x = { v R r { x C x ( v ) } } = x S o l C x . The solution space becomes
S o l S , Irr , C = C F C C = { v R r S v = 0 , v i 0 for i Irr , x C x ( v ) } = S o l S , Irr S o l x C x = x S o l S , Irr , C ( x ) = F C x S o l C x
or
S o l S , Irr , ILC , C = C F P C = { v R r S v = 0 , v i 0 for i Irr , G v h , x C x ( v ) } = S o l S , Irr , ILC S o l x C x = x S o l S , Irr , ILC , C ( x ) = F P x S o l C x = C F C C P ILC .
If some but not all of the values of quantities x are known, x concerns only those x whose value is unknown, the known values being integrated in C as parameters to simplify the notations.
We will obtain examples of such constraints x C x ( v ) from the cases of thermodynamic constraint TC and kinetic constraint KC .
Depending on the constraint considered, the solution space is in general no longer a cone or polyhedron and no longer convex. Nevertheless, the definitions of extreme vectors and extreme points, of elementary flux vectors and elementary flux points, and of elementary flux modes are still valid, as, respectively, non-decomposable and convex non-decomposable vectors, conformally non-decomposable and convex-conformally non-decomposable vectors, and support-minimal vectors, where decomposition and support minimality have to be understood w.r.t. vectors of the solution space. It is then clear that if a vector of F C (resp., F P ) satisfies one of those properties in F C (resp., F P ) and satisfies the given constraint, then it satisfies the same property in the solution space, i.e.,
E x V s ( F C ) S o l C x E x V s ( C F C C ( x ) ) E x V s ( F C ) S o l x C x E x V s ( C F C C ) E x V s ( C F P ) S o l C x E x V s ( C F P C ( x ) ) E x V s ( C F P ) S o l x C x E x V s ( C F P C ) E x P s ( F P ) S o l C x E x P s ( C F P C ( x ) ) E x P s ( F P ) S o l x C x E x P s ( C F P C ) E F V s ( F C ) S o l C x E F V s ( C F C C ( x ) ) E F V s ( F C ) S o l x C x E F V s ( C F C C ) E F V s ( C F P ) S o l C x E F V s ( C F P C ( x ) ) E F V s ( C F P ) S o l x C x E F V s ( C F P C ) E F P s ( F P ) S o l C x E F P s ( C F P C ( x ) ) E F P s ( F P ) S o l x C x E F P s ( C F P C ) E F M s ( F C ) S o l C x E F M s ( C F C C ( x ) ) E F M s ( F C ) S o l x C x E F M s ( C F C C ) E F M s ( F P ) S o l C x E F M s ( C F P C ( x ) ) E F M s ( F P ) S o l x C x E F M s ( C F P C )
For EFMs, the above formulas are valid whatever the structure of S o l C . It is also the case for ExPs and EFPs but, in practice, only really meaningful if C F P C ( x ) (resp., C F P C ) is a convex polyhedron. Finally, for ExVs and EFVs, this is only meaningful if C F C C ( x ) (resp., C F C C ) is a convex polyhedral cone and C F P C ( x ) (resp., C F P C ) is a convex polyhedron. This is the case when S o l C x (resp., S o l x C x ) is itself a convex polyhedral cone and we will see that, for almost all common biological constraints, this solution space is actually a union of such cones or of semi-open cones and thus the formulas will apply for each conical component. In this case we have C C F P C = C F C C C P ILC = C F P S o l C . This can be generalized when S o l C is a convex polyhedron, and thus C F C C and C F P C too, with C C F C C = F C C S o l C and C C F P C = C C F C C C P ILC = C F P C S o l C , by replacing S o l C by C S o l C in the formulas above regarding ExVs and EFVs.
However, in all the above cases, we generally do not have the reciprocal subset inclusion. This is what we will study for particular constraints.

2.2.1. Application to Thermodynamics

Assume first that the concentrations of the metabolites (resp., of the internal metabolites) are known and given. From Equations (51) and (52) and Lemma 1, we obtain
C F C TC ( M ) = F C ts M C F C TC ̲ ( M ̲ ) = F C ts M ̲
C F P TC ( M ) = F P ts M C F P TC ̲ ( M ̲ ) = F P ts M ̲ .
C F C TC ( M ) , C F C TC ̲ ( M ̲ ) are thus flux cones and C F P TC ( M ) , C F P TC ̲ ( M ̲ ) flux polyhedra (possibly equal to { 0 } or empty (for polyhedra) if the flux directions imposed by the concentrations of the metabolites and by the second law of thermodynamics are incompatible with the steady-state assumption). When nonempty, they consist of a single flux tope. From Equation (18), (36) and (39), it follows that
E F V s ( C F C TC ( M ) ) = E F V s ( F C ) ts M = E F M s ( C F C TC ( M ) ) = E F M s ( F C ) ts M = E x V s ( F C ts M )
E F P s ( C F P TC ( M ) ) = E F P s ( F P ) ts M = E x P s ( F P ts M ) E F V s ( C F P TC ( M ) ) = E F V s ( C F P ) ts M = E x V s ( C F P ts M ) E F M s ( C F P TC ( M ) ) = E F M s ( F P ) ts M
and the same for C F C TC ̲ ( M ̲ ) and C F P TC ̲ ( M ̲ ) . Now, considering the decomposition of F C (resp., F P ) into flux topes, we can proceed flux tope by flux tope for the computation of the sets above. Starting from a FT F C η (resp., F P η ), we have ( F C η ) ts M = F C η (resp., ( F P η ) ts M = F P η ) with η = inf ( η , ts M ) (where inf is defined component-wise with inf ( , + ) = 0 ) and the four sets above are equal to ExVs ( F C η ), ExPs ( F P η ), ExVs ( C F P η ) and EFMs ( F P η ) respectively. This is in particular the case if we split each reversible reaction into two irreversible ones as, in higher dimension, F C ˜ (resp., F P ˜ ) is reduced to a single flux tope defined by η = + , the all-positive sign vector, and thus η is obtained from ts M by changing each − into 0. Note also that F C (resp., F P ) consistent does not imply that C F C TC ( M ) (resp., C F P TC ( M ) ) is consistent ( η being full support does not imply that η is full support).
Proposition 1.
Given metabolite concentrations M , the space of flux vectors in F C (resp., F P ) satisfying the thermodynamic constraint TC M is a flux cone (resp., a flux polyhedron), made up of vectors of F C (resp., F P ) which conform to the thermodynamic sign vector ts M given in Equation (44). Its elementary flux vectors, identical to elementary flux modes, (resp., elementary flux points, elementary flux vectors and elementary flux modes) are exactly those of F C (resp., FP) that satisfy the constraint, i.e., that conform to ts M . The same result holds for internal metabolite concentrations M ̲ and thermodynamic constraint TC ̲ M ̲ with ts M ̲ .
Let us now consider the more usual case where the concentrations of the metabolites (at least those of internal metabolites) are unknown. From Equation (45), we obtain one or the other form for the thermodynamic constraint (existentially quantified on metabolite concentrations):
TC ( v ) = M TC M ( v ) = M s i g n ( v ) ts M TC ̲ ( v ) = M ̲ TC ̲ M ̲ ( v ) = M ̲ s i g n ( v ) ts M ̲ .
Though the metabolite concentrations M j are unknown, some lower bounds M j and upper bounds M j + on these concentrations are often known. They are thus added as additional constraints:
TC b ( v ) = M ( s i g n ( v ) ts M M M M + ) TC ̲ b ( v ) = M ̲ ( s i g n ( v ) ts M ̲ M ̲ M ̲ M ̲ + ) .
The solution space in R r of these constraints is thus S o l TC = M S o l TC M (idem for TC ̲ ) and S o l TC b = M M M M + S o l TC M (idem for TC ̲ b ). From the result of Lemma 1, it follows that S o l TC = M O ts M and S o l TC b = M M M M + O ts M (idem for TC ̲ and TC ̲ b with O ts M ̲ ). This is obviously enough to take the union on the maximal ts M ’s (resp., ts M ̲ ’s) when M varies. As these are at most 2 r , the union is finite.
Lemma 3.
The set S o l TC (resp., S o l TC b ) of vectors in R r satisfying the thermodynamic constraint TC (resp., TC b ) is a union of closed orthants. More precisely, it is the (finite) union for all M (resp., all bounded M ) of the sets of vectors that conform to ts M , i.e., O ts M ’s. The same result holds for TC ̲ and TC ̲ b with M and O ts M ̲ .
It follows from Equations (53), (54), (56) and (57) that
C F C TC = M F C ts M C F C TC b = M M M M + F C ts M
C F P TC = M F P ts M C F P TC b = M M M M + F P ts M
and the same for TC ̲ and TC ̲ b with M and ts M ̲ , where all unions are finite.
From this and Equations (18), (36) and (39), we get
E F V s ( C F C TC ) = M E F V s ( F C ts M ) = M E F V s ( F C ) ts M = E F M s ( C F C TC ) = M E F M s ( F C ts M ) = M E F M s ( F C ) ts M = M E x V s ( F C ts M )
E F P s ( C F P TC ) = M E F P s ( F P ts M ) = M E F P s ( F P ) ts M = M E x P s ( F P ts M ) E F V s ( C F P TC ) = M E F V s ( C F P ts M ) = M E F V s ( C F P ) ts M = M E x V s ( C F P ts M ) E F M s ( C F P TC ) = M E F M s ( F P ts M ) = M E F M s ( F P ) ts M
and the analog for TC b , TC ̲ and TC ̲ b . Now, considering the decomposition of F C (resp., F P ) into flux topes, we can proceed flux tope by flux tope for the computation of the sets above. Starting from a FT F C η (resp., F P η ), we have M ( F C η ) ts M = i F C η i (resp., M ( F P η ) ts M = i F P η i ) where { η i } are the maximal sign vectors in {inf ( η , ts M ) M R + * ( m + m ¯ ) } and the four sets above are equal to i ExVs ( F C η i ), i ExPs ( F P η i ), i ExVs ( C F P η i ) and i EFMs ( F P η i ), respectively. This is in particular the case if we split each reversible reaction into two irreversible ones as, in higher dimension, F C ˜ (resp., F P ˜ ) is reduced to a single flux tope defined by η = + . The analog holds for TC b , TC ̲ and TC ̲ b .
Proposition 2.
The space of flux vectors in F C (resp., F P ) satisfying the thermodynamic constraint TC , or TC b , is a finite union of flux cones (resp., flux polyhedra), obtained as those vectors of F C (resp., F P ) which conform to ts M , for a certain M , or bounded M . It is thus no longer convex but made up of particular faces F C η i (resp., F P η i ) of each flux tope F C η for F C (resp., F P η for F P ), with η i η . Its elementary flux vectors, identical to elementary flux modes, (resp., elementary flux points, elementary flux vectors and elementary flux modes) are exactly those of F C (resp., FP) that satisfy the constraint, i.e., that conform to a certain ts M , and coincide thus with the extreme vectors (resp., the extreme points, extreme vectors and elementary modes) of the F C η i ’s (resp., F P η i ’s). The same result holds for TC ̲ , or TC ̲ b , with M and ts M ̲ .
We will now refine this result in order to characterize the faces involved. Consider the case of a flux cone where the concentrations of external metabolites are given and assume first there are no bounds on the concentrations of the internal metabolites [37]. As C F C TC ̲ = F C M ̲ O ts M ̲ and C F P TC ̲ = F P M ̲ O ts M ̲ , let us begin by studying M ̲ O ts M ̲ . We have x M ̲ O ts M ̲ M ̲ R + * m ( s i g n ( x ) ts M ̲ ) M ̲ R + * m i , 1 i r , ( s i g n ( x i ) s i g n ( j M ̲ j S j i K ^ e q i ) ) from Equation (44). By applying the monotonic logarithm function and noting L M ̲ R m the vector whose components are given by L M ̲ j = l n ( M ̲ j ) , we obtain x M ̲ O ts M ̲ L M ̲ R m i , 1 i r , ( s i g n ( x i ) s i g n ( j S j i L M ̲ j l n ( K ^ e q i ) ) ) L M ̲ R m i s u p p ( x ) ( s i g n ( x i ) j S j i L M ̲ j < s i g n ( x i ) l n ( K ^ e q i ) ) . We will now make use of Gale’s theorem (or Kuhn–Fourier theorem), which is a form of Farkas’ duality lemma (for two vectors x , y , x < y means x i < y i for all i):
Gale’s theorem (a form of Farkas’ lemma). For any A R p × q and b R p , exactly one of the following statements holds:
(a) 
there exists y R q such that A y < b ;
(b) 
there exists z R p \ { 0 } such that z 0 , z T A = 0 and z T b 0 .
Applying this theorem in the formula above to A R s u p p ( x ) × m (where we note R s u p p ( x ) = R r i x i = 0 { z i = 0 } the subspace of vectors having a null component outside the support of x ), given by A i j = s i g n ( x i ) S j i , and b R s u p p ( x ) , given by b i = s i g n ( x i ) l n ( K ^ e q i ) , provides x M ̲ O ts M ̲ z R s u p p ( x ) \ { 0 } ( z 0 , j , 1 j m , i S j i s i g n ( x i ) z i = 0 i s i g n ( x i ) z i l n ( K ^ e q i ) > 0 ) z R r \ { 0 } ( s i g n ( z ) s i g n ( x ) , S z = 0 z T ln K ^ eq > 0 ) where ln K ^ eq R r is defined by: ( ln K ^ eq ) i = l n ( K ^ e q i ) .
We thus obtain C F C TC ̲ = { v F C z F C \ { 0 } ( s i g n ( z ) s i g n ( v ) z T ln K ^ eq > 0 ) } . This can be rewritten as C F C TC ̲ = { v F C F C s i g n ( v ) \ { 0 } H ln K ^ eq + } , where H ln K ^ eq + = { z R r z T ln K ^ eq > 0 } is an open half-space with the hyperplane H ln K ^ eq = { z R r z T ln K ^ eq = 0 } as a frontier. Considering the decomposition of F C into FTs F C η , this is equivalent to C F C TC ̲ η = { v F C η F C s i g n ( v ) \ { 0 } H ln K ^ eq + } for all η ’s maximal sign vectors in sign( F C ). Now, note that, in this formula, F C s i g n ( v is the minimal face (for set inclusion) of F C η containing v . We finally obtain finally:
C F C TC ̲ η = i F C η i for all maximal η i η s . t . F C η i \ { 0 } H ln K ^ eq + .
That is to say, C F C TC ̲ is the union of all the maximal faces F C η i of the FTs for F C such that F C η i \ { 0 } H ln K ^ eq + or, equivalently, the union of all maximal faces F i of the F C O ’s for all r-orthants O such that F i \ { 0 } H ln K ^ eq + (see Figure 3).
A topological consequence of this characterization is that C F C TC ̲ \ { 0 } is a connected set (actually arc-connected). Another consequence, regarding EFMs (or, equivalently, EFVs), is:
E F M s ( C F C TC ̲ ) = E F M s ( F C ) S o l TC ̲ = E F M s ( F C ) H ln K ^ eq +
i.e., the elementary flux modes in the (non-convex) cone of those flux vectors in F C satisfying the thermodynamic constraint TC ̲ are exactly those elementary flux modes in F C that satisfy TC ̲ , i.e., that belong to the open half-space H ln K ^ eq + . We can equivalently write
E F M s ( C F C TC ̲ ) = E F M s ( F C ln K ^ eq + ) \ H ln K ^ eq with F C ln K ^ eq + = { v R r S v = 0 , v i 0 for i Irr , v T ln K ^ eq 0 }
where F C ln K ^ eq + is the polyhedral cone obtained from F C merely by adding the homogeneous linear inequality v T ln K ^ eq 0 [38]. Each flux vector of C F C TC ̲ is a conformal conical sum of these EFMs, but the converse is false as C F C TC ̲ is not convex. Thus the set of EFMs, considered globally, does not characterize the solution space C F C TC ̲ . More precisely, we have
C F C TC ̲ η = i c o n e ( E i ) with E i = E F M s ( F C η i )
where F C η i is as in Equation (66), i.e., C F C TC ̲ is characterized by the decomposition of the set of EFMs into the (non-disjoint) subsets E i . Now, the E i ’s are exactly the maximal subsets of EFMs included in a given flux tope for F C (i.e., in a given r-orthant O) and whose conical hull is included in C F C TC ̲ , i.e., all vectors in this hull must satisfy the constraint TC ̲ : E i maximal such that E i EFMs ( C F C TC ̲ ) and E i O with O r-orthant and c o n e ( E i ) C F C TC ̲ . Such an E i is called a largest thermodynamically consistent set (LTCS) of EFMs [39] in O (or, equivalently, in the flux tope F C η = F C O defined by O).
We could want to estimate the ratio of thermodynamically feasible EFMs on all EFMs, i.e., the ratio of EFMs ( C F C TC ̲ ) on EFMs ( F C ) . If all reactions are reversible ( r I = 0 ), then the function v v maps the EFMs on one side of hyperplane H ln K ^ eq onto the EFMs on the other side. Thus, if we neglect the EFMs that might belong to this hyperplane, it means that 50% of the EFMs are thermodynamically feasible (and thus only 50% eliminated). Now, intuitively, irreversible reactions given in Irr come from an expert knowledge that can be seen as a form of compiled thermodynamic knowledge, as we saw that it is thermodynamics which imposes the direction in which a reaction may proceed. Therefore, we can assume, provided the adequacy of the model for some given environment (such as the concentrations of external metabolites, supposed here to be known), that any thermodynamically feasible EFM satisfies these irreversibility constraints. Under this assumption, the irreversibility constraints rule out only thermodynamically unfeasible EFMs so if we then apply the thermodynamic constraint TC ̲ we only eliminate at most than 50% of the remaining EFMs (from 50% when all reactions are reversible to 0% when all are irreversible, without splitting any reversible reaction into two irreversible ones).
If inhomogeneous linear constraints are added, we have C F P TC ̲ = C F C TC ̲ P ILC , with P ILC = { v R r G v h } Equation (25) and we obtain in the same way, for all η ’s maximal sign vectors in sign( F C )
C F P TC ̲ η = i F P η i for all maximal η i η s . t . F C η i \ { 0 } H ln K ^ eq + .
That is to say, C F P TC ̲ is the union of all the F P η i = F C η i P ILC such that F C η i is a maximal face of a FT for F C verifying F C η i \ { 0 } H ln K ^ eq + or, equivalently, the union of the F i P ILC for all maximal faces F i of the F C O ’s for all r-orthants O such that F i \ { 0 } H ln K ^ eq + . Take care because the F P η i ’s involved are faces of F P η included in H ln K ^ eq + { 0 } , but not any face of F P η included in H ln K ^ eq + { 0 } is included in C F P TC ̲ (as such a face may be defined by inhomogeneous equality constraints coming from ILC and not only by nullity constraints of the form v i = 0 ). C F P TC ̲ is no longer a connected set in general.
We can sum up these results as follows.
Proposition 3.
The space of flux vectors in F C satisfying the thermodynamic constraint TC ̲ is a finite union of flux cones, obtained as all the maximal faces of all the flux topes F C η that are entirely contained (except the null vector) in the open half-space H ln K ^ eq + = { z z T ln K ^ eq > 0 } . The thermodynamically feasible EFMs are thus those EFMs which belong to this half-space and can be simply computed as the EFMs of the flux cone obtained from F C by adding to it the homogeneous linear inequality v T ln K ^ eq 0 (and removing those EFMs that would belong to the frontier hyperplane of H ln K ^ eq + ). The set of these EFMs can be decomposed into (non-disjoint) maximal subsets of EFMs belonging to a same flux tope (i.e., a same r-orthant) and whose conical hull is made up of thermodynamically feasible vectors, each of these subsets thus representing the set of EFMs of one of the maximal faces above. At most 50% of the EFMs can thus be ruled out as thermodynamically infeasible, if we assume that no thermodynamically feasible flux vector violates given irreversibility constraints. In the presence of additional inhomogeneous linear constraints on flux vectors given by G v h , the space of flux vectors in FP satisfying TC is a finite union of flux polyhedra, obtained as intersections of the flux cones above with the polyhedron defined by G v h . All these results hold for constraint TC by just replacing the K ^ e q i ’s by the K e q i ’s.
This applies in particular to F C ˜ and F P ˜ (after splitting each reversible reaction into two irreversible ones) with the simplification, as F C ˜ (resp., F P ˜ ) is reduced to a single flux tope defined by η = + , that we must just consider the maximal faces of F C ˜ that are entirely contained (except the null vector) in the open half-space H ln K ^ eq + .
Consider now the case where certain bounds on the concentrations of internal metabolites are known, i.e., the case of the thermodynamic constraint TC ̲ b . We have C F C TC ̲ b = F C M ̲ M ̲ M ̲ M ̲ + O ts M ̲ . From what precedes, we obtain x M ̲ M ̲ M ̲ M ̲ + O ts M ̲ L M ̲ R m ( i s u p p ( x ) ( s i g n ( x i ) j S j i L M ̲ j < s i g n ( x i ) l n ( K ^ e q i ) ) L M ̲ l n ( M ̲ + ) L M ̲ l n ( M ̲ ) ) . Applying Gale’s theorem in this formula to A I m I m R ( s u p p ( x ) + 2 m ) × m , with A i j = s i g n ( x i ) S j i , and b l n M ̲ + l n M ̲ R s u p p ( x ) + 2 m , with b i = s i g n ( x i ) l n ( K ^ e q i ) , provides x M ̲ M ̲ M ̲ M ̲ + O ts M ̲ z z ̲ z ¯ R r + 2 m ( z 0 , z ̲ 0 , z ¯ 0 , s i g n ( z ) s i g n ( x ) , S z + z ̲ z ¯ = 0 ( z T z ̲ T z ¯ T ) ln K ^ eq b > 0 ) where ln K ^ eq b R r + 2 m is defined by ln K ^ eq b = ln K ^ eq l n M ̲ + l n ( 1 / M ̲ ) = l n K ^ eq M ̲ + 1 / M ̲ . Let F C b = { z z ̲ z ¯ R r + 2 m ( S I m I m ) z z ̲ z ¯ = 0 , z i 0 for i Irr , z ̲ 0 , z ¯ 0 } = { z w + ( S z ) w + ( S z ) + R r + 2 m z i 0 for i Irr , w 0 } , where ( S z ) j + = m a x ( ( S z ) j , 0 ) and ( S z ) j = m a x ( ( S z ) j , 0 ) . F C b is a flux cone in R r + 2 m of dimension r + m and F C b ( R r × { 0 } ) = F C . We thus obtain C F C TC ̲ b = { v F C ( z T z ̲ T z ¯ T ) T F C b ( z 0 , s i g n ( z ) s i g n ( v ) ( z T z ̲ T z ¯ T ) ln K ^ eq b > 0 ) } . This is equivalent to: C F C TC ̲ b η = { v F C η ( z T z ̲ T z ¯ T ) T F C b ( z O s i g n ( v ) \ { 0 } ( z T z ̲ T z ¯ T ) ln K ^ eq b > 0 ) } = { v F C η z R r ( z O s i g n ( v ) \ { 0 } ( z T ( S z ) T ( S z ) + T ) ln K ^ eq b > 0 ) } for all flux topes F C η for F C . This can be rewritten as C F C TC ̲ b η = { v F C η F C b ( O s i g n ( v ) \ { 0 } × R + 2 m ) H ln K ^ eq b + } , where we note H ln K ^ eq b + = { Z R r + 2 m Z T ln K ^ eq b > 0 } the open half-space with hyperplane H ln K ^ eq b = { Z R r + 2 m Z T ln K ^ eq b = 0 } as a frontier. We have H ln K ^ eq b + R r = H ln K ^ eq + and H ln K ^ eq b R r = H ln K ^ eq . Now, note that, in the formula above, F C b ( O s i g n ( v ) × R + 2 m ) is a face of the cone F C η + b (we note η + the sign vector of dimension r + 2 m obtained by concatenation of η of dimension r and + of dimension 2 m , thus F C η + b = F C b ( O η × R + 2 m ) ) whose intersection with R r × { 0 } is equal to F C s i g n ( v ) , the minimal face of the tope F C η containing v . Actually, among all the faces of F C η + b whose intersection with R r × { 0 } is equal to F C s i g n ( v ) , this is the maximal one (without any constraint on the R + 2 m factor). Thus, finally,
C F C TC ̲ b η = j F C η j for all maximal η j η s . t . F j b \ ( { 0 } × R 2 m ) H ln K ^ eq b + ,
where F j b is the maximal face of F C η + b whose intersection with R r × { 0 } is equal to F C η j .
Note that any F C η j given by Equation (71) is a face of a certain F C η i given by Equation (66) (i.e., η j η i η j η i ). C F C TC ̲ b \ { 0 } is no longer a connected set in general. A consequence, for EFMs (or, equivalently, EFVs), is:
E F M s ( C F C TC ̲ b η ) = E F M s ( F C η ) S o l TC ̲ b = { v E F M s ( F C η ) F v b \ ( { 0 } × R 2 m ) H ln K ^ eq b + } ,
where F v b is the maximal face of F C η + b s.t. F v b ( R r × { 0 } ) = R + v . This means that the elementary flux modes in the (non-convex) cone of those flux vectors in F C η satisfying the thermodynamic constraint TC ̲ b are exactly those elementary flux modes v in F C η (or in C F C TC ̲ η ) that satisfy TC ̲ b , i.e., such that the maximal face of F C η + b whose intersection with R r × { 0 } is equal to the ray R + v is included in H ln K ^ eq + ( { 0 } × R 2 m ) . Note that EFMs ( C F C TC ̲ b ) is thus given by
{ v E F M s ( F C ) z R r ( s i g n ( z ) s i g n ( v ) , ( z T ( S z ) T ( S z ) + T ) ln K ^ eq b 0 z = 0 ) } .
Each flux vector of C F C TC ̲ b is a conformal conical sum of these EFMs, but the converse is false as C F C TC ̲ b is not convex. More precisely, we have
C F C TC ̲ b η = j c o n e ( E j ) with E j = E F M s ( F C η j )
where F C η j is as in Equation (71), i.e., C F C TC ̲ b is characterized by the decomposition of the set of EFMs into the (non-disjoint) subsets E j . Now, the E j ’s are exactly the maximal subsets of EFMs included in a given flux tope for F C (i.e., in a given r-orthant O) and whose conical hull is included in C F C TC ̲ b , i.e., all vectors in this hull must satisfy the constraint TC ̲ b : E j maximal such that E j EFMs ( C F C TC ̲ b ) and E j O with O r-orthant and c o n e ( E j ) C F C TC ̲ b . Such an E j is called a largest thermodynamically consistent (for bounded concentrations of internal metabolites) set (LTCbS) of EFMs [39] in O (or, equivalently, in the flux tope F C η = F C O defined by O) and is included in a certain LTCS E i as in Equation (69).
If inhomogeneous linear constraints ILC are added, we have, for all η ’s maximal sign vectors in sign( F C ):
C F P TC ̲ b η = j F P η j
with η j as in Equation (71). That is to say, C F P TC ̲ b is the union of all the F P η j = F C η j P ILC such that F C η j is given as in Equation (71).
We can sum up these results as follows.
Proposition 4.
For F C a flux cone in R r , let us define its lift to R r + 2 m as the flux cone F C b = { Z R r × R + 2 m ) ( S I m I m ) Z = 0 , Z i 0 for i Irr , 1 i r } . Thus, F C b ( R r × { 0 } ) = F C . For any flux tope F C η for F C and any face F C of F C η , its lift to R r + 2 m is defined as the maximal face of F C η + b whose intersection with R r × { 0 } is equal to F C . The space of flux vectors in F C satisfying the thermodynamic constraint TC ̲ b is a finite union of flux cones, obtained as all the maximal faces of all the flux topes F C η whose lifts to R r + 2 m are entirely contained (except vectors from { 0 } × R 2 m ) in the open half-space H ln K ^ eq b + = { Z R r + 2 m Z T ln K ^ eq b > 0 } , where ln K ^ eq b = l n ( K ^ eq M ̲ + 1 / M ̲ ) T . The thermodynamically feasible EFMs in F C η are those EFMs in F C η whose lifts (as rays) are contained (except vectors from { 0 } × R 2 m ) in this half-space. The set of these EFMs can be decomposed into (non-disjoint) maximal subsets of E F M s belonging to a same flux tope (i.e., a same r-orthant) and whose conical hull is made up of thermodynamically feasible vectors, each of these subsets representing thus the set of EFMs of one of the maximal faces above. In the presence of additional inhomogeneous linear constraints on flux vectors given by G v h , the space of flux vectors in FP satisfying TCb is a finite union of flux polyhedra, obtained as intersections of the flux cones above with the polyhedron defined by G v h .

2.2.2. Application to Kinetics

In the same way, we obtain for the kinetic constraint in the absence of knowledge regarding values of concentrations of enzymes and metabolites:
KC ( v ) = E , M KC E , M ( v ) = E , M v = E κ ( M ) .
Once more we can add optional lower and upper bounds M j and M j + on metabolite concentrations and/or lower and upper bounds E i and E i + on enzyme concentrations if they are known, and also a global enzymatic resource constraint, which is often considered, as c T E W , where c is a constant positive vector of size r and W a positive constant:
KC b ( v ) = E , M ( v = E κ ( M ) c T E W E E E + M M M + ) .
We can also consider intermediate constraints, existentially quantified only on metabolite concentrations if enzyme concentrations are known:
KC E ( v ) = M KC E , M ( v )
possibly with bounds on metabolite concentrations in the quantification:
KC E b ( v ) = M ( KC E , M ( v ) M M M + )
or only on enzyme concentrations if metabolite concentrations are known:
KC M ( v ) = E KC E , M ( v )
possibly with enzymatic resource constraint and bounds on enzyme concentrations in the quantification:
KC M b ( v ) = E ( KC E , M ( v ) c T E W E E E + ) .
The solution space in R r of the constraint KC M is S o l KC M = E S o l KC E , M = E R + r { E κ ( M ) } = { v R r s i g n ( v ) s i g n ( κ ( M ) ) } = { v R r s i g n ( v ) ts M } = S o l TC M , as the sign of κ ( M ) is the thermodynamic sign vector ts M . It is thus the same as the solution space of the thermodynamic constraint TC M and is the closed orthant O ts M . It follows that the solution space of the constraint KC , given by S o l KC = M S o l KC M = M S o l TC M = S o l TC , is the same as the solution space of the thermodynamic constraint TC and is the finite union of the closed orthants O ts M . Similarly, S o l KC b = S o l TC b if the bounds are only on the metabolite concentrations, i.e., there are no bounds on enzyme concentrations. From S o l KC M b = E E E + c T E W { E κ ( M ) } , it follows that the flux vectors in S o l KC M b are defined by the linear inequalities: κ i ( M ) E i v i κ i ( M ) E i + for those i’s such that κ i ( M ) > 0 , κ i ( M ) E i + v i κ i ( M ) E i for those i’s such that κ i ( M ) < 0 , v i = 0 for those i’s such that κ i ( M ) = 0 and i κ i ( M ) 0 c i κ i ( M ) 1 v i W i κ i ( M ) = 0 c i E i . Thus S o l KC M b is a convex polyhedron (a parallelepiped contained in the closed orthant O ts M , truncated by a hyperplane). Now, if E = 0 , i.e., in the absence of positive lower bounds on enzyme concentrations, this polyhedron has 0 as a vertex and is nothing else than S o l KC M , i.e., the closed orthant O ts M , truncated as a parallelepiped by the faces given by v i = κ i ( M ) E i + for those i’s such that κ i ( M ) 0 and by the hyperplane i κ i ( M ) 0 c i κ i ( M ) 1 v i = W . Consequently, S o l KC b = M S o l KC M b is equal to a certain truncation of S o l KC = S o l TC , defined as the finite union of truncations of the closed orthants O ts M (each O ts becoming a parallelepiped, after being truncated by hyperplanes according to equations v i = s u p ( κ i ( M ) ) E i + if ts i = + (resp., v i = i n f ( κ i ( M ) ) E i + if ts i = ), where s u p (resp., i n f ) applies to those M such that s i g n ( κ ( M ) ) = ts , which gives a polyhedron, but also, in the presence of an enzymatic resource constraint, by an algebraic, nonlinear, surface, which gives in this case a local solution space in O ts that is no longer a polyhedron and is not necessarily convex).
Proposition 5.
Given metabolite concentrations M , the kinetic constraint KC M is identical to the thermodynamic constraint TC M and thus the set S o l KC M of vectors in R r satisfying KC M is the closed orthant defined by the thermodynamic sign vector ts M . The kinetic constraint KC is identical to the thermodynamic constraint TC and thus the set S o l KC of vectors in R r satisfying KC is a union of closed orthants. In the presence of bounds only on metabolite concentrations (and not on enzyme concentrations), the kinetic constraint KC b is identical to the thermodynamic constraint TC b and thus the set S o l KC b of vectors in R r satisfying KC b is a union of closed orthants. Therefore, these kinetic constraints boil down to thermodynamic constraints and the results regarding the geometrical structure of the corresponding spaces of flux vectors and the characterization of elementary flux vectors (or elementary flux modes) given by Propositions 1–4 apply: in particular, C F C KC ( M ) = C F C TC ( M ) is a flux cone and C F C KC = C F C TC and C F C KC b = C F C TC b are finite unions of flux cones. For a given M and in the presence of bounds on enzyme concentrations, the set of vectors in R r satisfying KC M b is a convex polyhedron and C F C KC b ( M ) is thus a flux polyhedron. In the particular case where positive lower bounds on enzyme concentrations are absent, C F C KC b ( M ) is just the parallelepiped obtained by truncating the flux cone C F C KC ( M ) = C F C TC ( M ) by hyperplanes originating from the upper bounds on enzyme concentrations and by a hyperplane originating from the enzymatic resource constraint and coincides thus with the said flux cone in a certain adequate neighborhood of 0 , i.e., for sufficiently small values of the fluxes. In this case, C F C KC b is thus a truncation (by an algebraic surface) of C F C TC and coincides with this union of flux cones in a certain adequate neighborhood of 0 and the characterization of elementary flux vectors remains valid in this neighborhood. These results extend in the presence of additional inhomogeneous linear constraints on flux vectors given by G v h by intersecting the solution spaces above with the polyhedron defined by G v h , giving rise to flux polyhedra (results in a neighborhood of 0 holding only if h < 0).
However, the geometric structure of S o l KC E , of S o l KC E b and of S o l KC b in the presence of positive lower bounds on enzyme concentrations, depends on the kinetic function κ ( M ) and C F C KC ( E ) , C F C KC b ( E ) and C F C KC b are generally neither polyhedra nor convex.
Proposition 5 has important consequences on constrained enzyme allocation problems in kinetic metabolic networks. Considering a kinetic metabolic network, with possible bounds on metabolite concentrations, but not on enzyme concentrations, i.e., with solution space C F C KC = C F C TC , or C F C KC b = C F C TC b , which is a finite union, for M varying, of the flux cones F C ts M , the generic enzyme allocation problem consists in maximizing the specific flux (or specific rate, i.e., rate expressed per unit of biomass amount) of a given reaction, say k, defined as v k / E T , where v is a flux vector in C F C KC or C F C KC b and E T denotes the total protein content in the system. E T is expressed in all generality as a fixed weighted sum of the enzyme concentrations, E T = i = 1 r w i E i (the w i ’s being given positive weights that denote the fraction of the resource used per unit of enzyme), able to encode different enzymatic constraints (such as limited cellular or membrane surface space) as well as other constraints regarding the abundance of certain enzymes. Likewise, the steady-state flux component v k > 0 may stand for diverse metabolic processes, ranging from the synthesis rate of a particular product within a specific pathway to the rate of overall cellular growth. The formation rate of a metabolic product expressed per gram of biomass and the specific growth rate of a cell are both examples of such specific rates. We look for maximization in the solution space by varying the metabolite concentrations (inside their bounds, if any) and the enzyme concentrations (without bounds), which gives rise to a complex non-convex optimization problem. Now, maximizing v k / E T is equivalent to fixing the rate v k to a positive value, e.g., to 1, and minimizing the E T needed to attain this level of v k . This means minimizing the function i = 1 r ( w i / κ i ( M ) ) v i by varying M (with possible bounds) and, for each given M , the flux vector v in F P M = F C ts M { v v k = 1 } (without bounds on v as we assume there are no bounds on enzyme concentrations). If all F P M ’s are empty, the problem is unsolvable, i.e., v k > 0 is incompatible with the kinetics. Otherwise, i.e., when the problem is solvable, we consider successively each nonempty F P M . Such an F P M is a flux polyhedron whose elementary flux points (which are equal to the extreme points or vertices) correspond to the intersections of the hyperplane { v v k = 1 } with extreme rays (edges) of F C ts M , i.e., EFMs of C F C KC or C F C KC b , and whose elementary vectors (equal to extreme vectors), if any, correspond to the extreme rays of F C ts M { v v k = 0 } . As, for M fixed, the function to minimize is linear in v , its minimum on F P M is reached on a lower-dimensional face of F P M (as v i / κ i ( M ) 0 , w i > 0 and F P M is included in a closed orthant, this face is necessarily a convex hull of certain extreme points of F P M even if it is not a polytope, i.e., no extreme vector may occur as one of the generators of this face), and thus reached in particular at at least one extreme point, i.e., at an EFM. Now, as the total number of EFMs is finite, so is the number of those for which the minimum of the function occurs for any given M , considering all nonempty F P M ’s. Therefore, when M varies in its domain, we obtain the result that the maximum of the specific flux v k / E T occurs (at least) at an EFM of C F C KC or C F C KC b . In the case of C F C KC b and assuming the κ i ( M ) ’s are continuous, this maximum is attained at an EFM at finite metabolite concentrations as M then varies in the compact set j [ M j , M j + ] . In the case of C F C KC , the maximum might not be attained at finite metabolite concentrations. This is the result already obtained in [40,41] and we followed a similar proof, but relying this time on a precise characterization of the solution space C F C KC or C F C KC b and of the EFMs given by the Proposition 5, which was not the case in the above-quoted works. Finally, the enzyme allocation problem can theoretically be solved by computing all the thermodynamically feasible EFMs having k in their support (i.e., satisfying the Boolean constraint B c = k , see next subsection), say { e l } (all components of which are fixed by e k l = 1 ), and, for each one, by computing the minimum (if it exists) of i s u p p ( e l ) ( w i e i l ) / κ i ( M ) for M varying in its domain, such that s i g n ( κ i ( M ) ) = s i g n ( e i l ) for all i s u p p ( e l ) , which is a much simpler optimization problem than the initial one. The global minimum, if it exists, is then the smallest of these local minima, for all e l ’s. We then obtain the maximum value of the specific flux v k / E T , an EFM where this maximum occurs and the values of the metabolite concentrations for which it occurs.
Proposition 6.
Given a kinetic metabolic network with possible bounds on metabolite concentrations, but not on enzyme concentrations, i.e., with solution space C F C KC or C F C KC b , if the enzyme allocation problem, which consists in maximizing the specific flux (rate per unit of biomass amount) in a given reaction k, i.e., v k / E T , where v is a flux vector in C F C KC or C F C KC b and E T = i = 1 r w i E i is the total protein content in the system (the w i ’s being fixed positive weights), has an optimal solution (which is always the case for C F C KC b if the problem is solvable), then this solution is reached in particular at an EFM of C F C KC or C F C KC b .

2.2.3. Application to Regulatory Constraints

From Lemma 2, we obtain that C F C RC B c (resp., C F P RC B c ) is the disjoint union of F C O ˚ ’s (resp., F P O ˚ ’s), for certain open orthants O ˚ and, after merging, the disjoint union of F C O ’s (resp., F P O ’s), for certain semi-open orthants O (orthants without some of their faces), without any further merging possible.
We will call open polyhedral cone (resp., open polyhedron) in dimension r > 0 an r-polyhedral cone C (resp., r-polyhedron P) without its facets, and we will note it C ˚ (resp., P ˚ ) as it coincides with the topological interior of C (resp., P) in the affine span of C (resp., P). Reciprocally, C (resp., P) is the topological closure of C ˚ (resp., P ˚ ) in R r . C ˚ (resp., P ˚ ) is defined as in Equation (7) (resp., in Equation (26)) but with strict inequalities (and equalities for defining its affine space). We will in particular consider open faces of a cone C or polyhedron P as F ˚ for F a face of C or P. C (resp., P) is the disjoint union of its open faces (by convention, a face of dimension 0, i.e., a vertex, is equal to its open face).
It follows from these definitions that, for any closed r-orthant O and any open orthant O ˚ O , F C O ˚ (resp., F P O ˚ ) is an open face of F C O (resp., a disjoint union of open faces of F P O , as we have to keep faces corresponding to G v = h ). In all, we obtain that C F C R C B c O (resp., C F C R C B c O ) is a disjoint union of open faces of the flux cone F C O (resp., the flux polyhedron F P O ) or, equivalently, that, for any flux tope FC for FC(resp., FP for FP), C F C R C B c η (resp., C F P R C B c η ) is a disjoint union of open faces of FC (resp., FP). As we grouped together open orthants into semi-open orthants in Lemma 2, we can also group together with such an open face F ° all those other open faces F ° in question where F is a face of F to obtain thus a (minimal) disjoint union of semi-open polyhedral cones (resp., semi-open polyhedra). Here, we call semi-open polyhedral cone C ° (resp., semi-open polyhedron P ° a polyhedral cone C (resp., polyhedron P) without certain (between zero and all) of its faces of lesser dimension, that can be thus expressed as a disjoint union of certain (between all and only C ° , resp., P ° ) of the open faces of C (resp., P). We have: C ° C ° C , resp., P ° P ° P ).
Proposition 7.
Given an arbitrary Boolean constraint Bc, the solution space C F C R C B c (resp., C F P R C B c ) for the regulatory constraint R C B c is a finite disjoint union of open polyhedral cones (resp., open polyhedra), which are certain open faces of all the flux topes F c η (resp., F P η ). Grouping together certain of these open faces according to Lemma 2, we obtain a disjoint union of semi-open (i.e., without certain of their faces of lesser dimension) faces of the F c η ’s (resp., F c η ’s), without any possible further merging of two of them to make up a bigger semi-open face.
Note that the rules presented in the proof of Lemma 2 to group together open faces into semi-open faces are applied globally to all flux topes. If we choose to apply them separately for each flux tope, then the union of semi-open faces obtained is no longer disjoint in general as two such semi-open faces for two different flux topes may have a nonempty intersection (as C F C R C B c O (resp., C F P R C B c O ) and C F C R C B c O (resp., C F P R C B c O ), for two different closed r-orthants O and O′ may have open faces in common). Note also that, even after this merging, there may exist a strict inclusion relationship between the closures of two semi-open faces.
From this geometrical structure of the solution space, we will now deduce what its EFMs and its EFVs (or EFPs) are. Let us begin with a preliminary remark. We know that, by definition of EFVs, we have, for a flux cone F C , EFVs ( F C ) = η maximal in s i g n ( F C ) EFVs ( F C η ) (and idem for a flux polyhedron F P with EFPs and EFVs), which remains true for any constrained flux cone subset C F C C (or constrained flux polyhedron subset C F P C ) whatever the biological constraint C is. We saw that this equality was also satisfied by EFMs for F C : EFMs ( F C ) = η maximal in s i g n ( F C ) EFMs ( F C η ) (and idem for F P ) and remained true for C F C C for C any thermodynamic constraint or any kinetic constraint of the form KC M , KC or KC b (in the absence of bounds on enzyme concentrations), thus allowing in these cases to decompose the identification of EFMs flux tope by flux tope, as for EFVs. However this is no longer true for regulatory constraints. Actually, counter-examples can be found in the presence of reversible reactions, used in a given direction in an EFM local to a certain flux tope and in the other direction in an EFM local to another flux tope, choosing the Boolean constraint such that it imposes a strict set inclusion between the supports of these two local EFMs.
Example 1.
Consider the simple network comprising one reaction R : A B , where A and B are the two internal metabolites, and four transfer reactions T 1 : A , T 2 : B , T 3 : A , T 4 : B , and assume that R and T 1 are reversible (we will note par rev() the reversed reaction), the three other reactions being irreversible (see Figure 4). It is obvious that e 1 = ( 1 , 1 , 0 , 0 , 1 ) T , made up of T 1 , R , T 4 , e 4 = ( 1 , 1 , 1 , 0 , 0 ) T , made up of T 2 , r e v ( R ) , r e v ( T 1 ) , and e 3 = ( 0 , 0 , 1 , 0 , 1 ) T , made up of T 2 , T 4 , are EFMs of F C . Consider now the Boolean constraint B c = T 1 T 4 . Then, e 1 is an EFM of C F C RC B c , belonging to C F C RC B c ( + , + , + , + , + ) T . Moreover, any e = λ e 3 + e 4 with λ > 0 is an EFM of C F C RC B c ( , , + , + , + ) T (obviously, the sub-pathways e 3 and e 4 of e do not belong to C F C RC B c ). As s u p p ( e 1 ) = { R , T 1 , T 4 } and s u p p ( e ) = { R , T 1 , T 2 , T 4 } , we have s u p p ( e 1 ) s u p p ( e ) and thus e is not an EFM of C F C RC B c .
Now, from a biological point of view, it is not relevant to compare supports of two pathways with a certain reaction in a given direction in the first support and in the other direction in the second support (case of R and T 1 in the example). This means that the useful concept concerning minimality is not support-minimality, but sign-minimality (exactly in the same way as, concerning decomposition, we saw that the useful concept is not non-decomposability but conform non-decomposability), which is equivalent to comparing supports separately for each closed r-orthant, i.e., for each flux tope. We will thus identify the EFMs flux tope by flux tope (note that this is analog to distinguishing a positive flux from a negative flux in a regulatory constraint, e.g., distinguishing the constraints T 1 + T 4 and T 1 T 4 in the example above, which could be done by adopting a tri-valued logic instead of a Boolean logic to represent these constraints; this is obviously done automatically when splitting each reversible reaction into two irreversible ones, where only the positive r-orthant has to be considered).
Therefore, we will consider in the following an arbitrary closed r-orthant O given by a full support sign vector η (with η i = + for i Irr ) and consider the part of the solution space limited to this orthant, i.e., C F C RC B c η (resp., C F P RC B c η ), thus we can limit ourselves to sign vectors η that are maximal in s i g n ( C F C RC B c ) (resp., s i g n ( C F P RC B c ) ). We saw that we could rewrite the Boolean constraint as a disjunction of two by two exclusive disjuncts, B c = k D k , decomposing thus the solution space C F C RC B c η (resp., C F P RC B c η ) into the disjoint solution spaces for each disjunct, C F C RC D k η (resp., C F P RC D k η ). The elementary flux vectors (i.e., faces of dimension one) of F C η that satisfy the constraint RC D k are obviously elementary flux vectors of C F C RC D k η and the reciprocal is also true: if a flux vector of the semi-open polyhedral cone C F C RC D k η is not an elementary flux vector of F C η , i.e., is not a face of dimension one, then it belongs to the interior of a face of C F C RC D k η of dimension at least two and is thus conformally decomposable in this face, i.e., is not elementary in C F C RC D k η . It follows that the elementary flux vectors of C F C RC B c η are made up of all the elementary flux vectors of the C F C RC D k η ’s. We can sum up the results regarding EFVs as:
For B c = k D k , E F V s ( C F C RC B c ) = η maximal in s i g n ( C F C RC B c ) E F V s ( C F C RC B c η ) = η k E F V s ( C F C RC D k η ) = η k E F V s ( F C η ) S o l RC D k .
This means that EFVs can be computed flux tope by flux tope and constraint-disjunct by constraint-disjunct. Moreover, the result holds also for EFPs (by considering vertices of F P η ) and EFVs of C F P RC B c .
However, for EFMs, we have to take care that a phenomenon similar to that described in the example above still arises and that, even in a given flux tope, an EFM of C F C RC D k η is not necessarily an EFM of C F C RC B c η .
Example 2.
Consider the network comprising three irreversible reactions and one internal metabolite A: R 1 : A , R 2 : A , R 3 : A , and the Boolean constraint B c = ¬ R 1 R 2 , decomposed as B c = D 1 D 2 , with D 1 = ¬ R 1 and D 2 = R 1 R 2 (see Figure 5). Take η = + . Then, e 1 = ( 0 , 1 , 1 ) T , made up of R 2 , R 3 , which is the only EFM of C F C RC D 1 η , is also the only EFM of C F C RC B c η . However any e 2 = ( λ , 1 , 1 + λ ) T with λ > 0 , made up of R 1 , R 2 , R 3 , is an EFM of C F C RC D 2 η and is not an EFM of C F C RC B c η . Note that the way the constraint is decomposed matters. With the decomposition B c = D 1 D 2 with D 1 = R 2 and D 2 = ¬ R 1 ¬ R 2 , the result for D 1 is identical to that for D 1 , but C F C RC D 2 η = { 0 } .
This means that, if it is natural to study each C F C RC D k η (resp., C F P RC D k η ) separately in order to characterize the solution space and if EFVs are obtained in this way by collecting all the local EFVs, it is not the case for EFMs and, after collecting all local EFMs, we must only keep those with minimal support:
E F M s ( C F C RC B c η ) = Min supp { E F M s ( C F C RC D k η ) } k .
Proposition 8.
Given an arbitrary regulatory constraint RC B c with B c = k D k , where the disjuncts D k are taken two by two inconsistent, and the associated constrained flux cone subset C F C RC B c (resp., constrained flux polyhedron subset C F P RC B c ), its EFVs (resp., EFPs and EFVs) are obtained by collecting these for each flux tope (i.e., in each r-orthant O η ) and for each disjunct, i.e., for each C F C RC D k η (resp., C F P RC D k η ), and are nothing else than the EFVs (resp., EFPs and EFVs) of F C (resp., F P ) satisfying the constraint RC B c . This is not the case for EFMs. First, an EFM of C F C RC B c η is not necessarily an EFM of C F C RC B c , but actually the biologically relevant minimality concept being sign-minimality and not support-minimality, we will consider EFMs for each flux tope C F C RC B c η . Second, an EFM of C F C RC D k η is not necessarily an EFM of C F C RC B c η . EFMs of C F C RC B c η are actually obtained by collecting EFMs of C F C RC D k η for all k and keeping only those with minimal support.
This being said, we will now focus on an arbitrary disjunct D of the form i I v i j J ¬ v j with I , J { 1 , , r } , I J = . Thus, C F C RC D η (resp., C F P RC D η ) is the semi-open face F of F C η (resp., F P η ), obtained from the face F defined by { v j = 0 , j J } (i.e., F = F C η j J { v j = 0 } , idem with F P ) by removing its facets { v i = 0 } for all i I . Let’s note EFMs RC D = EFMs ( F C η ) S o l RC D those EFMs of F C η that satisfy the constraint RC D and EFMs RC D = EFMs ( C F C RC D η ) the EFMs of the part of the solution space in O η . Obviously, EFMs RC D EFMs RC D . If I = , F = F and EFMs RC D = EFMs RC D , so we will consider the case I . In this case, and contrary to what happens for thermodynamic and kinetic (as described in Proposition 5) constraints, there is generally no longer identity between EFMs RC D and EFMs RC D [42].
Example 3.
Consider the network of Example 1 (thus D = T 1 T 4 ) and let η = + and e 2 = ( 0 , 1 , 0 , 1 , 0 ) T be the EFM of F C η made up of T 1 and T 3 . Then, for any λ > 0 , e 2 + λ e 3 , positive conformal combination of the two EFMs e 2 and e 3 of F C η , has support { T 1 , T 2 , T 3 , T 4 } and belongs to E F M s RC D \ E F M s RC D (see Figure 6).
EFMs RC D correspond to the faces of dimension one (edges or extreme rays), if any, of the semi-open face F , i.e., the edges of the face F that have not been removed, which means that they are not included in any hyperplane of equation v i = 0 , for a certain i I , or equivalently that their supports contain I. Now, consider a face F of F, of dimension at least two, such that F ˚ F but no facet of F has its interior included in F (i.e., any facet of F is included in an hyperplane of equation v i = 0 , for a certain i I ), if any. Note that several such F may exist, but none can be included in another one, i.e., they are not comparable for inclusion in the lattice of the faces of F. Then, all vectors of F ˚ have the same minimal support, i.e., F ˚ EFMs RC D \ EFMs RC D . If { e k } k K (with | K | 2 ) are representatives of the extreme vectors of F , we have F = c o n e ( { e k } ) (which is actually the same as c o n e ( { e k } ) as we are in orthant O η ) and thus F ˚ = { k K β k e k β k > 0 } = c o n e + ( { e k } ) and the common minimal support of all vectors of F ˚ is s u p p ( F ˚ ) = k K s u p p ( e k ) . Conversely, if v EFMs RC D , let F be the minimal face of F containing v . If F has dimension one (extreme ray), then F \ { 0 } F and v EFMs RC D . If F has dimension at least two, then no facet of F has its interior included in F , because if it were the case for one facet, then any vector of its interior would have its support strictly included in the support of v , which would contradict the minimality of the latter. Finally, v F ˚ F and all vectors of F ˚ have the same support as v and belong to EFMs RC D \ EFMs RC D . We have thus characterized both EFMs RC D and EFMs RC D \ EFMs RC D . We stipulate now, for the latter one, the decomposition of its vectors into e k EFMs ( F ) \ EFMs RC D , in order to compute their supports, which is generally the only useful information (the precise decomposition into fluxes not often being very relevant). Therefore, we consider a face F of F, of dimension at least two, such that F ˚ F but no facet of F has its interior included in F and { e k } 1 k N a minimal set of vectors in EFMs ( F ) such that s u p p ( F ˚ ) = 1 k N s u p p ( e k ) . Note that, for any k, s u p p ( e k ) I (and, as we have seen, I s u p p ( e k ) ), because if for a certain k 0 we had s u p p ( e k 0 ) I = , then e k 0 would belong to all hyperplanes of equation v i = 0 for i I , thus to all facets of F , which is impossible for a non-null vector. Let us note S 1 = s u p p ( e 1 ) I and, for any k, 2 k N , S k = ( s u p p ( e k ) I ) \ 1 j k 1 S j . Then, for any k, S k , because if for a certain k 0 we had S k 0 = , then the vectors of c o n e + ( { e k } k k 0 ) would verify the constraint RC D and have their supports included in s u p p ( F ˚ ) , thus equal to it as it is minimal for vectors in C F C RC D η , which would contradict the minimality of { e k } . Finally, as by construction I = 1 k N S k and S k S j = for k j , we obtain that { S k } 1 k N constitutes a partition of I and { e k } 1 k N is a set of vectors in EFMs ( F ) EFMs ( F ) \ EFMs RC D such that s u p p ( e k ) S k by construction and s u p p ( e k ) S j for any j k , otherwise, by the same reasoning as above, e j could be suppressed from the set { e k } , contradicting its minimality. Finally, we obtain that the support of any vector in EFMs RC D \ EFMs RC D can be written as 1 k N s u p p ( e k ) , where { e k } 1 k N , N 2 , are vectors in EFMs ( F ) verifying s u p p ( e k ) S k and s u p p ( e k ) S j for all k and j k , where { S k } 1 k N is a partition of I (note that we have the same result for EFMs RC D by taking N = 1 ). Now, a given 1 k N s u p p ( e k ) is not necessarily minimal among the whole collection when we vary N , { e k } and { S k } . It is also possible that it strictly contains the support of a vector in EFMs RC D . Therefore, to obtain exactly the supports of vectors in EFMs RC D \ EFMs RC D , we must only keep the minimal elements for inclusion w.r.t. the whole collection extended by s u p p ( EFMs RC D ) .
Example 4.
Let us continue with the network of Examples 1 and 3, so D = T 1 T 4 and η = + . We have E F M s ( F ) = { e 1 , e 2 , e 3 } and E F M s RC D = { e 1 } (see Figure 4). The only partition of { T 1 , T 4 } with a size 2 is given by: S 1 = { T 1 } and S 2 = { T 4 } . The only vector in E F M s ( F ) whose support contains S 1 and not S 2 is e 2 and the only one whose support contains S 2 and not S 1 is e 3 . Thus, s u p p ( E F M s RC D \ E F M s RC D ) = { s u p p ( e 2 ) s u p p ( e 3 ) } = { { T 1 , T 3 } { T 2 , T 4 } } = { { T 1 , T 2 , T 3 , T 4 } } . Actually, we have from Example 3: E F M s RC D \ E F M s RC D = { e 2 + λ e 3 λ > 0 } (see Figure 6).
Consider now the following network comprising two internal metabolites and seven irreversible reactions, six of which are transfer reactions, R : A B , T 1 : A , T 2 : B , T 3 : A , T 4 : B , T 5 : A , T 6 : B (see Figure 7). Let D = T 1 T 2 and η = + . Let e 1 = ( 1 , 1 , 1 , 0 , 0 , 0 , 0 ) T made up of T 1 , R , T 2 , e 2 = ( 1 , 0 , 0 , 1 , 1 , 0 , 0 ) T made up of T 3 , R , T 4 , e 3 = ( 1 , 1 , 0 , 0 , 1 , 0 , 0 ) T made up of T 1 , R , T 4 , e 4 = ( 1 , 0 , 1 , 1 , 0 , 0 , 0 ) T made up of T 3 , R , T 2 , e 5 = ( 0 , 1 , 0 , 0 , 0 , 1 , 0 ) T made up of T 1 , T 5 , e 6 = ( 0 , 0 , 0 , 1 , 0 , 1 , 0 ) T made up of T 3 , T 5 , e 7 = ( 0 , 0 , 1 , 0 , 0 , 0 , 1 ) T made up of T 6 , T 2 and e 8 = ( 0 , 0 , 0 , 0 , 1 , 0 , 1 ) T made up of T 6 , T 4 . We have E F M s ( F ) = { e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 , e 8 } and E F M s RC D = { e 1 } . The only partition of { T 1 , T 2 } with a size 2 is given by S 1 = { T 1 } and S 2 = { T 2 } . The vectors in E F M s ( F ) whose support contains S 1 and not S 2 are e 3 and e 5 and those whose support contains S 2 and not S 1 are e 4 and e 7 . We have s u p p ( e 3 ) s u p p ( e 4 ) = { R , T 1 , T 2 , T 3 , T 4 } , s u p p ( e 3 ) s u p p ( e 7 ) = { R , T 1 , T 2 , T 4 , T 6 } , s u p p ( e 5 ) s u p p ( e 4 ) = { R , T 1 , T 2 , T 3 , T 5 } and s u p p ( e 5 ) s u p p ( e 7 ) = { T 1 , T 2 , T 5 , T 6 } . Each one of the four supports obtained is minimal in this collection, but the first three contain s u p p ( e 1 ) = { R , T 1 , T 2 } . Thus, s u p p ( E F M s RC D \ E F M s RC D ) = { { T 1 , T 2 , T 5 , T 6 } } . Actually, E F M s RC D \ E F M s RC D = { e 5 + λ e 7 λ > 0 } .
Finally, consider another network comprising three internal metabolites and seven irreversible reactions, five of which are transfer reactions, R 1 : B A , R 2 : C A , T 1 : B + C , T 2 : A , T 3 : A , T 4 : B , T 5 : C . Let D = R 1 R 2 T 2 and η = + (see Figure 8). Let e 1 = ( 1 , 1 , 1 , 2 , 0 , 0 , 0 ) T made up of T 1 , R 1 , R 2 , T 2 , e 2 = ( 1 , 1 , 1 , 0 , 2 , 0 , 0 ) T made up of T 1 , R 1 , R 2 , T 3 , e 3 = ( 1 , 0 , 0 , 1 , 0 , 1 , 0 ) T made up of T 4 , R 1 , T 2 , e 4 = ( 1 , 0 , 0 , 0 , 1 , 1 , 0 ) T made up of T 4 , R 1 , T 3 , e 5 = ( 0 , 1 , 0 , 1 , 0 , 0 , 1 ) T made up of T 5 , R 2 , T 2 and e 6 = ( 0 , 1 , 0 , 0 , 1 , 0 , 1 ) T made up of T 5 , R 2 , T 3 . We have E F M s ( F ) = { e 1 , e 2 , e 3 , e 4 , e 5 , e 6 } and E F M s RC D = { e 1 } . There are four partitions of { R 1 , R 2 , T 2 } with a size 2 . The partition { { R 1 } , { R 2 } , { T 2 } } gives nothing because there is no vector in E F M s ( F ) whose support contains { T 2 } but neither { R 1 } nor { R 2 } . For the partition S 1 = { R 1 , T 2 } and S 2 = { R 2 } , the vector in E F M s ( F ) whose support contains S 1 and not S 2 is e 3 and those whose support contains S 2 and not S 1 are e 2 , e 5 and e 6 , providing s u p p ( e 3 ) s u p p ( e 2 ) = { R 1 , R 2 , T 1 , T 2 , T 3 , T 4 } , s u p p ( e 3 ) s u p p ( e 5 ) = { R 1 , R 2 , T 2 , T 4 , T 5 } and s u p p ( e 3 ) s u p p ( e 6 ) = { R 1 , R 2 , T 2 , T 3 , T 4 , T 5 } . For the partition S 1 = { R 2 , T 2 } and S 2 = { R 1 } , the vector in E F M s ( F ) whose support contains S 1 and not S 2 is e 5 and those whose support contains S 2 and not S 1 are e 2 , e 3 and e 4 , providing s u p p ( e 5 ) s u p p ( e 2 ) = { R 1 , R 2 , T 1 , T 2 , T 3 , T 5 } , s u p p ( e 5 ) s u p p ( e 3 ) = { R 1 , R 2 , T 2 , T 4 , T 5 } and s u p p ( e 5 ) s u p p ( e 4 ) = { R 1 , R 2 , T 2 , T 3 , T 4 , T 5 } . For the partition S 1 = { R 1 , R 2 } and S 2 = { T 2 } , the vector in E F M s ( F ) whose support contains S 1 and not S 2 is e 2 and those whose support contains S 2 and not S 1 are e 3 and e 5 , providing s u p p ( e 2 ) s u p p ( e 3 ) = { R 1 , R 2 , T 1 , T 2 , T 3 , T 4 } and s u p p ( e 2 ) s u p p ( e 5 ) = { R 1 , R 2 , T 1 , T 2 , T 3 , T 5 } . The minimal elements of this collection of supports are { R 1 , R 2 , T 1 , T 2 , T 3 , T 4 } , { R 1 , R 2 , T 1 , T 2 , T 3 , T 5 } and { R 1 , R 2 , T 2 , T 4 , T 5 } . Minimizing also w.r.t. s u p p ( e 1 ) = { R 1 , R 2 , T 1 , T 2 } gives thus s u p p ( E F M s RC D \ E F M s RC D ) = { { R 1 , R 2 , T 2 , T 4 , T 5 } } . Actually, E F M s RC D \ E F M s RC D = { e 3 + λ e 5 λ > 0 } .
We have thus proved the following result.
Proposition 9.
Given a Boolean constraint D as a conjunction of literals i I v i j J ¬ v j and an arbitrary closed r-orthant O η , let C F C RC D η be the constrained flux cone subset for the regulatory constraint RC D in O η . Let F be the face of the flux cone F C η defined as F C η j J { v j = 0 } , then C F C RC D η is the semi-open flux cone F obtained from F by removing its facets { v i = 0 } for all i I . Let E F M s RC D = E F M s ( F C η ) S o l RC D be the EFMs of F C η that satisfy RC D and E F M s RC D = E F M s ( C F C RC D η ) be the EFMs of C F C RC D η . We get E F M s RC D E F M s RC D but, unlike the case of thermodynamic and kinetic (as in Proposition 5) constraints, there is generally no longer identity between E F M s RC D and E F M s RC D . E F M s RC D correspond to the extreme rays (faces of dimension one) of F , i.e., the edges of F whose support contains I. E F M s RC D \ E F M s RC D are the vectors of the F ˚ ’s for all F faces of F of dimension at least two, such that F ˚ F (i.e., F is not included in any hyperplane { v i = 0 } with i I ) but no facet of F has its interior included in F (i.e., any facet of F is included in a certain hyperplane { v i = 0 } with i I ). The supports of those vectors in E F M s RC D \ E F M s RC D are obtained as 1 k N s u p p ( e k ) , where { e k } 1 k N , N 2 , are vectors in E F M s ( F ) verifying s u p p ( e k ) S k and s u p p ( e k ) S j for all k and j k , where { S k } 1 k N is a partition of I (they are actually the minimal elements, for subset inclusion and including the supports of E F M s RC D when checking minimality, obtained like this, for any partition { S k } of I, of size at least two, and any choice of { e k } as above).
This characterization of EFMs for flux cones in the presence of regulatory constraints does not hold as simply for flux polyhedra, because the added inhomogeneous linear constraints ILC , given by G v h , may not respect the structure of C F C RC B c at all and may cut the interior of flux cones; we already mentioned that there is no direct relationship between EFMs and extreme or elementary fluxes for a flux polyhedron. Nevertheless, the ideas developed above for flux cones can be applied to flux polyhedra F P η , where certain of their faces have to be removed due to the constraint D because they are included in certain hyperplanes {vi = 0}, in order to determine EFMS RC D \ EFMS RC D . Moreover, in practice, ILC is generally only used to bound (below and/or above) fluxes, in which case each extreme ray + e of F C η is still partly present in F P η as for example an edge [α+]e, defined by the two vertices αe and α+e. The results above can then be transposed by using convex-conformal decomposition into vertices and conformal decomposition into elementary vectors to characterize EFMS RC D \ EFMS RC D .
We are now interested in looking for vectors in the solution space that are (in some sense) non-decomposable while not being support-minimal, and in characterizing them, if any. We could think of EFVs (or EFPs) but, from the study above regarding EFMs, we note that, for a flux cone constrained by the regulatory constraint D, the vectors in EFMs RC D \ EFMs RC D are necessarily conformally decomposable, as they can be described as the interiors of polyhedral cones in O η of dimension at least two (for example, e 2 + e 3 in Example 3 can be decomposed as ( 0.7 e 2 + 0.3 e 3 ) + ( 0.3 e 2 + 0.7 e 3 ) ). More straightforwardly, we can notice that EFMs RC D is equal to EFVs ( C F C RC D η ) , i.e., precisely the conformally non-decomposable vectors. As already pointed out, what matters more than the precise decomposition into fluxes is the decomposition of the supports (in the decomposition above of e 2 + e 3 , the supports of the components are unchanged, i.e., equal to s u p p ( e 2 ) s u p p ( e 3 ) = { T 1 , T 2 , T 3 , T 4 } ). It follows that a relevant concept for a solution vector is not to be conformally non-decomposable but, less strictly, to be (conformally) support-wise non-decomposable, in the sense that the vector cannot be conformally decomposed into two vectors of different (necessarily not greater) supports. Now, for all faces F of F of dimension at least two such that F ˚ F , not covered by Proposition 9, i.e., owning at least one facet F with interior included in F , it so happens that all vectors of F ˚ are actually support-wise decomposable, as each such vector can always be decomposed into a vector of F ˚ of same support and into a vector of F ˚ of smaller support. Therefore, in the decomposition, we do not authorize the same support for one of the component vectors. Thus the really relevant (less strict) concept is that the vector cannot be conformally decomposed into two vectors of smaller supports and we call a nonzero vector x of a convex polyhedral cone C as (conformally) support-wise non-strictly-decomposable, if
x = x 1 x 2 , with nonzero x 1 , x 2 C , implies s u p p ( x 1 ) = s u p p ( x ) or s u p p ( x 2 ) = s u p p ( x ) .
The (conformally) support-wise non-strictly-decomposable vectors in C (resp., flux vectors in a flux cone F C ) will be noted swNSDVs (resp., swNSDFVs). Note that we could have defined support-wise non-strictly decomposability more generally without imposing a conformal decomposition, but, as already pointed out, this is not relevant for biological fluxes and we will only use this concept in a certain tope. It is obvious that, in such a tope (resp., flux tope), EMs, ExVs = EVs ⊆ swNSDVs (resp., ExVs = EFVs = EFMs = swNSDFVs, so that all four definitions coincide). We will introduce a similar definition with a convex combination for a polyhedron P (the previous definition and the associated relationships above being valid for its recession cone C P ). A vector x of a polyhedron P is called support-wise convex(-conformally) non-strictly-decomposable, if
x = λ x 1 ( 1 λ ) x 2 , with x 1 , x 2 P , 0 < λ < 1 , implies s u p p ( x 1 ) = s u p p ( x ) or s u p p ( x 2 ) = s u p p ( x ) .
The support-wise convex(-conformally) non-strictly-decomposable vectors in P (resp., flux vectors in a flux polyhedron F P ) will be called support-wise non-strictly-decomposable points (resp., support-wise non-strictly-decomposable flux points) and noted swNSDPs (resp., swNSDFPs). It is obvious that, in any given tope (resp., flux tope), EMs, ExPs = EPs ⊆ swNSDPs (resp., EFMs, ExPs = EFPs ⊆ swNSDFPs).
With the notations of Proposition 9, let swNSDFVs RC D = swNSDFVs ( F C η ) S o l RC D and swNSDFVs RC D = swNSDFVs ( C F C RC D η ) . We have thus EFMs RC D = swNSDFVs RC D EFMs RC D swNSDFVs RC D but, unlike the case of thermodynamic and kinetic (as in Proposition 5) constraints for flux cones, not only there is no longer identity between swNSDFVs RC D and swNSDFVs RC D (consequence of the non-identity between EFMs RC D and EFMs RC D ), but we will now see that there is generally no longer identity between EFMs RC D and swNSDFVs RC D .
Example 5.
Continuing with the network of Example 3, e 1 + e 3 s w N S D F V s RC D \ E F M s RC D . More precisely, we have (by considering only a representative for each ray) E F M s RC D = s w N S D F V s RC D = { e 1 } E F M s RC D = E F M s RC D { e 2 + λ e 3 λ > 0 } s w N S D F V s RC D = E F M s RC D { e 1 + λ e 2 λ > 0 } { e 1 + λ e 3 λ > 0 } (see Figure 9).
Therefore, we extend the results of Proposition 9 by refining the structure in F of swNSDFVs RC D . Let { e m } m M be representatives of the extreme vectors of F, thus F = c o n e ( { e m } m M ) . We get the following structure for F regarding EFMs and swNSDFVs, given in the form of an algorithm.
  • Let R = { m M e m F } = { m M i I e i m 0 } . Then EFMs RC D = { e m } m R . Note that R can vary from to M, thus EFMs RC D from to EFMs ( F ) . If R = M , then EFMs RC D = EFMs RC D = swNSDFVs RC D = EFMs ( F ) and the analysis is done (if η has been chosen maximal in s i g n ( C F C RC D ) , this case corresponds to I = ). We consider the case R M here below.
  • Let us consider successively all faces F of F of dimension at least two, such that F ˚ F , i.e., F is not included in any hyperplane { v i = 0 } with i I (the lattice of faces of F can be explored for example in a way such that a sub-face is visited before a super-face; once such an F has been found, all faces of F containing it are also suitable). Let { e k } k K , with K M , be representatives of the extreme vectors of F . Thus, F ˚ = c o n e + ( { e k } k K ) and I k K s u p p ( e k ) . For each of these F , three exclusive cases can now be distinguished.
  • If no facet of F has its interior included in F , i.e., any facet of F is included in a certain hyperplane { v i = 0 } with i I (a necessary but insufficient condition is K M \ R ), then F ˚ EFMs RC D \ EFMs RC D and k K s u p p ( e k ) is a minimal support for vectors in F .
  • If exactly one facet of F has its interior included in F , i.e., it is the only facet not included in any hyperplane { v i = 0 } with i I , then F ˚ swNSDFVs RC D \ EFMs RC D and k K s u p p ( e k ) is a non-strictly-decomposable non-minimal support for vectors in F (non-strictly-decomposable support means that any vector which is a conical sum of the e k ’s having this support is support-wise non-strictly-decomposable, independently of the choice of the non-negative coefficients fixing the contribution of each e k in the distribution of the fluxes). This result follows immediately from the facts that one facet is not enough to decompose a certain vector in F ˚ strictly in terms of supports and that the support of the vectors in the interior of the facet in question is strictly included in the support of the vectors in F ˚ .
  • If at least two facets of F have their interior included in F , i.e., these facets are not included in any hyperplane { v i = 0 } with i I , then let { e l } l L , with L K , be representatives of the extreme vectors of all these facets (note that we have then necessarily K R L , i.e., K \ L K \ R ). Thus, the strict conical sum of the interiors of these facets, which is equal to c o n e + ( { e l } l L ) , is not empty in F ˚ (as there are at least two such facets) and is made up of the support-wise strictly-decomposable vectors of F ˚ (by construction): c o n e + ( { e l } l L ) = F ˚ \ swNSDFVs RC D . Two subcases must therefore be distinguished.
  • If L = K , i.e., the strict conical sum of the interiors of these facets is equal to F ˚ , then F ˚ F \ swNSDFVs RC D and k K s u p p ( e k ) is a strictly-decomposable support for vectors in F (which means that any vector which is a conical sum of the e k ’s having this support is support-wise strictly-decomposable, independently of the choice of the non-negative coefficients fixing the contribution of each e k in the distribution of the fluxes).
  • If L K , then F ˚ is split into two nonempty subsets: c o n e + ( { e l } l L ) F \ swNSDFVs RC D and F ˚ \ c o n e + ( { e l } l L ) swNSDFV RC D \ EFMs RC D (note that F ˚ \ c o n e + ( { e l } l L ) is made up of the vectors of c o n e + ( { e k } k K ) the conical decomposition of which on the e k ’s requires at least one e k with k K \ L ). This means that part of the vectors of F ˚ are support-wise non-strictly-decomposable and part are support-wise strictly-decomposable, while having the same non-minimal support k K s u p p ( e k ) (still equal to l L s u p p ( e l ) ). This proves that support-wise strict decomposability generally depends not only on the support, but also on the positive values of the fluxes and that, unlike the particular cases above, we cannot speak of a strictly-decomposable or non-strictly-decomposable support. Note that, in this subcase, the support-wise non-strictly-decomposable vectors of F ˚ constitute the complementary, in the open cone c o n e + ( { e k } k K ) , of the open sub-cone c o n e + ( { e l } l L ) which is thus a finite disjoint union of semi-open cones, each of which is conically generated (with positive or non-negative coefficients according to faces that are present or not) by extreme vectors e k ’s with either k K \ L K \ R or k L \ R , thus in any case with k K \ R , i.e., by extreme vectors e k EFMs ( F ) \ EFMs RC D .
We have finally proved the following result (keeping the notations of Proposition 9).
Proposition 10.
Let s w N S D F V s RC D = s w N S D F V s ( C F C RC D η ) be the support-wise non-strictly-decomposable vectors of C F C RC D η . We obtain E F M s RC D s w N S D F V s RC D but, unlike the case of thermodynamic and kinetic (as in Proposition 5) constraints, there is in general no longer identity between E F M s RC D and s w N S D F V s RC D . Consider all faces F of F of dimension at least two, such that F ˚ F (i.e., F is not included in any hyperplane { v i = 0 } with i I ), and let D ( F ) = c o n e + ( { F ˚ F facet of F with F ˚ F } ) . Result of Proposition 9 can be stated as E F M s RC D \ E F M s RC D = { F D ( F ) = } F ˚ . Now, we have s w N S D F V s RC D \ E F M s RC D = { F D ( F ) F ˚ } F ˚ \ D ( F ) . Note that the F ’s considered here necessarily own at least one facet F that is not included in any hyperplane { v i = 0 } with i I . There are actually two cases: For those F ’s which own exactly one such facet F (thus D ( F ) = F ˚ F \ F ˚ ), we obtain F ˚ s w N S D F V s RC D \ E F M s RC D and the common support of vectors in F ˚ is thus a non-strictly-decomposable (independently of the choice of the distribution of the fluxes) non-minimal support for vectors in F . For those F ’s which own at least two such facets F , we obtain F ˚ \ D ( F ) s w N S D F V s RC D \ E F M s RC D and D ( F ) F \ s w N S D F V s RC D , thus part of the vectors of F ˚ (consisting of a finite disjoint union of semi-open cones) are support-wise non-strictly-decomposable and part (consisting of an open cone) are support-wise strictly-decomposable (depending on the choice of the distribution of the fluxes), while having the same non-minimal support.
Example 6.
Let us continue with the network of Examples 3 and 5 (see Figure 9). F C is a pointed cone of dimension 3 included in the positive orthant of R 5 with axes { R , T 3 , T 2 , T 1 , T 4 } in this order. We have F C = { ( x y z x + y x + z ) T x , y , z 0 } . F C has 3 extreme rays (EFMs e 1 , e 2 , e 3 seen above) and 3 facets (see Figure 10). When we add the Boolean constraint D = T 1 T 4 , then F = F C and the solution space is the semi-open cone C F C RC D = F = { ( x y z x + y x + z ) T x , y , z 0 , x + y > 0 , x + z > 0 } , i.e., F without its two edges corresponding to e 2 and e 3 . The only edge remaining in F provides E F M s RC D = { e 1 } with support { T 1 , R , T 4 } . There are three faces of F of dimension two with their interiors included in F : F 1 = c o n e ( { e 2 , e 3 } ) , F 2 = c o n e ( { e 1 , e 3 } ) and F 3 = c o n e ( { e 1 , e 2 } ) . F 1 is the only one to have no facet with its interior included in F , i.e., such that D ( F 1 ) = , thus E F M s RC D \ E F M s RC D = F 1 ˚ = c o n e + ( { e 2 , e 3 } ) . This means that the new EFMs that appear are all the positive combinations of e 2 and e 3 , all with the same minimal support { T 1 , T 3 , T 2 , T 4 } . Regarding F 2 and F 3 , they both have R + e 1 as their only facet with interior included in F , thus D ( F 2 ) = D ( F 3 ) = R + * e 1 and F 2 ˚ = c o n e + ( { e 1 , e 3 } ) and F 3 ˚ = c o n e + ( { e 1 , e 2 } ) are both included in s w N S D F V s RC D \ E F M s RC D . { T 1 , R , T 2 , T 4 } and { T 1 , R , T 3 , T 4 } are thus non-strictly-decomposable (independently of the respective values of the fluxes in T 1 and T 2 or in T 3 and T 4 , respectively) non-minimal supports. The last face of F with its interior included in F is F itself, of dimension three, with F ˚ = c o n e + ( { e 1 , e 2 , e 3 } ) . As all the facets F 1 , F 2 , F 3 of F have their interior included in F , we obtain D ( F ) = F ˚ . Thus, the pathways with nonzero fluxes in all five reactions, i.e., with support { T 1 , T 3 , R , T 2 , T 4 } , are exactly the support-wise strictly-decomposable ones.
Example 7.
Let us consider the simple network comprising one reaction R : A B , where A and B are the two internal metabolites, and four transfer reactions T 1 : A , T 2 : A , T 3 : B , T 4 : B , and assume the five reactions irreversible (see Figure 11). F C is a pointed cone of dimension 3 included in the positive orthant of R 5 with axes { T 1 , T 2 , T 3 , R , T 4 } in this order. We have: F C = { ( x y z x + y x + y z ) T x , y , z 0 , x + y z } . F C has 4 facets and 4 extreme rays. Representatives of these extreme rays ( E F M s ) are e 1 = ( 1 0 1 1 0 ) T , e 2 = ( 1 0 0 1 1 ) T , e 3 = ( 0 1 1 1 0 ) T and e 4 = ( 0 1 0 1 1 ) T , defined by their supports: { T 1 , R , T 3 } , { T 1 , R , T 4 } , { T 2 , R , T 3 } and { T 2 , R , T 4 } , respectively.
Let us consider the Boolean constraint D = T 1 T 3 . Then, F = F C and the solution space is the semi-open cone C F C RC D = F = { ( x y z x + y x + y z ) T x , z > 0 , y 0 , x + y z } , i.e., F without its two facets { x = 0 } and { z = 0 } .
The only EFM still present in F is e 1 , thus E F M s RC D = { e 1 } with support { T 1 , R , T 3 } .
There are two faces of F of dimension two with their interiors included in F : F 1 = c o n e ( { e 1 , e 2 } ) and F 2 = c o n e ( { e 1 , e 3 } ) , both with R + e 1 as their only facet with interior included in F . Thus, D ( F 1 ) = D ( F 2 ) = R + * e 1 , E F M s RC D \ E F M s RC D = , F 1 ˚ = c o n e + ( { e 1 , e 2 } ) s w N S D F V s RC D \ EFMs RC D and F 2 ˚ = c o n e + ( { e 1 , e 3 } ) s w N S D F V s RC D \ E F M s RC D . { T 1 , R , T 3 , T 4 } and { T 1 , T 2 , R , T 3 } are thus non-strictly-decomposable (independently of the respective values of the fluxes in T 3 and T 4 or in T 1 and T 2 , respectively) non-minimal supports.
The last face of F with its interior included in F is F itself, of dimension three, with F ˚ = c o n e + ( { e 1 , e 2 , e 3 , e 4 } ) defined in F by { y > 0 , x + y > z } . F has exactly two facets with their interiors included in F , namely, F 1 and F 2 and thus D ( F ) = c o n e + ( { e 1 , e 2 , e 3 } ) = F \ s w N S D F V s RC D . The support-wise strictly-decomposable vectors constitute the open sub-cone of F ˚ defined by { y < z } : F \ s w N S D F V s RC D = { ( x y z x + y x + y z ) T x , y , z > 0 , y < z < x + y } . They are the pathways of support { T 1 , T 2 , R , T 3 , T 4 } such that the flux in T 2 is smaller than the flux in T 3 ; actually any vector ( x y z x + y x + y z ) T with x , y , z > 0 , y < z < x + y can be decomposed as ( k ( z y ) e 1 + ( x + y z ) e 2 ) ( ( 1 k ) ( z y ) e 1 + y e 3 ) , with arbitrary k, 0 < k < 1 , i.e., into two support-wise non-strictly-decomposable vectors respectively in F 1 ˚ with support { T 1 , R , T 3 , T 4 } and in F 2 ˚ with support { T 1 , T 2 , R , T 3 } .
We thus have F ˚ \ D ( F ) = c o n e + ( { e 2 , e 3 , e 4 } ) c o n e + ( { e 2 , e 3 } ) s w N S D F V s RC D \ E F M s RC D . The support-wise non-strictly-decomposable vectors constitute the semi-open sub-cone of F ˚ defined by { z y } : { ( x y z x + y x + y z ) T x , y , z > 0 , z y } . They are the pathways of support { T 1 , T 2 , R , T 3 , T 4 } such that the flux in T 2 is not smaller than the flux in T 3 . Any vector ( x y z x + y x + y z ) T with x , y , z > 0 , z y is equal to x e 2 + z e 3 + ( y z ) e 4 and belongs to c o n e + ( { e 2 , e 3 , e 4 } ) if z < y and to c o n e + ( { e 2 , e 3 } ) if z = y .
However, in the case of a general constraint B c = k D k , take care that, if the support-wise strictly-decomposable vectors of each C F C RC D k η are support-wise strictly-decomposable in C F C RC B c η , it is not true for support-wise non-strictly-decomposable vectors and some of them for C F C RC D k η may be decomposable in C F C RC B c η .
This characterization of support-wise non-strictly-decomposable vectors in the presence of regulatory constraints does not extend directly from flux cones to flux polyhedra. The reason is that, due to inhomogeneous linear constraints ILC , certain faces of F P η are not defined by equalities of the form v k = 0 and thus play no role in the definition of the support of the vectors they contain. Consequently, the basis of the reasoning above, namely, that the support of vectors of the interior of any face is larger than the supports of vectors of the interior of any facet of this face, does not hold. Nevertheless, the ideas developed for flux cones for dealing with faces included in a certain hyperplane { v i = 0 } with i I , and thus removed from the solution space due to the Boolean constraint considered, can be applied. However, it is necessary at each step, when considering an arbitrary face, to distinguish its facets resulting from ILC , not involved in the definition of the support, its facets included in a certain hyperplane { v k = 0 } with k I that are still present and contribute to the definition of the support, and its facets included in a certain hyperplane { v i = 0 } with i I that are removed.

2.2.4. Case of Several Types of Constraints

In general, when analyzing a metabolic pathway, all known biological constraints will have to be taken into account together, typically kinetic constraints and regulatory constraints. For two such constraints C 1 (say, a kinetic constraint KC , equivalent to TC , or KC b in the absence of bounds on enzyme concentrations, equivalent to TC b ) and C 2 (say a regulatory constraint RC B c ), the solution space, in the case of a flux cone F C (the reasoning would be similar for a flux polyhedron F P ) is given by C F C C 1 C 2 = F C S o l C 1 C 2 = F C S o l C 1 S o l C 2 = C F C C 1 S o l C 2 . Now, from Propositions 2 and 5, C F C C 1 is a finite union of flux cones F C i , which are certain particular faces of each flux tope of F C and the constraint C 2 can be applied to each one as to an original flux cone: C F C C 1 C 2 = i F C i S o l C 2 = i C F C i C 2 . From Proposition 7, each C F C i C 2 is a disjoint union of particular open faces of F C i . In all, the solution space is a disjoint union of particular open faces of the flux topes of F C . From propositions above and Proposition 8, we can also conclude that the EFVs of C F C C 1 C 2 are exactly the EFVs of F C that satisfy both constraints C 1 and C 2 .
Proposition 11.
The space of flux vectors in F C (resp., F P ) satisfying both the kinetic constraint KC (or KC b in the absence of bounds on enzyme concentrations) and the regulatory constraint RC B c is a finite disjoint union of open polyhedral cones (resp., open polyhedra) which are certain open faces of the flux topes of F C (resp., F P ). The elementary flux vectors (resp., elementary flux points and vectors) are those of F C (resp., F P ) that satisfy both constraints.
Note that Proposition 10 applies also, by starting from each F C i instead of each flux tope of F C , to determine elementary flux modes and support-wise non-strictly-decomposable vectors.

3. General Case of Sign-Compatible Constraints

We are interested in determining, for general biological constraints C , what can be said about the structure of the solution spaces C F C C or C F P C . More precisely, in identifying certain general features regarding constraints C (some of which are present in particular in biological constraints Equations (42), (47) and (50) we have considered so far), allowing us to clarify the mathematical and geometrical structure of the solution space, to determine the EFs (conformal non-decomposable fluxes) or EFMs (support-minimal fluxes) for this space and whether they characterize it and lastly how to compute them efficiently, in particular by checking the integration of this computation into the DD method. We will focus on deducing pertinent geometrical characteristics of the solution space only from general properties of the compatibility of the constraints with vector signs (i.e., with vector supports in each closed r-orthant or flux tope) and we will see that part of the results obtained above for thermodynamic, kinetic and regulatory constraints actually depends only on these general global properties.

3.1. Sign-Invariant Constraints

A constraint C x ( v ) (resp., x C x ( v ) ), is said to be sign-invariant if it depends only on s i g n ( v ) , i.e., on the signs of the v i ’s but not on their values (in the formulas below, free variables are assumed universally quantified):
C x ( v ) , s i g n ( v ) = s i g n ( v ) C x ( v )
x C x ( v ) , s i g n ( v ) = s i g n ( v ) x C x ( v ) .
It follows from these definitions that, if the C x ( v ) ’s are sign-invariant for all x , then x C x ( v ) is sign-invariant. Note that the property for a constraint to be support-invariant (i.e., to depend only on s u p p ( v ) ) is stronger than the property to be sign-invariant. Actually, in any flux tope (or in any given closed orthant), having the same sign for two vectors is equivalent to having the same support and thus being sign-invariant for a constraint means being support-invariant in each flux tope independently.
Example 8.
It follows then from definitions Equations (41), (42) or (43) that the thermodynamic constraint TC M ( v ) or TC ̲ M ̲ ( v ) is sign-invariant. Thus this is also the case for TC ( v ) and TC ̲ ( v ) Equation (60) and for TC b ( v ) and TC ̲ b ( v ) Equation (61).
On the other hand, the kinetic constraint KC E , M ( v ) Equation (47) is not sign-invariant, and this is also the case of KC E ( v ) Equation (78), because κ i ( M ) is bounded below and above (e.g., for Michaelis–Menten kinetics, k i < κ i < k i + ) and thus, for a given nonzero enzyme concentration E i , there is no metabolite concentrations vector M allowing an arbitrary flux value v i in reaction i (as soon as the value of v i / E i is outside the κ i bounds).
The kinetic constraint KC M ( v ) Equation (80) is however sign-invariant. This follows from the linear dependency of v i on E i . Actually, if KC M ( v ) is satisfied (for a certain E ) and s i g n ( v ) = s i g n ( v ) , then, for all i s u p p ( v ) , v i and κ i ( M ) have the same sign; thus, it is also the case for v i and κ i ( M ) . Therefore, by taking E i = v i / κ i ( M ) = E i v i / v i when v i 0 , and E i = 0 when v i = 0 , KC M ( v ) is satisfied (for E ). This sign-invariant property thus holds also for KC ( v ) Equation (76) and for KC b ( v ) Equation (77) if the only bounds are on metabolite concentrations. However, it does not hold for KC M b ( v ) (81) and for KC b ( v ) (77) in the presence of bounds on enzyme concentrations.
Finally, by definition, the regulatory constraint RC B c ( v ) Equation (50) depends only on s u p p ( v ) , so is support-invariant and thus sign-invariant.
Lemma 4.
All thermodynamic constraints are sign-invariant. Only kinetic constraints KC M and KC are sign-invariant, as well as KC b with bounds only on metabolite concentrations (but not on enzyme concentrations, thus in particular without enzymatic resource constraint). The regulatory constraints are support-invariant, thus sign-invariant.
The structure of the solution space of a sign-invariant constraint follows directly from its definition. Actually, if v satisfies a given sign-invariant constraint C = C x ( v ) (resp., C = x C x ( v ) ), then O ˚ s i g n ( v ) = { x R r s i g n ( x ) = s i g n ( v ) } is included in S o l C and thus S o l C = { v C ( v ) } O ˚ s i g n ( v ) is a disjoint union of open orthants. This result can be compared to Lemma 2. More precisely, if we consider a support-invariant constraint C , the same reasoning applies and gives S o l C = { v C ( v ) } S s u p p ( v ) , where S ρ = { x R r s u p p ( x ) = ρ } denotes the set of vectors of support ρ , where we code a support as a binary vector ρ (1 coding support membership and 0 non-membership). A support-invariant constraint C is thus equivalent (in extension) to a family { S ρ i } i of such support sets and, as there are 2 r possible binary vectors, there are thus 2 2 r different support-invariant constraints. Now, a given S ρ is equivalent to the Boolean constraint D ρ = { i ρ i = 1 } v i { j ρ j = 0 } ¬ v j and thus an arbitrary support-invariant constraint C is equivalent to the Boolean constraint B c = i D ρ i , thus to the regulatory constraint RC B c . Thus support-invariant constraints identify with regulatory constraints and all properties we have demonstrated for the latter apply to the first. For example, by associating with any sign vector η its support ρ (i.e., the support of any vector having the given sign), given by ρ i = 1 if η i = + or − and ρ i = 0 if η i = 0 , we deduce that, for any binary vector ρ , there are 2 | s u p p ( ρ ) | sign vectors η with support ρ , given by η i = + or − if ρ i = 1 and 0 else, and we obtain S ρ = { η i s u p p ( η i ) = ρ } O ˚ η i and S o l C = { η i s u p p ( η i ) = s u p p ( v ) C ( v ) } O ˚ η i , i.e., that S o l C is a disjoint union of open orthants, which is precisely the statement of Lemma 2. Note that a support-invariant constraint C is entirely defined by its restriction to an arbitrary closed r-orthant O η , i.e., for η an arbitrary full support sign vector, as S o l C can be reconstituted from S o l C O η . Moreover, as to be sign-invariant and to be support-invariant coincide in any closed r-orthant, a sign-invariant constraint is nothing but a constraint whose restriction to any closed r-orthant is support-invariant, i.e., which coincides in any closed r-orthant with a well-defined regulatory constraint (but such regulatory constraints differ in general from one r-orthant to another, while obviously coinciding on their intersection). A sign-invariant constraint C is equivalent to a family { O ˚ η i } i of open orthants and, as there are 3 r possible sign vectors, there are thus 2 3 r different sign-invariant constraints. Applying the merging method described in the proof of Lemma 2, these open orthants can be grouped together in order to obtain a family of semi-open orthants, without any possible further merging between any two of them.
Proposition 12.
Support-invariant constraints coincide with Boolean constraints, i.e., regulatory constraints. They are completely characterized by their restriction to any closed r-orthant and number 2 2 r . Sign-invariant constraints coincide with constraints whose restriction to any closed r-orthant is given by a regulatory constraint (with identity of such regulatory constraints on the intersections of any two of these orthants). They number 2 3 r . The set S o l C of vectors in R r satisfying the sign-invariant constraint C = C x ( v ) (resp., C = x C x ( v ) ) is a disjoint union of open orthants, which can be grouped together according to Lemma 2 to provide a disjoint union of semi-open orthants without any possible further merging between any two of them.
An important consequence is that, if we reason for each closed r-orthant separately, i.e., flux tope by flux tope in the solution space, then all results demonstrated in Section 2.2.3 for regulatory constraints apply to sign-invariant constraints. With the previous definitions of open or semi-open polyhedral cones and open or semi-open polyhedra and, more generally, of open or semi-open faces of polyhedral cones or polyhedra, we get the following result.
Theorem 1.
Let C = C x ( v ) (resp., C = x C x ( v ) ) be a sign-invariant constraint and C F C C (resp., C F P C ) be the associated constrained flux cone subset (resp., the associated constrained flux polyhedron subset). Then, C F C C (resp., C F P C ) is a finite disjoint union of open polyhedral cones (resp., open polyhedra), which are certain open faces of the flux topes F C η (resp., F P η ) for all η maximal sign vectors in s i g n ( F C ) (resp., s i g n ( F P ) ). They can be grouped together according to Lemma 2 to provide a disjoint union of semi-open faces of the F C η ’s (resp., F P η ’s) without any possible further merging between any two of them. Elementary fluxes are obtained by collecting those for each flux tope (which is not the case for elementary flux modes where, after collecting them, only those with minimal support have to be kept) and are nothing but the elementary fluxes of F C (resp., F P ) that satisfy the constraint C (which again is not the case for elementary flux modes):
E F V s ( C F C C ) = E F V s ( F C ) S o l C = E F M s ( F C ) S o l C E F M s ( C F C C )
E F P s ( C F P C ) = E F P s ( F P ) S o l C E F V s ( C F P C ) = E F V s ( C F P ) S o l C .
Proof of Theorem 1.
The disjoint decomposition of the solution space into open faces of the flux topes F C η (resp., F P η ) results from Proposition 12 or directly from the fact that all vectors of the strict conical hull c o n e + ( v 1 , , v n ) = { β 1 v 1 + + β n v n β 1 , , β n > 0 } and of the strict convex hull c o n v + ( v 1 , , v n ) = { α 1 v 1 + + α n v n α 1 , , α n > 0 , α 1 + + α n = 1 } of given vectors v 1 , , v n in a flux tope, thus in a closed orthant (so the sum is conformal), have the same support s u p p ( v 1 ) s u p p ( v n ) , thus the same sign, and from the fact that an open face of a polyhedron is the Minkowski sum of the strict convex hull of its vertices and the strict conical hull of its extreme vectors. Thus, an open face of a flux tope F C η (resp., F P η ) either does not intersect the solution space or is completely included in it. We have seen (55) that any elementary flux or any elementary flux mode of F C (resp., F P ) that satisfies the constraint is an elementary flux or an elementary flux mode of the solution space (the conditions of validity of (55) for elementary fluxes are satisfied from the structure of the solution space). These elementary fluxes or elementary flux modes actually identify with the extreme vectors of the F C η ’s (resp., extreme points of the F P η ’s and extreme vectors of the C F P η ’s) that satisfy the constraint, i.e., whose sign satisfies the constraint. Reciprocally, any vector of the solution space that is not in this case necessarily belongs to an open face of a certain F C η (or C F P η ) of dimension at least two (resp., of a certain F P η of dimension at least one) and is thus conformally (resp., convex-conformally) decomposable in this open face and consequently cannot be an elementary vector (resp., elementary point) of the solution space (obviously, the decomposition involves two vectors of same support, thus nothing can be deduced regarding elementary flux modes). This gives the result for elementary fluxes. □
Theorem 1 applies in particular to all thermodynamic constraints, regulatory constraints and those kinetic constraints described in Lemma 4. Especially, we directly obtain the results of Propositions 7 and 8 for regulatory constraints.
Note that our only knowledge of the EFVs for C F C C (resp., the EFPs and EFVs for C F P C ) is really a long way from characterizing the solution space C F C C (resp., C F P C ): actually they are just the ExVs of each tope C F C C η (resp., the ExPs and ExVs of each C F P C η ), i.e., the only one-dimension open polyhedral cones, i.e., edges (resp., zero-dimension polyhedra, i.e., vertices, and edges for the recession cones) among all the open cones (resp., open polyhedra) of any dimension that constitute C F C C (resp., C F P C ).
As in each closed r-orthant, i.e., in each flux tope, a sign-invariant constraint identifies with a regulatory constraint, Propositions 9 and 10 demonstrated for regulatory constraints apply thus to sign-invariant constraints. Nevertheless, they are stated for constraints that are conjunctions of literals and, even if a decomposition into such disjuncts is always possible from the identity above, it is not necessarily natural for an arbitrary sign-invariant constraint and a global characterization for each flux tope is preferable, in particular for what concerns support-wise non-strictly-decomposable vectors. The difference is that one has to deal for C F C C η with an arbitrary family of open faces of F C η instead of a single semi-open face of F C η , obtained specifically as a face of F C η without certain of its facets, for C F C RC D η .
Let us adopt notations similar to those used in Propositions 9 and 10, i.e., let EFMs C = EFMs ( F C η ) S o l C be the elementary flux modes of F C η that satisfy C (i.e., the elementary flux vectors of C F C C η ), EFMs C = EFMs ( C F C C η ) be the elementary flux modes of C F C C η and swNSDFVs C = swNSDFVs ( C F C C η ) be the support-wise non-strictly-decomposable vectors of C F C C η . We have thus EFMs C EFMs C swNSDFVs C . The proof of Proposition 9 adapts straightforwardly: EFMs C correspond to the open faces of dimension one (extreme rays) of F C η in C F C C η and EFMs C \ EFMs C are made up of all the open faces F ˚ in C F C C η where F is a face of dimension at least two of F C η but not any proper (i.e., with positive dimension, less than the dimension of F) face F of F is such that F ˚ is in C F C C η . The result is based on the properties that all vectors of an open face F ˚ in C F C C η have the same support, and the common support of vectors of F ˚ , where F is a face of F, is strictly included in that of vectors of F ˚ . In the same way, the proof of Proposition 10 is easily adapted for what concerns swNSDFVs C \ EFMs C . For this we consider successively all the open faces F ˚ in C F C C η where F is a face of dimension at least two of F C η which owns at least one proper face F such that F ˚ is in C F C C η . For each such given F ˚ , let { F j } j J be the family of all such F . Let { e k } k K be representatives of the extreme vectors of F and { e l } l L , with L K , be representatives of the extreme vectors of all these F j . Thus, c o n e + ( { e k } k K ) = F ˚ and c o n e + ( { e l } l L ) = c o n e + ( { F ˚ j } j J ) and we obtain F ˚ c o n e + ( { e l } l L ) = F ˚ \ swNSDFVs C , as by construction those vectors of F ˚ which belong to c o n e + ( { F ˚ j } j J ) are precisely those in F ˚ which are support-wise strictly-decomposable in C F C C η , the decomposition being achieved along vectors in the F ˚ j ’s, whose support is strictly included in the common support of vectors in F ˚ . Therefore, three different exclusive cases can be distinguished for F ˚ :
  • c o n e + ( { e l } l L ) F ˚ = , i.e., the F j ’s are all included in a same facet of F, that is in a same hyperplane { v i = 0 } for a certain coordinate i of the affine span of F, which is still equivalent to l L s u p p ( e l ) k K s u p p ( e k ) (this is for example the case if | J | = 1 ). Then, F ˚ swNSDFVs C \ EFMs C is entirely made up of support-wise non-strictly-decomposable vectors and k K s u p p ( e k ) is a non-strictly-decomposable non-minimal support for vectors in C F C C η .
  • c o n e + ( { e l } l L ) = F ˚ , i.e., L = K . Then, F ˚ C F C C η \ swNSDFVs C is entirely made up of support-wise strictly-decomposable vectors and k K s u p p ( e k ) is a strictly-decomposable support for vectors in C F C C η .
  • c o n e + ( { e l } l L ) F ˚ , i.e., L K and not all F j ’s are included in a same facet of F, which means that, for all i coordinates of the affine span of F, there exists j J such that F j is not included in the hyperplane { v i = 0 } , or equivalently, there exists l L such that e i l 0 , which is still equivalent to l L s u p p ( e l ) = k K s u p p ( e k ) . Then, F ˚ is split into two nonempty subsets: c o n e + ( { e l } l L ) = F ˚ \ swNSDFVs C , made up of support-wise strictly-decomposable vectors, and F ˚ \ c o n e + ( { e l } l L ) swNSDFVs C \ EFMs C , made up of support-wise non-strictly-decomposable vectors (note that F ˚ \ c o n e + ( { e l } l L ) is made up of those vectors of F ˚ the conical decomposition of which on the e k ’s requires at least one e k with k K \ L ). As all vectors of F ˚ have the same non-minimal support k K s u p p ( e k ) , this proves that support-wise strict decomposability does not generally only depend on the support but also on the positive values of the fluxes and that, unlike the two particular cases above, we cannot speak of a strictly-decomposable or non-strictly-decomposable support. Note that the support-wise non-strictly-decomposable vectors of F ˚ are obtained as the complementary, in this open cone, of the open sub-cone c o n e + ( { e l } l L ) , which is thus a disjoint union of semi-open cones (its connected components), each of which being conically generated both by certain extreme vectors e k with k K \ L (that necessarily do not belong to EFMs C ), with non-negative coefficients, and by certain extreme vectors e k with k L (that may belong or not to EFMs C ), with positive coefficients (in order to keep the concerned facets common with c o n e ( { e l } l L ) ).
We have finally proved the following result (while keeping the notations of Theorem 1).
Theorem 2.
For C a sign-invariant constraint, F C η a flux tope and C F C C η the associated subset of the solution space (disjoint union of open faces of F C η from Theorem 1), let E F M s C = E F M s ( F C η ) S o l C be the elementary flux modes of F C η that satisfy C (equal to E F V s ( C F C C η ) ), p p l E F M s C = E F M s ( C F C C η ) be the elementary flux modes of C F C C η , and p p l s w N S D F V s C = p p l s w N S D F V s ( C F C C η ) be the support-wise non-strictly-decomposable vectors of C F C C η , with E F M s C E F M s C s w N S D F V s C . p p l E F M s C correspond to the open faces of dimension one (extreme rays) of F C η belonging to C F C C η . Consider now successively all the open faces F ˚ in C F C C η where F is a face of dimension at least two of F C η , and for each such given F, let { e k } k K be representatives of the extreme vectors of F and { e l } l L , with L K , be representatives of the extreme vectors of all proper (i.e., with positive dimension, less than the dimension of F) faces F of F such that F ˚ belongs to C F C C η . Thus, c o n e + ( { e k } k K ) = F ˚ and let D ( F ) = c o n e + ( { e l } l L ) . Consequently, if L = , then F ˚ p p l E F M s C \ E F M s C ; if L and l L s u p p ( e l ) k K s u p p ( e k ) , then F ˚ s w N S D F V s C \ E F M s C and k K s u p p ( e k ) is a non-strictly-decomposable non-minimal support for vectors in C F C C η ; if L = K , then F ˚ C F C C η \ s w N S D F V s C and k K s u p p ( e k ) is a strictly-decomposable support for vectors in C F C C η ; if L K and l L s u p p ( e l ) = k K s u p p ( e k ) , then D ( F ) = F ˚ \ s w N S D F V s C , a non-empty open cone, constitutes the support-wise strictly-decomposable vectors in F ˚ and F ˚ \ D ( F ) s w N S D F V s C \ p p l E F M s C , a non-empty finite disjoint union of semi-open cones, constitutes the support-wise non-strictly-decomposable vectors in F ˚ , all having the same non-minimal support (decomposability depending in this case not only on the support but also on the distribution of the fluxes). We finally obtain p p l E F M s C \ E F M s C and p p l s w N S D F V s C \ E F M s C by collecting these vectors for each F ˚ .
Note that even the knowledge of swNSDFVs C is not enough to completely reconstruct the tope C F C C η of the solution space.
In order to deal with the usual enzymatic resource constraint (more generally capacity constraints) in the kinetic constraints, we generalize the sign-invariance criterion. A constraint C x ( v ) (resp., x C x ( v ) ), is said to be contracting-sign-invariant if, when satisfied by one vector, it is satisfied by any vector having the same sign which is not greater (on each component), i.e., that belongs to the open rectangle parallelepiped defined by the null vector and the given vector:
C x ( v ) , s i g n ( v ) = s i g n ( v ) , i | v i | | v i | C x ( v )
x C x ( v ) , s i g n ( v ) = s i g n ( v ) , i | v i | | v i | x C x ( v ) .
Obviously a sign-invariant constraint is contracting-sign-invariant and, if the C x ( v ) ’s are contracting-sign-invariant for all x , then x C x ( v ) is contracting-sign-invariant.
Example 9.
The kinetic constraints KC M b ( v ) Equation (81) and KC b ( v ) Equation (77) are contracting-sign-invariant in the absence of positive lower bounds on enzyme concentrations (i.e., when E = 0 ). Actually, if KC M b ( v ) is satisfied (for a certain E verifying c T E W and E E + ) and s i g n ( v ) = s i g n ( v ) with | v i | | v i | for all i, we saw that, by taking E i = v i / κ i ( M ) = E i v i / v i when v i 0 , and E i = 0 when v i = 0 , then KC M ( v ) is satisfied (for E ). Now, as | v i | | v i | , we have E E and thus c T E c T E W and E E + , so KC M b ( v ) is satisfied.
Lemma 5.
The kinetic constraints KC M b and KC b are contracting-sign-invariant in the absence of positive lower bounds on enzyme concentrations (i.e., when E = 0 ).
We will call 0 -star domain a subset S D of R r which has the property that the whole open segment joining 0 to any element of S D is included in S D : v S D , 0 < λ < 1 λ v S D . Any cone is a 0 -star domain.
Theorem 3.
If C = C x ( v ) (or C = x C x ( v ) ) is contracting-sign-invariant, then C F C C is a 0 -star domain and there exists a neighborhood N = δ , + δ r of 0 for a certain δ > 0 such that C F C C N is a finite disjoint union of N-truncated open polyhedral cones, which are the intersection with N of open faces of the flux topes F C η . In addition, results (88) regarding E F V s and p p l E F M s hold in N, i.e., for sufficiently small fluxes. The same holds for C F P C with the flux topes F P η if 0 is an interior point of F P , and in this case results (89) regarding E F V s and E F P s hold in N.
Proof of Theorem 3.
The property of being a 0 -star domain is a direct consequence of the definition: if a contracting-sign-invariant constraint is satisfied by a vector v , it is satisfied by λ v for all 0 < λ < 1 . From the proof of Theorem 1, it follows that if a nonzero vector v F C η verifies a given contracting-sign-invariant constraint, then any vector v of the minimal open face of F C η containing v (thus having the same sign as v ) and belonging to δ v , + δ v r with δ v = m i n i s u p p ( v ) | v i | (contracting condition) also verifies the constraint. The result is obtained by taking for δ the minimum of the δ v ’s on all open faces of the F C η ’s (whose number is smaller than the number of possible signs, i.e., 3 r ). The proof is still valid for F P but the result is meaningful only if N is included in F P , i.e., if 0 is an interior point of F P , which means that the inhomogeneous linear constraints defining F P , given by G v h , verify h < 0 . □
Theorem 3 tells us that the result of Theorem 1 for sign-invariant constraints regarding the geometrical structure of the solution space applies identically for contracting-sign-invariant constraints locally in a neighborhood of 0 , i.e., when considering only pathways with sufficiently small amounts of fluxes. This applies in particular to those kinetic constraints described in Lemma 5.

3.2. Sign-Monotone Constraints

We now consider a property of compatibility of a constraint with signs that is stronger than sign-invariance. A constraint C x ( v ) (resp., x C x ( v ) ), is said to be sign-monotone if, when satisfied by a vector, it is satisfied by any other vector with a smaller or equal sign (for the partial order on signs):
C x ( v ) , s i g n ( v ) s i g n ( v ) C x ( v )
x C x ( v ) , s i g n ( v ) s i g n ( v ) x C x ( v ) .
It follows from these definitions that, if the C x ( v ) ’s are sign-monotone for all x , then x C x ( v ) is sign-monotone and that any sign-monotone constraint is sign-invariant. Note that, in any given closed orthant, thus in any flux tope, s i g n ( v ) s i g n ( v ) is equivalent to s u p p ( v ) s u p p ( v ) and thus being sign-monotone for a constraint means being support-monotone in each flux tope. Note also that if a sign-monotone constraint has a solution, then the null vector 0 is a solution.
Example 10.
It then follows from definitions (41), (42) or (43) that the thermodynamic constraint TC M ( v ) or TC ̲ M ̲ ( v ) is sign-monotone. This is thus also the case for TC ( v ) and TC ̲ ( v ) Equation (60) and for TC b ( v ) and TC ̲ b ( v ) Equation (61).
The argument given previously to establish that the kinetic constraint KC M ( v ) Equation (80) is sign-invariant proves in fact that it is sign-monotone. This sign-monotone property holds thus also for KC ( v ) Equation (76) and for KC b ( v ) Equation (77) in the absence of bounds on enzyme concentrations.
The regulatory constraint RC B c ( v ) Equation (50) is not generally sign-monotone: consider for example B c reduced to a positive literal. Actually, RC B c ( v ) is sign-monotone if and only if B c in DNF contains no positive literal, which is a very special case, requiring only certain fluxes to be zero but unable to express that a given flux is nonzero (in this particular case, it is support-monotone, which is stronger than sign-monotone).
Lemma 6.
All thermodynamic constraints are sign-monotone. Only those kinetic constraints KC M and KC are sign-monotone, as well as KC b in the absence of bounds on enzyme concentrations. The regulatory constraint RC B c is not sign-monotone, except if B c in DNF contains only negative literals (i.e., assigns only zero fluxes, but no nonzero fluxes), in which case it is support-monotone.
The structure of the solution space of a sign-monotone constraint follows directly from its definition.
Lemma 7.
The set S o l C of vectors in R r satisfying the sign-monotone constraint C = C x ( v ) (resp., C = x C x ( v ) ) is a union of closed orthants.
Actually, if v satisfies C , then O s i g n ( v ) = { x R r s i g n ( x ) s i g n ( v ) } is included in S o l C and thus S o l C = { v C ( v ) } O s i g n ( v ) . Obviously, we can keep only those v ’s for which s i g n ( v ) is maximal, in order to avoid any inclusion between the closed orthants. This result can be compared to Lemmas 1, 3 and first part of Proposition 5, which are obviously more precise but it shows that the mathematical structure of the solution space of these thermodynamic and kinetic constraints is mainly the only consequence of their sign-monotonicity.
Theorem 4.
Consider indifferently a sign-monotone constraint C x ( v ) or x C x ( v ) (e.g., if C x ( v ) is sign-monotone for any x ), and let us name it C and the associated constrained flux cone subset C F C C (resp., the associated constrained flux polyhedron subset C F P C ). Then, C F C C (resp., C F P C ) is a finite union of polyhedral cones (resp., polyhedra), which are faces of the flux topes F C η (resp., F P η ) for all η maximal sign vectors in s i g n ( F C ) (resp., s i g n ( F P ) ) and we obtain
E F V s ( C F C C ) = E F M s ( C F C C ) = E F V s ( F C ) S o l C = E F M s ( F C ) S o l C
E F P s ( C F P C ) = E F P s ( F P ) S o l C E F V s ( C F P C ) = E F V s ( C F P ) S o l C .
Proof of Theorem 4.
The (non-disjoint) decomposition of the solution space into faces of the flux topes F C η (resp., F P η ) follows from Lemma 7 or directly from the fact that all vectors of the conical hull c o n e ( v 1 , , v n ) = { β 1 v 1 + + β n v n β 1 , , β n 0 } and of the convex hull c o n v ( v 1 , , v n ) = { α 1 v 1 + + α n v n α 1 , , α n 0 , α 1 + + α n = 1 } of given vectors v 1 , , v n in a flux tope, thus in a closed orthant (so the sum is conformal), have their supports included in s u p p ( v 1 ) s u p p ( v n ) , which is the support of any vector v in c o n e + ( v 1 , , v n ) or c o n v + ( v 1 , , v n ) , thus their signs being less than or equal to the sign of v , and from the fact that a face of a polyhedron is the Minkowski sum of the convex hull of its vertices and the conical hull of its extreme vectors. Thus, if a vector v of a flux tope F C η (resp., F P η ) belongs to the solution space, the minimal face of F C η (resp., F P η ) containing v is completely included in it. As a sign-monotone constraint is sign-invariant, the results of Theorem 1 apply and provide formulas for EFVs and EFPs. Thus remains only the case of EFMs. Now, if an elementary flux mode v of the solution space were not support-minimal in F C , a nonzero vector v would exist in F C with s u p p ( v ) s u p p ( v ) . Moreover, we could choose λ R such that v = v λ v is nonzero, belongs to F C and verifies s u p p ( v ) s u p p ( v ) and s i g n ( v ) s i g n ( v ) . From the sign-monotonicity property of the constraint, v would also satisfy the constraint and would thus belong to the solution space, which would contradict the fact that v is an elementary flux mode in this space. Therefore, v is an elementary flux mode in F C and we get the result for EFMs. □
Of course, in the decomposition of C F C C ( x ) or C F C C (resp., C F P C ( x ) or C F P C ) as a union of certain faces of the F C η ’s (resp., F P η ’s), we can only keep those faces which are maximal. Theorem 4 applies in particular to all thermodynamic constraints, to those kinetic constraints and to those very few regulatory constraints described in Lemma 6. In particular, we directly obtain (except of course the reference to ts M ) Proposition 2 for thermodynamic constraints TC and TC b , as well as the structure of the solution space as a union of flux cones (resp., flux polyhedra) for kinetic constraints KC and KC b (in the absence of bounds on enzyme concentrations) and the characterization of elementary fluxes and elementary flux modes given by Proposition 5, which proves that these results depend only on the fact that these constraints are sign-monotone.
Note that the knowledge of EFVs, equal to EFMs, of C F C C (resp., and of EFPs of C F P C ) is not good enough to reconstruct the structure of the solution space C F C C (resp., C F P C ) as a union of polyhedral cones (resp., polyhedra). For this, it is necessary to know the (non-disjoint) decomposition { E i } of these EFVs (resp., and EFPs) as extreme vectors (resp., and extreme points) of the faces F i of the flux topes F C η (resp., F P η ) that constitute the solution space. The E i ’s are exactly the maximal subsets of EFVs (resp., and EFPs) whose conformal conical hull (resp., conformal convex hull) is entirely contained in the solution space: c o n e ( E i ) C F C C (resp., and c o n v ( E i ) C F P C ). By analogy with LTCS, we will call each such E i a largest C -consistent set of EFVs or EFMs (resp., and of EFPs), noted LCS C .
In order to deal with enzymatic capacity constraints, we generalize the sign-monotonicity criterion exactly in the same way we generalized the sign-invariance criterion. A constraint, C x ( v ) (resp., x C x ( v ) ), is said to be contracting-sign-monotone if, when satisfied by one vector, it is satisfied by any vector with a smaller or equal sign which is not greater (on each component) than the vector itself, i.e., by any vector that belongs to the rectangle parallelepiped defined by the null vector and the given vector:
C x ( v ) , s i g n ( v ) s i g n ( v ) , i | v i | | v i | C x ( v )
x C x ( v ) , s i g n ( v ) s i g n ( v ) , i | v i | | v i | x C x ( v ) .
Obviously a sign-monotone constraint is contracting-sign-monotone and, if the C x ( v ) ’s are contracting-sign-monotone for all x , then x C x ( v ) is contracting-sign-monotone. Furthermore, any contracting-sign-monotone constraint is contracting-sign-invariant.
Example 11.
The argument given in Example 9 to establish that the kinetic constraints KC M b ( v ) Equation (81) and KC b ( v ) Equation (77) are contracting-sign-invariant in the absence of positive lower bounds on enzyme concentrations (i.e., when E = 0 ), proves in fact that they are contracting-sign-monotone.
Lemma 8.
The kinetic constraints KC M b and KC b are contracting-sign-monotone in the absence of positive lower bounds on enzyme concentrations (i.e., when E = 0 ).
Theorem 5.
If C = C x ( v ) (or C = x C x ( v ) ) is contracting-sign-monotone, then C F C C is a 0 -star domain and there exists a neighborhood N = δ , + δ r of 0 for a certain δ > 0 such that C F C C N is a finite union of N-truncated polyhedral cones, which are the intersection with N of faces of the flux topes F C η . In addition, results (94) regarding E F V s and E F M s hold in N, i.e., for sufficiently small fluxes. The same holds for C F P C with flux topes F P η if 0 is an interior point of F P , and in this case results (95) regarding E F V s and E F P s hold in N.
Proof of Theorem 5.
The proof becomes identical to that of Theorem 3 by using the proof of Theorem 4. □
Theorem 5 tells us that the result of Theorem 4 for sign-monotone constraints regarding the geometrical structure of the solution space and the determination of elementary fluxes and elementary flux modes applies identically for contracting-sign-monotone constraints locally in a neighborhood of 0 , i.e., when considering only pathways with sufficiently small amounts of fluxes. It applies in particular to those kinetic constraints described in Lemma 8, providing the result for the structure of C F C KC M b and C F C KC b quoted in the last part of Proposition 5.
Finally, Proposition 11 extends straightforwardly to the case where we have to deal with both one sign-monotone constraint and one sign-invariant constraint. As a sign-monotone constraint is sign-invariant the conjunction of the two constraints is sign-invariant and Theorem 1 applies.
Theorem 6.
The space of flux vectors in F C (resp., F P ) satisfying both a sign-monotone constraint and a sign-invariant constraint is a finite disjoint union of open polyhedral cones (resp., open polyhedra) which are certain open faces of the flux topes of F C (resp., F P ). The elementary flux vectors (resp., elementary flux points and vectors) are those of F C (resp., F P ) that satisfy both constraints.
Note that, in the case of F C , Theorem 2 applies also to determine elementary flux modes and support-wise non-strictly-decomposable vectors.

3.3. Consequences on the Computation of Elementary Fluxes

From the results above, we will see that the computation of elementary fluxes in the presence of a sign-monotone constraint can be efficiently performed with the Double Description (DD) method.
First let us briefly remember the principle of the DD method [11]. This is an algorithm that takes as input the implicit description of a pointed convex polyhedral cone C as its representation matrix, i.e., a finite set of homogeneous linear inequalities defining C as the intersection of the corresponding vector half-spaces, and produces as output the explicit description of C as a (minimal) generating matrix, i.e., the set of extreme rays of C. More generally, it can deal in the same way with a pointed convex polyhedron P producing, from a finite set of linear inequalities defining P as the intersection of the corresponding affine half-spaces, the explicit description of P as two generating matrices providing respectively the vertices of P and the extreme rays of C P . As the first are obtained as the extreme rays of the pointed cone obtained from P by adding one dimension to the space and considering the conical hull of P from an origin of this extended vector space, it is enough to explain how the DD method works on pointed convex polyhedral cones. The DD method is an incremental algorithm that processes one by one each of the n homogeneous linear inequalities defining C. At each step i , 1 i n , it builds the intermediate extreme rays of the intermediate current cone C i defined by the i first linear inequalities from the knowledge of the extreme rays of C i 1 built at previous steps and of the ith linear inequality. At the end, for i = n , C n = C and the extreme rays are thus obtained. The ith linear inequality defines a half-space H i + with a vector hyperplane H i as a frontier. All the extreme rays of C i 1 that are on the right side of H i , i.e., belong to H i + , are extreme rays of C i and thus kept; all those that are on the wrong side do not belong to C i and will be discarded; the new extreme rays at step i appear from the intersection of C i 1 with H i and are obtained as the intersection with H i of the 2-faces of C i 1 defined by two adjacent extreme rays, one on each side of H i . We therefore just need to keep and update at each step this list of pairs of adjacent extreme rays and, when the ith linear inequality is processed, each such pair with elements on both sides of H i determines with H i an extreme ray for C i . The key fact is that each new extreme ray built at step i is the conical sum of two extreme rays existing at step i 1 and, as we will proceed by flux tope, i.e., in a given closed r-orthant, this sum is conformal and the support of this new extreme ray is the union of the supports of two previously existing extreme rays. Moreover, similarly, this new extreme ray, if involved in the next steps as a member of a relevant pair of adjacent rays, in order to build a new extreme ray, will have its support included in the support of the latter one. Finally, all we need to do is initialize C 1 . For a flux cone F C defined by Equation (6), what has been shown as the most efficient initialization is to start from a basis of the nullspace of S (the authors of [43] proposed a method to compute EFMs directly as linear combinations of the vectors of a basis of this nullspace but it turned out to be less efficient than DD), and this basis (of size r m by assuming S of full rank) can be chosen so that r m of the r linear inequalities are satisfied (we consider here the worst case corresponding to the highest number of such inequalities, i.e., r I = r , which happens in particular if each reversible reaction is split into two irreversible ones). There thus remain m linear inequalities to satisfy, i.e., there are n = m steps in the DD algorithm.
Let us now consider a sign-monotone constraint C . From Theorem 4, the EFVs (identical to EFMs) of the constrained flux cone subset C F C C are simply those of F C , i.e., of all its flux topes F C η successively, that satisfy constraint C (and the same is true for EFPs of C F P C ). Therefore, we just have to build the extreme rays of each F C η by using the DD algorithm as presented above and filter those that satisfy C . Obviously a filtering at the end, once all extreme rays built, will not at all improve the efficiency. However, by exploiting the key fact that any newly built extreme ray at a certain step, if used in other next steps to build other new extreme rays will have necessarily its support included in each support of the latter ones, and exploiting the sign-monotonicity of C , we can conclude that, each time a newly built extreme ray does not satisfy the constraint C , it can be immediately discarded because not only it has no use at the step of its discovery but also no use at the future steps as all extreme rays, in the building of which it could be involved, would have a larger support and thus would thus not satisfy C either. Thus, the computation of elementary fluxes or elementary flux modes in the presence of a sign-monotone constraint can be fully integrated into the incremental DD algorithm, the filtering of the extreme rays satisfying the constraint being achieved at each step when they are newly created. The extra-cost is having to check all intermediate extreme rays built for satisfiability by C (for most constraints, this can be practically instantaneous) and the gain is that all intermediate extreme rays for which this checking gives a negative answer are discarded.
Proposition 13.
Given a sign-monotone constraint C , the E F V s (or E F M s ) of the constrained flux cone subset C F C C are obtained by the DD algorithm as the extreme rays of each flux tope F C η by filtering out at each step all those newly built extreme rays that do not satisfy C .
This result applies in particular to thermodynamic constraints TC and TC b [44,45] and to kinetic constraints KC and KC b (in the absence of bounds on enzyme concentrations). Actually, for TC (or KC which is identical from Proposition 5), Proposition 3 demonstrated that the criterion for an extreme ray to be thermodynamically satisfiable is to belong to a given fixed open half-space delimited by a vector hyperplane. Checking this satisfiability thus boils down to computing the scalar product of the extreme ray with a fixed vector (normal to the hyperplane and oriented towards the half-space) and verifying it is positive. In this case, as it has been shown, it is much simpler to integrate this thermodynamic constraint as one supplementary ( n + 1 ) th homogeneous linear inequality and we can choose for example that it is processed first by the DD algorithm. However, as it has been already pointed out, providing the list of the extreme flux vectors of C F C C does not provide its structure in terms of a union of polyhedral cones. For this, we have to identify all the LCS C ’s, i.e., the maximal subsets of those extreme flux vectors whose conformal conical hull is entirely contained in C F C C . For each flux tope F C η , this amounts to computing all maximal upper bounds for sets of signs of those extreme rays built from this flux tope such that any vector of this tope with a sign equal to such an upper bound satisfies constraint C , as the vectors of the strict conical hull of a subset of extreme rays have as a common sign the upper bound for the signs of these extreme rays. To each such upper bound identified is associated the LCS C made up of all the extreme rays whose signs are smaller than, or equal to, this upper bound.
For sign-invariant constraints, and thus for regulatory constraints, no such incremental filtering is possible as a newly built intermediate extreme ray may have a sign that does not satisfy the constraint but may participate at some later stage in the construction of other new extreme rays (with necessarily larger supports) whose sign will satisfy the constraint. Thus, such an extreme ray cannot be discarded. In addition, we saw in any case there is no longer identity between those EFMs of F C that satisfy the constraint, the EFMs of C F C C and the support-wise non-strictly-decomposable vectors in C F C C and, as only the last two are biologically relevant, computing only the first ones would be of limited interest.

4. Conclusions

The analysis of metabolic networks in a steady state takes classically into account only stoichiometric and reactions irreversibility constraints. In this context, well-known pathways have been introduced, such as elementary flux modes EFMs, expressing the property of minimality, or elementary flux vectors EFVs, expressing the property of non-decomposability (which coincide in this standard framework), that are biologically relevant in their own and from which all pathways, whose structure is a convex polyhedral flux cone F C , can be reconstructed. However, first, the computation of these EFMs or EFVs is hampered by the combinatorial explosion of their number in genome-scale metabolic models, and, second, most of them are biologically irrelevant, because other important biological constraints are not taken into account. With the objective of both enumerating only biologically feasible EFMs or EFVs (without filtering them afterwards) and, as there are considerably fewer of them, improving the scalability of this computation, we took into consideration in this paper on one side thermodynamic and, more generally, kinetic constraints and on the other side regulatory constraints, and we tackled the problem of revisiting in this new extended framework the concept of EFMs and EFVs and, more largely, of characterizing the geometry of the solution space. Actually, we considered a more general conceptual framework for constraints, namely their property to be compatible with vector signs (i.e., with vector supports separately in each closed r-orthant), because most properties of the geometrical structure of the solution space happen to depend only on this very general sign-compatibility of constraints. This is how we demonstrated, for constraints which are sign-monotone (i.e., once satisfied by a certain vector are satisfied by any vector with smaller or equal sign), which is the case of thermodynamic constraints and of kinetic constraints in the absence of bounds on enzyme concentrations, that the solution space is a union of convex polyhedral cones, which are certain faces of the flux topes of F C , and that the EFMs, which still coincide with the EFVs, are simply those of F C that satisfy the constraint. In addition, we showed that their computation can be efficiently integrated into the classical double description algorithm, as each newly built intermediate extreme ray that does not satisfy the constraint can be filtered out at each incremental step. For the specific case of thermodynamic constraints or of kinetic constraints in the absence of bounds on enzyme concentrations, we demonstrated that their solution spaces are identical and, when there are also no bounds on metabolite concentrations, made up of those maximal faces of the flux topes of F C which are entirely contained in a fixed open vector half-space and thus that the EFMs are simply those of F C belonging to this half-space and thus computable by adding one single homogeneous linear inequality to the initial (stoichiometric and reactions irreversibility) ones. The situation is more complex, and challenging from a computational perspective, for constraints which are only sign-invariant (i.e., once satisfied by a certain vector are satisfied by any vector with the same sign), which is the case of regulatory constraints. We showed that in fact sign-invariant constraints are constraints that are support-invariant in each closed r-orthant separately and that support-invariant constraints coincide with regulatory (Boolean) constraints, therefore regulatory constraints are prototypical of sign-invariant constraints. For such sign-invariant constraints, we demonstrated that the solution space is a finite disjoint union of semi-open convex polyhedral cones, obtained from certain faces of the flux topes of F C by removing certain of their own faces of lesser dimension, and that the EFVs are simply those of F C (equal to the EFMs of F C ) that satisfy the constraint. However, these cannot be efficiently computed by the double description algorithm because they cannot be filtered out during the incremental process. In addition it is no longer true that EFVs identify with EFMs, as it exists in general minimal solutions that are decomposable, and that EFMs identify with support-wise non-strictly-decomposable vectors (a new property that we introduced and which is the proper one to express non-decomposability in this general framework, in terms of support non-decomposability), as there are in general such support-wise non-strictly-decomposable solutions that are not minimal: there are usually strict inclusions between these three sets and we provided again a complete characterization of the two latter ones. Finally, we extended all these results to the case where inhomogeneous linear constraints (expressing for example capacity constraints or bounds on fluxes) exist, dealing thus with a convex flux polyhedron F P instead of a flux cone F C . Basically, most of the results regarding the geometrical structure of the solution space in the presence of the above biological constraints remain the same with cones replaced by polyhedra.
Future work will be carried out along two paths. First, the present theoretical work will be extended to minimal cut sets (MCSs). Such MCSs are defined as minimal (for set inclusion) sets of reactions whose deletion will block the operation of given objectives or target reactions (as, e.g., those producing some toxic or undesirable product), i.e., removal of an MCS (that is, the knockout of its reactions) implies a zero flux for the target reactions in a steady state. MCSs are important for computing intervention strategies, e.g., for metabolic engineering. It is obvious from this definition that, for a metabolic network modeled in a steady state by a flux cone F C , the MCSs are the minimal hitting sets of the set of target EFMs (identical to EFVs of the given metabolic network (i.e., the set of EFMs that comprise at least one of the target reactions), where a hitting set of the target EFMs is a set (of reactions) that has a nonempty intersection with each one of these EFMs. Moreover, this generalizes to a metabolic network modeled by a flux polyhedron F P by using EFPs. This gives an indirect method for computing MCSs [46,47] from the preliminary computation of EFVs (or EFPs) that can be applied in the presence of biological constraints by using the results obtained in this paper. However, it is known that there is also a method for computing MCSs directly as the EFVs of a dual network [48,49], obtained basically by transposing the stoichiometric matrix of the original matrix (thus in some sense, exchanging reactions and metabolites). We will therefore study how to define this dual operation properly for metabolic networks in the presence of biological constraints in order to preserve this result (or most of it). Second, the algorithms described in this paper will be implemented and testing conducted on metabolic networks described in the literature for which biological (thermodynamic, kinetic or regulatory) data are known. As it has been shown, the computation of elementary flux vectors (or elementary flux points) boils down to filtering those of extreme rays (or vertices) in each flux tope which satisfy the biological constraints. Moreover, for sign-monotone constraints as thermodynamic constraints or kinetic constraints in the absence of bounds on enzyme concentrations, this filtering can be easily integrated at each step of the incremental double description method, so it is natural to rely first on this method which benefits from very efficient implementations. In the absence of bounds on metabolite concentrations, the advantage is obvious as it is just necessary to add one step (i.e., one homogeneous linear constraint to deal with) in the DD algorithm. It has nevertheless been shown that at most 50% of the EFVs can be ruled out that way, which is clearly insufficient to scale up to GSMMs. Adding bounds on metabolite concentrations has already proved capable of ruling out a higher percentage of EFVs. This is achieved by checking the extreme rays at each step of the DD algorithm thanks to a call to a linear programming solver and it will be interesting to compare its efficiency with the method given by Proposition 4 which does not need using an LP solver but has the disadvantage of requiring reasoning in a higher dimension. As the structure of the solution space cannot be deduced from merely the knowledge of the EFVs but requires the identification of the largest consistent (w.r.t. the constraints) sets of EFVs, efficient computation of these LCS C ’s will be studied. We can reasonably think that scaling up will be obtained only by dealing with all biological constraints together. However, as it has been shown, handling regulatory constraints poses serious problems as these constraints are not sign-monotone and thus filtering of the EFVs that satisfy them cannot be integrated incrementally into the DD algorithm. In addition the structure of the solution space as union of semi-open polyhedral cones (or semi-open polyhedra) is more complex and pathways of interest for biologists, such as elementary flux modes or support-wise non-strictly-decomposable vectors, no longer coincide with EFVs and require more complex computations to be identified. Algorithms will be carefully studied for maximal efficiency and novel ways of using the DD method or the use of other methods recently proposed such as local reverse search or satisfiability based methods [42,50] will be investigated.

Funding

This research received no external funding.

Data Availability Statement

Data sharing not applicable.

Acknowledgments

I would like to express my sincere thanks to my colleagues Sabine Peres, who introduced me a few years ago to bioinformatics and more particularly to metabolic pathway analysis, that was a field completely unknown to me, and Antoine Deza, who explained to me all I had to know about convex geometry and the double description algorithm. Without their help and the huge number of hours we spent working and discussing together, this work could never have been accomplished. I also warmly thank Rosemary Patricot for having proofread the English.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Urbanczik, R.; Wagner, C. An improved algorithm for stoichiometric network analysis: Theory and applications. Bioinformatics 2005, 21, 1203–1210. [Google Scholar] [CrossRef] [PubMed]
  2. Wagner, C.; Urbanczik, R. The Geometry of the Flux Cone of a Metabolic Network. Biophys. J. 2005, 89, 3837–3845. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Trinh, C.T.; Wlaschin, A.; Srienc, F. Elementary mode analysis: A useful metabolic pathway analysis tool for characterizing cellular metabolism. Appl. Microbiol. Biotechnol. 2009, 81, 813–826. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Llaneras, F.; Picó, J. Which Metabolic Pathways Generate and Characterize the Flux Space? A Comparison among Elementary Modes, Extreme Pathways and Minimal Generators. J. Biomed. Biotechnol. 2010, 2010, 753904. [Google Scholar] [CrossRef] [Green Version]
  5. Zanghellini, J.; Ruckerbauer, D.E.; Hanscho, M.; Jungreuthmayer, C. Elementary flux modes in a nutshell: Properties, calculation and applications. Biotechnol. J. 2013, 8, 1009–1016. [Google Scholar] [CrossRef]
  6. Müller, S.; Regensburger, G. Elementary Vectors and Conformal Sums in Polyhedral Geometry and their Relevance for Metabolic Pathway Analysis. Front. Genet. 2016, 7, 90. [Google Scholar] [CrossRef] [Green Version]
  7. Klamt, S.; Regensburger, G.; Gerstl, M.P.; Jungreuthmayer, C.; Schuster, S.; Mahadevan, R.; Zanghellini, J.; Müller, S. From elementary flux modes to elementary flux vectors: Metabolic pathway analysis with arbitrary linear flux constraints. PLoS Comput. Biol. 2017, 13, e1005409. [Google Scholar] [CrossRef]
  8. Schilling, C.H.; Palsson, B.O. The underlying pathway structure of biochemical reaction networks. Proc. Natl. Acad. Sci. USA 1998, 95, 4193–4198. [Google Scholar] [CrossRef] [Green Version]
  9. Schilling, C.H.; Letscher, D.; Palsson, B.O. Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective. J. Theor. Biol. 2000, 203, 229–248. [Google Scholar] [CrossRef] [Green Version]
  10. Jevremovic, D.; Trinh, C.T.; Srienc, F.; Boley, D. On Algebraic Properties of Extreme Pathways in Metabolic Networks. J. Comput. Biol. 2010, 17, 107–119. [Google Scholar] [CrossRef]
  11. Motzkin, T.S.; Raiffa, H.; Thompson, G.L.; Thrall, R.M. The double description method. In Contributions to the Theory of Games II, Annals of Mathematics Studies; Kuhn, H.W., Tucker, A.W., Eds.; Princeton University Press: Princeton, NJ, USA, 1953; Volume 28. [Google Scholar]
  12. Fukuda, K.; Prodon, A. Double Description Method Revisited. In Combinatorics and Computer Science; Deza, M., Euler, R., Manoussakis, I., Eds.; Springer: Berlin/Heidelberg, Germany, 1996; Volume 1120, pp. 91–111. [Google Scholar]
  13. Gerstl, M.P.; Müller, S.; Regensburger, G.; Zanghellini, J. Flux tope analysis: Studying the coordination of reaction directions in metabolic networks. Bioinformatics 2019, 35, 266–273. [Google Scholar] [CrossRef]
  14. Gudmundsson, S.; Thiele, I. Computationally efficient flux variability analysis. BMC Bioinform. 2010, 11, 489. [Google Scholar] [CrossRef] [Green Version]
  15. Clarke, B.L. Stoichiometry network analysis. Cell Biophys. 1988, 12, 237–253. [Google Scholar] [CrossRef]
  16. Rockafellar, R.T. The elementary vectors of a subspace of RN. In Combinatorial Mathematics and Its Applications; University of North Carolina Press: Chapel Hill, NC, USA, 1969; pp. 104–127. [Google Scholar]
  17. Schuster, S.; Hilgetag, C. On elementary flux modes in biochemical reaction systems at steady state. J. Biol. Syst. 1994, 2, 165–182. [Google Scholar] [CrossRef]
  18. Gagneur, J.; Klamt, S. Computation of elementary modes: A unifying framework and the new binary approach. BMC Bioinform. 2004, 5, 175. [Google Scholar] [CrossRef] [Green Version]
  19. de Figueiredo, L.F.; Podhorski, A.; Rubio, A.; Kaleta, C.; Beasley, J.E.; Schuster, S.; Planes, F.J. Computing the shortest elementary flux modes in genome-scale metabolic networks. Bioinformatics 2009, 25, 3158–3165. [Google Scholar] [CrossRef]
  20. Khachiyan, L.; Boros, E.; Borys, K.; Elbassioni, K.; Gurvich, V. Generating all vertices of a polyhedron is hard. Discret. Comput. Geom. 2008, 39, 174–190. [Google Scholar] [CrossRef]
  21. Acuña, V.; Marchetti-Spaccamela, A.; Sagot, M.F.; Stougie, L. A note on the complexity of finding and enumerating elementary modes. BioSystems 2010, 99, 210–214. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Acuña, V.; Chierichetti, F.; Lacroix, V.; Marchetti-Spaccamela, A.; Sagot, M.F.; Stougie, L. Modes and cuts in metabolic networks: Complexity and algorithms. BioSystems 2009, 95, 51–60. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Dyer, M.E. The complexity of vertex enumeration methods. Math. Oper. Res. 1983, 8, 381–402. [Google Scholar] [CrossRef]
  24. McMullen, P. The maximum numbers of faces of a convex polytope. Mathematika 1970, 17, 179–184. [Google Scholar] [CrossRef] [Green Version]
  25. Terzer, M.; Stelling, J. Large-scale computation of elementary flux modes with bit pattern trees. Bioinformatics 2008, 24, 2229–2235. [Google Scholar] [CrossRef]
  26. van Klinken, J.B.; Willems van Dijk, K. FluxModeCalculator: An efficient tool for large-scale flux mode computation. Bioinformatics 2016, 32, 1265–1266. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Atkins, P.; de Paula, J. Physical Chemistry, 10th ed.; Freeman: Alexandria, VA, USA, 2014. [Google Scholar]
  28. Henry, C.S.; Broadbelt, L.J.; Hatzimanikatis, V. Thermodynamics-based metabolic flux analysis. Biophys. J. 2007, 92, 1792–1805. [Google Scholar] [CrossRef] [Green Version]
  29. Liebermeister, W.; Uhlendorf, J.; Klipp, E. Modular rate laws for enzymatic reactions: Thermodynamics, elasticities and implementation. Bioinformatics 2010, 26, 1528–1534. [Google Scholar] [CrossRef]
  30. Noor, E.; Flamholz, A.; Liebermeister, W.; Bar-Even, A.; Milo, R. A note on the kinetics of enzyme action: A decomposition that highlights thermodynamic effects. FEBS Lett. 2013, 587, 2772–2777. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Jungreuthmayer, C.; Ruckerbauer, D.E.; Zanghellini, J. Utilizing gene regulatory information to speed up the calculation of elementary flux modes. arXiv 2012, arXiv:1208.1853. [Google Scholar]
  32. Jungreuthmayer, C.; Ruckerbauer, D.E.; Zanghellini, J. regEfmtool: Speeding up elementary flux mode calculation using transcriptional regulatory rules in the form of three-state logic. BioSystems 2013, 113, 37–39. [Google Scholar] [CrossRef] [PubMed]
  33. Jungreuthmayer, C.; Ruckerbauer, D.E.; Gerstl, M.P.; Hanscho, M.; Zanghellini, J. Avoiding the enumeration of infeasible elementary flux modes by including transcriptional regulatory rules in the enumeration process saves computational costs. PLoS ONE 2015, 10, e0129840. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Burgard, A.P.; Nikolaev, E.V.; Schilling, C.H.; Maranas, C.D. Flux coupling analysis of genome-scale metabolic network reconstructions. Genome Res. 2004, 14, 301–312. [Google Scholar] [CrossRef] [Green Version]
  35. Larhlimi, A.; David, L.; Selbig, J.; Bockmayr, A. F2C2: A fast tool for the computation of flux coupling in genome-scale metabolic networks. BMC Bioinform. 2012, 13, 57. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Röhl, A.; Goldstein, Y.; Bockmayr, A. EFM–Recorder—Faster Elementary Mode Enumeration via Reaction Coupling Order. Adv. Syst. and Synth. Biol. 2015, 1, 91–99. [Google Scholar]
  37. Peres, S.; Jolicœur, M.; Moulin, C.; Dague, P.; Schuster, S. How important is thermodynamics for identifying elementary flux modes? PLoS ONE 2017, 12, e0171440. [Google Scholar] [CrossRef] [PubMed]
  38. Peres, S.; Schuster, S.; Dague, P. Thermodynamic constraints for identifying the elementary flux modes. Biochem. Soc. Trans. 2018, 46, 641–647. [Google Scholar] [CrossRef] [PubMed]
  39. Gerstl, M.P.; Jungreuthmayer, C.; Muller, S.; Zanghellini, J. Which sets of elementary flux modes form thermodynamically feasible flux distributions? FEBS J. 2016, 283, 1782–1794. [Google Scholar] [CrossRef]
  40. Müller, S.; Regensburger, G.; Steuer, R. Enzyme allocation problems in kinetic metabolic networks: Optimal solutions are elementary flux modes. J. Theor. Biol. 2014, 347, 182–190. [Google Scholar] [CrossRef] [Green Version]
  41. Wortel, M.T.; Peters, H.; Hulshof, J.; Teusink, B.; Bruggeman, F.J. Metabolic states with maximal specific rate carry flux through an elementary flux mode. FEBS J. 2014, 281, 1547–1555. [Google Scholar] [CrossRef] [Green Version]
  42. Morterol, M.; Dague, P.; Peres, S.; Simon, L. Minimality of Metabolic Flux Modes under Boolean Regulation Constraints. In Proceedings of the 12th International Workshop on Constraint-Based Methods for Bioinformatics (WCB’16), Toulouse, France, 5 September 2016. [Google Scholar]
  43. Wagner, C. Nullspace Approach to Determine the Elementary Modes of Chemical Reaction Systems. J. Phys. Chem. B 2004, 108, 2425–2431. [Google Scholar] [CrossRef]
  44. Gerstl, M.P.; Jungreuthmayer, C.; Zanghellini, J. tEFMA: Computing thermodynamically feasible elementary flux modes in metabolic networks. Bioinformatics 2015, 31, 2232–2234. [Google Scholar] [CrossRef] [Green Version]
  45. Gerstl, M.P.; Ruckerbauer, D.E.; Mattanovich, D.; Jungreuthmayer, C.; Zanghellini, J. Metabolomics integrated elementary flux mode analysis in large metabolic networks. Sci. Rep. 2015, 5, 8930. [Google Scholar] [CrossRef] [Green Version]
  46. Haus, U.U.; Klamt, S.; Stephen, T. Computing knock-out strategies in metabolic networks. J. Comput. Biol. 2008, 15, 259–268. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Jungreuthmayer, C.; Nair, G.; Klamt, S.; Zanghellini, J. Comparison and improvement of algorithms for computing minimal cut sets. BMC Bioinform. 2013, 14, 318. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Ballerstein, K.; von Kamp, A.; Klamt, S.; Haus, U. Minimal cut sets in a metabolic network are elementary modes in a dual network. Bioinformatics 2012, 28, 381–387. [Google Scholar] [CrossRef]
  49. von Kamp, A.; Klamt, S. Enumeration of Smallest Intervention Strategies in Genome-Scale Metabolic Networks. PLoS Comput. Biol. 2014, 10, e1003378. [Google Scholar] [CrossRef] [PubMed]
  50. Mahout, M.; Carlson, R.P.; Peres, S. Answer Set Programming for Computing Constraints-Based Elementary Flux Modes: Application to Escherichia coli Core Metabolism. Processes 2020, 8, 1649. [Google Scholar] [CrossRef]
Figure 1. On the left, we consider the convex polyhedral cone shaded in gray. The extreme vectors are the blue vectors. The elementary vectors are all the blue or red vectors and coincide with the extreme vectors for all topes (here we have two topes: the subcone with x 0 and the subcone with x 0 ). The elementary mode is the red vector (common to both topes). Flux cones being defined by linear inequalities of the form v i 0 , elementary modes coincide with elementary vectors in this case (at right).
Figure 1. On the left, we consider the convex polyhedral cone shaded in gray. The extreme vectors are the blue vectors. The elementary vectors are all the blue or red vectors and coincide with the extreme vectors for all topes (here we have two topes: the subcone with x 0 and the subcone with x 0 ). The elementary mode is the red vector (common to both topes). Flux cones being defined by linear inequalities of the form v i 0 , elementary modes coincide with elementary vectors in this case (at right).
Computation 09 00111 g001
Figure 2. We consider the convex polyhedron in gray and its associated recession cone delimited by dashed lines. The extreme points are the vertices in blue and the extreme vectors are the blue vectors. The elementary points are all the blue or red points, and the elementary vectors are all the blue or red vectors, and they coincide, respectively, with the extreme points and extreme vectors for all topes (here we have four topes, corresponding to the subpolyhedra in the four quadrants defined by x and y). The elementary modes are the purple segment and half-line and they do not have a direct relationship with elementary elements, but their two possible supports, { x } and { y } , are the minimal nonempty supports of the elementary elements.
Figure 2. We consider the convex polyhedron in gray and its associated recession cone delimited by dashed lines. The extreme points are the vertices in blue and the extreme vectors are the blue vectors. The elementary points are all the blue or red points, and the elementary vectors are all the blue or red vectors, and they coincide, respectively, with the extreme points and extreme vectors for all topes (here we have four topes, corresponding to the subpolyhedra in the four quadrants defined by x and y). The elementary modes are the purple segment and half-line and they do not have a direct relationship with elementary elements, but their two possible supports, { x } and { y } , are the minimal nonempty supports of the elementary elements.
Computation 09 00111 g002
Figure 3. Structure of the solution space in the presence of the thermodynamic constraint (without bounds on the metabolite concentrations). We consider a given flux tope and represent the section of the flux cone with an affine hyperplane (polytope shaded in pink). Vertices e 1 to e 8 thus represent the EFMs. When we add the thermodynamic constraint, the solution space is the union of all the maximal faces of the flux cone that are entirely contained in the open half-space { z R r z T ln K ^ eq > 0 } (left side of frontier H in the section), represented in red. In particular, the thermodynamically feasible EFMs are given by the vertices e 1 to e 5 that belong to this half-space. The set of these EFMs is decomposed into LTCSs, made up of the EFMs of the maximal faces above (in the 2D shape, these LTCSs are thus given by { { e 1 , e 2 } , { e 2 , e 3 } , { e 3 , e 4 } , { e 4 , e 5 } } ).
Figure 3. Structure of the solution space in the presence of the thermodynamic constraint (without bounds on the metabolite concentrations). We consider a given flux tope and represent the section of the flux cone with an affine hyperplane (polytope shaded in pink). Vertices e 1 to e 8 thus represent the EFMs. When we add the thermodynamic constraint, the solution space is the union of all the maximal faces of the flux cone that are entirely contained in the open half-space { z R r z T ln K ^ eq > 0 } (left side of frontier H in the section), represented in red. In particular, the thermodynamically feasible EFMs are given by the vertices e 1 to e 5 that belong to this half-space. The set of these EFMs is decomposed into LTCSs, made up of the EFMs of the maximal faces above (in the 2D shape, these LTCSs are thus given by { { e 1 , e 2 } , { e 2 , e 3 } , { e 3 , e 4 } , { e 4 , e 5 } } ).
Computation 09 00111 g003
Figure 4. Simple network of Example 1 with R and T 1 reversible. e 1 , e 4 , e 3 are EFMs. When we add the regulatory constraint T 1 T 4 , e 1 is still an EFM, which is not the case of any positive combination of e 3 and e 4 , which has a larger support. Nevertheless, the latter is an EFM in the flux tope defined by reverse fluxes in R and T 1 .
Figure 4. Simple network of Example 1 with R and T 1 reversible. e 1 , e 4 , e 3 are EFMs. When we add the regulatory constraint T 1 T 4 , e 1 is still an EFM, which is not the case of any positive combination of e 3 and e 4 , which has a larger support. Nevertheless, the latter is an EFM in the flux tope defined by reverse fluxes in R and T 1 .
Computation 09 00111 g004
Figure 5. Simple network of Example 2. When we add the regulatory constraint ¬ R 1 R 2 , e 1 is the only EFM, and also the only one for the disjunct D 1 = ¬ R 1 of the constraint, but no e 2 combination is an EFM, while it is an EFM for the disjunct D 2 = R 1 R 2 .
Figure 5. Simple network of Example 2. When we add the regulatory constraint ¬ R 1 R 2 , e 1 is the only EFM, and also the only one for the disjunct D 1 = ¬ R 1 of the constraint, but no e 2 combination is an EFM, while it is an EFM for the disjunct D 2 = R 1 R 2 .
Computation 09 00111 g005
Figure 6. Simple network of Example 3. e 2 and e 3 are EFMs. When we add the regulatory constraint T 1 T 4 , any positive combination of e 2 and e 3 is now an EFM, that is not obtained as an original EFM satisfying the constraint.
Figure 6. Simple network of Example 3. e 2 and e 3 are EFMs. When we add the regulatory constraint T 1 T 4 , any positive combination of e 2 and e 3 is now an EFM, that is not obtained as an original EFM satisfying the constraint.
Computation 09 00111 g006
Figure 7. Second network of Example 4. e 1 to e 8 are all the EFMs. When we add the regulatory constraint T 1 T 2 , only e 1 remains an EFM, with support { T 1 , T 2 , R } , but new EFMs appear, made up of any positive combination of e 5 and e 7 , all with the same support { T 1 , T 2 , T 5 , T 6 } .
Figure 7. Second network of Example 4. e 1 to e 8 are all the EFMs. When we add the regulatory constraint T 1 T 2 , only e 1 remains an EFM, with support { T 1 , T 2 , R } , but new EFMs appear, made up of any positive combination of e 5 and e 7 , all with the same support { T 1 , T 2 , T 5 , T 6 } .
Computation 09 00111 g007
Figure 8. Third network of Example 4. e 1 to e 6 are all the EFMs. When we add the regulatory constraint R 1 R 2 T 2 , only e 1 remains an EFM, with support { T 1 , R 1 , R 2 , T 2 } , but new EFMs appear, made up of any positive combination of e 3 and e 5 , all with the same support { T 4 , T 5 , R 1 , R 2 , T 2 } .
Figure 8. Third network of Example 4. e 1 to e 6 are all the EFMs. When we add the regulatory constraint R 1 R 2 T 2 , only e 1 remains an EFM, with support { T 1 , R 1 , R 2 , T 2 } , but new EFMs appear, made up of any positive combination of e 3 and e 5 , all with the same support { T 4 , T 5 , R 1 , R 2 , T 2 } .
Computation 09 00111 g008
Figure 9. Simple network of Examples 3 and 5. The EFMs are e 1 , e 2 , e 3 . When we add the regulatory constraint T 1 T 4 , only e 1 remains an EFM, with support { T 1 , R , T 4 } , but new EFMs appear, made up of any positive combination of e 2 and e 3 , all with the same support { T 1 , T 2 , T 3 , T 4 } . Moreover, new swNSDFVs also appear, which are not EFMs: any positive combination of e 1 and e 2 , all with the same support { T 1 , R , T 3 , T 4 } , and any positive combination of e 1 and e 3 , all with the same support { T 1 , R , T 2 , T 4 } .
Figure 9. Simple network of Examples 3 and 5. The EFMs are e 1 , e 2 , e 3 . When we add the regulatory constraint T 1 T 4 , only e 1 remains an EFM, with support { T 1 , R , T 4 } , but new EFMs appear, made up of any positive combination of e 2 and e 3 , all with the same support { T 1 , T 2 , T 3 , T 4 } . Moreover, new swNSDFVs also appear, which are not EFMs: any positive combination of e 1 and e 2 , all with the same support { T 1 , R , T 3 , T 4 } , and any positive combination of e 1 and e 3 , all with the same support { T 1 , R , T 2 , T 4 } .
Computation 09 00111 g009
Figure 10. Simple network of Examples 3, 5 and 6. The flux cone F C of dimension 3 has 3 extreme rays, which are EFMs e 1 , e 2 , e 3 , and 3 facets (in color in the picture that represents the projection of R 5 on the subspace spanned by x , y , z , the fluxes in R , T 3 , T 2 , respectively). When we add the regulatory constraint T 1 T 4 , the solution space becomes F C without its two edges corresponding to e 2 and e 3 (in gray), thus a disjoint union of the following open cones: extreme ray e 1 , the interiors of the red, blue and green facets, and the interior of F C . Thus, only e 1 remains an EFM, with support { T 1 , R , T 4 } , but there are new EFMs, namely all positive combinations of e 2 and e 3 (the interior of the red facet), all with the same minimal support { T 1 , T 3 , T 2 , T 4 } . Moreover, there are new swNSDFVs with non-minimal support: all positive combinations of e 1 and e 3 (the interior of the blue facet), all with the same non-strictly-decomposable support { T 1 , R , T 2 , T 4 } , and all positive combinations of e 1 and e 2 (the interior of the green facet), all with the same non-strictly-decomposable support { T 1 , R , T 3 , T 4 } . All vectors in the interior of F C , which have { T 1 , R , T 3 , T 2 , T 4 } as a support, are support-wise strictly-decomposable.
Figure 10. Simple network of Examples 3, 5 and 6. The flux cone F C of dimension 3 has 3 extreme rays, which are EFMs e 1 , e 2 , e 3 , and 3 facets (in color in the picture that represents the projection of R 5 on the subspace spanned by x , y , z , the fluxes in R , T 3 , T 2 , respectively). When we add the regulatory constraint T 1 T 4 , the solution space becomes F C without its two edges corresponding to e 2 and e 3 (in gray), thus a disjoint union of the following open cones: extreme ray e 1 , the interiors of the red, blue and green facets, and the interior of F C . Thus, only e 1 remains an EFM, with support { T 1 , R , T 4 } , but there are new EFMs, namely all positive combinations of e 2 and e 3 (the interior of the red facet), all with the same minimal support { T 1 , T 3 , T 2 , T 4 } . Moreover, there are new swNSDFVs with non-minimal support: all positive combinations of e 1 and e 3 (the interior of the blue facet), all with the same non-strictly-decomposable support { T 1 , R , T 2 , T 4 } , and all positive combinations of e 1 and e 2 (the interior of the green facet), all with the same non-strictly-decomposable support { T 1 , R , T 3 , T 4 } . All vectors in the interior of F C , which have { T 1 , R , T 3 , T 2 , T 4 } as a support, are support-wise strictly-decomposable.
Computation 09 00111 g010
Figure 11. Simple network of Example 7. Flux cone F C of dimension 3 has 4 extreme rays, which are EFMs e 1 , e 2 , e 3 , e 4 , and 4 facets (in color in the picture that represents the projection of R 5 on the subspace spanned by x , y , z , the fluxes in T 1 , T 2 , T 3 , respectively). When we add the regulatory constraint T 1 T 3 , the solution space becomes F C without its two facets in yellow and orange, thus a disjoint union of the following open cones: the extreme ray e 1 , the interiors of the blue and purple facets, and the interior of F C . Thus, only e 1 remains an EFM, with support { T 1 , R , T 3 } and there are no new EFMs, but the following new swNSDFVs occur. First, any positive combination of e 1 and e 2 (the interior of the blue facet), all with the same support { T 1 , R , T 3 , T 4 } , and any positive combination of e 1 and e 3 (the interior of the purple facet), all with the same support { T 1 , T 2 , R , T 3 } , these supports being thus non-strictly-decomposable (independently of the positive values of the fluxes) though not minimal. Second, the complementary in the interior of F C of the open sub-cone generated by the positive combinations of e 1 , e 2 , e 3 (i.e., the interior of the cone with blue, purple and grey facets), that is the original network with nonzero fluxes in all five reactions, thus with support { T 1 , T 2 , R , T 3 , T 4 } , in the case where input flux y in T 2 is not smaller than output flux z in T 3 , obtained as a positive combination of e 2 , e 3 , e 4 (the interior of the cone with yellow, orange and grey facets) if y > z , or of e 2 , e 3 (the interior of the gray cone) if y = z . On the other hand, when y is smaller than z, the global pathway can be decomposed into two swNSDFVs of non-strictly-decomposable supports { T 1 , R , T 3 , T 4 } and { T 1 , T 2 , R , T 3 } in an infinite number of ways depending on a parameter k , 0 < k < 1 (decomposition of a vector in the interior of the cone with blue, purple and gray facets into two vectors in the interiors of the blue and purple facets respectively). This shows that in general support-wise strict-decomposability does not depend only on the support of the pathway but on the distribution of the fluxes.
Figure 11. Simple network of Example 7. Flux cone F C of dimension 3 has 4 extreme rays, which are EFMs e 1 , e 2 , e 3 , e 4 , and 4 facets (in color in the picture that represents the projection of R 5 on the subspace spanned by x , y , z , the fluxes in T 1 , T 2 , T 3 , respectively). When we add the regulatory constraint T 1 T 3 , the solution space becomes F C without its two facets in yellow and orange, thus a disjoint union of the following open cones: the extreme ray e 1 , the interiors of the blue and purple facets, and the interior of F C . Thus, only e 1 remains an EFM, with support { T 1 , R , T 3 } and there are no new EFMs, but the following new swNSDFVs occur. First, any positive combination of e 1 and e 2 (the interior of the blue facet), all with the same support { T 1 , R , T 3 , T 4 } , and any positive combination of e 1 and e 3 (the interior of the purple facet), all with the same support { T 1 , T 2 , R , T 3 } , these supports being thus non-strictly-decomposable (independently of the positive values of the fluxes) though not minimal. Second, the complementary in the interior of F C of the open sub-cone generated by the positive combinations of e 1 , e 2 , e 3 (i.e., the interior of the cone with blue, purple and grey facets), that is the original network with nonzero fluxes in all five reactions, thus with support { T 1 , T 2 , R , T 3 , T 4 } , in the case where input flux y in T 2 is not smaller than output flux z in T 3 , obtained as a positive combination of e 2 , e 3 , e 4 (the interior of the cone with yellow, orange and grey facets) if y > z , or of e 2 , e 3 (the interior of the gray cone) if y = z . On the other hand, when y is smaller than z, the global pathway can be decomposed into two swNSDFVs of non-strictly-decomposable supports { T 1 , R , T 3 , T 4 } and { T 1 , T 2 , R , T 3 } in an infinite number of ways depending on a parameter k , 0 < k < 1 (decomposition of a vector in the interior of the cone with blue, purple and gray facets into two vectors in the interiors of the blue and purple facets respectively). This shows that in general support-wise strict-decomposability does not depend only on the support of the pathway but on the distribution of the fluxes.
Computation 09 00111 g011
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Dague, P. Metabolic Pathway Analysis in the Presence of Biological Constraints. Computation 2021, 9, 111. https://doi.org/10.3390/computation9100111

AMA Style

Dague P. Metabolic Pathway Analysis in the Presence of Biological Constraints. Computation. 2021; 9(10):111. https://doi.org/10.3390/computation9100111

Chicago/Turabian Style

Dague, Philippe. 2021. "Metabolic Pathway Analysis in the Presence of Biological Constraints" Computation 9, no. 10: 111. https://doi.org/10.3390/computation9100111

APA Style

Dague, P. (2021). Metabolic Pathway Analysis in the Presence of Biological Constraints. Computation, 9(10), 111. https://doi.org/10.3390/computation9100111

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop