Metabolic pathway analysis in the presence of biological constraints

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, genomescale, 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 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 the sign function: either signinvariance for regulatory constraints or sign-monotonicity (a stronger property) for thermodynamic and kinetic constraints. We show in particular that EFs for sign-monotone constraints 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 EFMs computation.


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 [42,46,41,25,48,30,22], in addition to the references in the text, for the proofs or historical 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 ji 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: dM(t) dt = Sv. (1)

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: With only this assumption, the solution space Sol 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 [38], of S: The dimension of the flux subspace is given by dim(F S) = r − rank(S) ≥ r − m. Often, possibly after a pre-processing to eliminate its linearly dependent rows, S is assumed to be of full rank and then dim(F S) = r − m.
The support of a vector x ∈ R r is defined by: 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).

Irreversible reactions, flux cones and polyhedral cones
In addition to the homogeneous equality constraints provided by the steady-state assumption, the flux vectors have also in general to satisfy homogeneous inequality constraints corresponding to reactions known as irreversible, whose set will be noted Irr: This means that fluxes in irreversible reactions are constrained to be nonnegative (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 Sol S,Irr is the intersection of the linear subspace F S with r I nonnegative 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: A (general) convex polyhedral cone is defined implicitly (or by intension) by finitely many homogeneous linear inequalities: 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 dim(C), is defined as the dimension of its affine span. A flux cone corresponds thus to a particular matrix A ∈ R (2m+r I )×r given by: 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 nonnegativity 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: 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 | Ax = 0}. For a flux cone F C given by (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.

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 [37,17] 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., nonnegative linear) combination: where the y k 's, k ∈ K (finite index set), are representatives of the extreme rays (unique up to positive scalar multiplication) and cone 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∈Kx β k y k with |K x | ≤ min(dim(C), |supp(x)|+|supp(Ax)|). 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 | ≤ min(dim(C), |supp(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 [29,8], known as Fourier-Motzkin, is an incremental algorithm (which processes one by one each linear inequality (Ax) 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: 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's 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).

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 cancelations is biochemically meaningful, since 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 cancelations, are biochemically admissible. A sum v = v 1 + v 2 of vectors is called conformal if, for all i ∈ {1, . . . , r}: An equivalent definition is to define a sum v = v 1 + v 2 as conformal if: where the sign vector sign(v) ∈ {−, 0, +} r is defined by applying the sign function component-wise, i.e., sign(v) i = sign(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 | sign(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 sign(v) ≤ η. We call two vectors v 1 , v 2 conformal if sign(v 1 ), sign(v 2 ) ≤ η for a certain sign vector η or, equivalently, if v 1 i v 2 i ≥ 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 (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. (15) A vector x (resp., flux vector v) of a convex polyhedral cone C (resp., a flux cone F C) is called elementary [22] 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: where the e g 's, g ∈ V G (finite index set), are representatives of the elementary rays (unique up to positive scalar multiplication) and cone ⊕ 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 Gx β g e g with |V G x | ≤ min(dim(C), |supp(x)| + |supp(Ax)|). 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 | ≤ min(dim(C), |supp(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 since 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: 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 (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 sign(C) = {sign(x) | x ∈ C}.
It is then enough in (17) to consider the closed orthants O i = {x ∈ R r | sign(x) ≤ η i } and there is no equality or inclusion between the The C ≤η i 's are called topes, noted Ts (flux topes [12] noted FTs for a flux cone F C). C is thus decomposed into topes and (17) can be rewritten as: ExVs(C ≤η i ).
Note that, for a flux cone F C (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 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 , wherer = 2r − r I is the new number of reactions and allr 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 [6]). F C is included in the positiver-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 ≤η .

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 (13), (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: 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 [35], for a subspace of R r , i.e., for a flux linear subspace F S (3) and then [39] 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 sign(x) is a minimal (for the partial order defined above on {−, 0, +} r ) nonzero element of sign(F C) and, conversely, it is shown that a minimal nonzero sign vector σ ∈ sign(F C) determines an EFM x with sign(x) = σ by: There is thus a one-to-one mapping between EFMs and minimal nonzero sign vectors of sign(F C). Comparing with the one-to-one mapping between FTs and maximal sign vectors of sign(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 (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 (16) for flux cones F C is to prove it with EFMs as a conformal generating set and to prove that EFMs = EFVs). From (18), for any maximal sign vector η of sign(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 [9,7]. EFMs are thus decomposed into subsets according to the decomposition of F C into flux topes [12]: the EFMs of the FT F C ≤η correspond to the F C ≤σ 's, where the σ's are the minimal nonzero sign vectors of sign(F C) such that σ ≤ η.
For a general polyhedral cone C, there is no direct relationship between elementary modes and elementary vectors but 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 : And 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 supp(x ) ⊂ supp(x), then a vector x exists in this convex set with sign(x ) ≤ sign(x) and supp(x ) ⊂ supp(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:

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: or, more generally, any set of inhomogeneous linear constraints, noted ILC, that can be written in the general form: 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 | Gv ≥ h}. The solution space Sol S,Irr,ILC is thus now a s-polyhedron or flux polyhedron noted F P and defined by: This is a particular case of (general) convex polyhedron that is defined implicitly by finitely many linear inequalities: with a suitable matrix A ∈ R n×r and vector b ∈ R n and is thus the intersection of n (affine) half-spaces. 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 (2m+r I +l)×r and vector b ∈ R 2m+r I +l given by: meaning that the inequalities that are homogeneous actually divide into m equalities defining a lowerdimensional subspace given by the nullspace of the stoichiometric matrix S and into r I nonnegativity 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 (25) is associated its so-called recession cone C P = {x ∈ R r | Ax ≥ 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 | Ax = 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.

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 distinct vectors of P : 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", resp.) generators of P for (convex and conical, resp.) combination: 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 conv is the convex hull (if P is a pointed polyhedral cone, then it has only one vertex, which is the zero vector, and the formula (28) reduces to (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 : . 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 (25), we can associate a flux cone F C P in a higher dimension R r+1+n defined by: This introduces a correspondence between vectors x of P and vectors 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 Ax 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 (24), the correspondence (29) simplifies by associating to F P the flux cone F C F P in dimension R r+1+l defined by: 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 (25), i.e., enumerates its vertices ExPs and extreme vectors ExVs.
If P is not pointed, it is still finitely generated: 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 formula (31) reduces to (12)). In fact, Minkowski-Weyl's theorem for polyhedra states that it is equivalent for a set P to be a polyhedron (25) 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).

Elementary points and vectors and conformal generating sets
As was the case for a flux cone, however, such decomposition into a finite set of generators is not in general satisfactory for a flux polyhedron as only a decomposition without cancelations is biochemically meaningful. In the same way 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: 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) [22]. 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: 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 conv ⊕ 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 Gx α g e g ⊕ g∈V Gx β g e g with . This result can be demonstrated from result (16) for a flux cone by using the correspondence (29) 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 (27) and (32) that an extreme point (vertex) of P is an elementary point of P , i.e., ExPs ⊆ EPs, but the converse is generally false since 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): 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: If O i is the closed orthant defined by η i , then the corresponding tope for P is P ∩ O i = P ≤η i = {x ∈ P | sign(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 ) = max{η ∈ sign(C P ) | η ≤ η i } is the greatest (it is necessarily unique) sign vector in sign(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 (24), 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 , wherer = 2r − r I is the new number of reactions and allr 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 positiver-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. And 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 ≤η .

Elementary modes
For a general polyhedron P (and even for a flux polyhedron F P ), there is no direct relationship between support-minimal (19) vectors and elementary elements as it was the case for flux cones F C. Nevertheless, from (33), 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, all minimal support patterns of vectors appear in the set of supports of elementary elements and are actually the minimal elements in this set for subset inclusion: 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 (29) 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): with F C P given by (29). The same formula holds for the particular case of a flux polyhedron F P (24) with F C F P given by (30).
And, as seen in the proof of (21), 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:

Complexity results
Enumerating the extreme vectors of a polyhedral cone or the vertices of a polytope has been proven to be a #P-complete problem. A consequence is that enumerating EFs or EFMs of a metabolic network is #P-complete too [1,2]. 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) is an open problem. So the possibility of enumerating EFs or EFMs of a metabolic network with a polynomial delay is an open problem (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 [35]). 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 [1,2]. The McMullen's upper bound theorem [27] 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 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 = 2r − r I (sor varies between r and 2r). If the number m of internal metabolites is small compared to the total numberr 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 tor, this number is bounded above by a quantity of order Θ(r (r−m)/2 ). The worst case occurs when m is close tor/3 with an upper bound approximately equal to 2 2m m . We obtain the same results with r instead ofr 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 ).

Metabolic pathways in the presence of biological constraints
Although the current improved implementations of the DD method [40,43] 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.

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 [3], given by: where: ∆ r G 0 i 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 ji 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: sign(v i ) ≤ −sign(∆ r G i ). And that a flux vector v is thermodynamically feasible [16] if and only if all its components v i satisfy such a constraint (it is enough to consider those i ∈ supp(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: As at equilibrium ∆ r G i is null, we obtain: Sji jeq is the equilibrium constant of reaction i. Thus, ∆ r G i can be rewritten as ∆ r G i = RT ln( j M Sji j /K i eq ) and the thermodynamic constraint TC M (v) as: 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 = RT ln( j M Sji j /K i eq ) can be rewritten, by dividing the numerator and denominator of the fraction by the terms dealing with external metabolites, as is the apparent equilibrium constant of the reaction i and M (resp., M) the vector of internal (resp., external) metabolite concentrations. The thermodynamic constraint (41) can thus be rewritten as: 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: Then the thermodynamic constraint TC M (resp., TC M ) can be rewritten as: Thus the set Sol TC M (resp., Sol TC M ) of vectors v satisfying the constraint TC M (v) (resp., TC M (v)), given by {v ∈ R r | sign(v) ≤ ts M } (resp., {v ∈ R r | sign(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 2.1. Given metabolite concentrations M (resp., internal metabolite concentrations M), the set Sol TCM (resp., Sol 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 ) (43), i.e., the closed orthant O tsM (resp., O tsM ) defined by ts M (resp., ts M ).

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 [24]. What is important for our study is that these modular rate laws share the general form: where E i is the (nonnegative) 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 reg i and D reg i , 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: where • is the component-wise product of vectors: 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 [32]. 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: where k + i and k − i are the maximal forward and backward rates in reaction i per unit of enzyme and the K M ij 's are the Michaelis constants. Equating the numerator to zero at equilibrium, we obtain: This gives: 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 ∆rGi/RT 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 sign(κ(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:

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 [20,21,19]. 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 Bc, to a flux vector v, the positive literal i is interpreted as v i = 0, i.e., i ∈ supp(v) (4), and the negative literal ¬i as v i = 0, i.e., i / ∈ supp(v). The regulatory constraint RC Bc for v induced by Bc can thus be defined as: Note that coupled reactions as used by Flux Coupling Analysis (FCA) [5,23,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 Bc = i ∨ ¬j, and i partially coupled to j, meaning that zero flux through i is equivalent to zero flux through j, is expressed by Bc = (i ∧ j) ∨ (¬i ∧ ¬j).
By rewriting the Boolean constraint Bc in DNF (Disjunctive Normal Form), Bc = k D k , the set Sol RC Bc = {v ∈ R r | RC Bc (v)} of vectors v satisfying the constraint RC Bc 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 . Thus, the solution space for a disjunct is itself the disjoint union of subspaces each one defined by constraints of type We can iteratively process like this by considering each time the open orthants in Sol RC Bc 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 orthantO rsj , 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 Sol RC Bc 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 η, Sol 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.

Characterizing the solution space
We start from a flux cone (6) (or more generally a flux polyhedron (24)) 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 Sol C , defined as: For given values of quantities x, the solution space is thus the constrained flux cone subset (not necessarily a cone, depending on C x ): or, more generally, the constrained flux polyhedron subset (not necessarily a polyhedron, depending on C x ): 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 ∃xC x whose solution space The solution space becomes: 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 ∃xC 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 convexconformally 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.: For EFMs, the above formulas are valid whatever the structure of Sol C . It is also the case for ExPs and EFPs but, in practice, only really meaningful if CF P C (x) (resp., CF P C ) is a convex polyhedron. Finally, for ExVs and EFVs, this is only meaningful if CF C C (x) (resp., CF C C ) is a convex polyhedral cone and CF P C (x) (resp., CF P C ) is a convex polyhedron. This is the case when Sol Cx (resp., Sol ∃xCx ) 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 CF P C = CF C C ∩C P ILC = C F P ∩Sol C . This can be generalized when Sol C is a convex polyhedron, and thus CF C C and CF P C too, with C CF C C = F C ∩ C Sol C and C CF P C = C CF C C ∩ C P ILC = C F P ∩ C Sol C , by replacing Sol C by C Sol 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.

Application to thermodynamics
Assume first that the concentrations of the metabolites (resp., of the internal metabolites) are known and given. From (50), (51) and Lemma 2.1, we obtain:  (18), (35) and (38), it follows that: EFMs(CF P TC (M)) = EFMs(F P ) ≤ts M (58) and the same for CF C TC (M) and CF 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 CF C TC (M) (resp., CF P TC (M)) is consistent (η being full support does not imply that η is full support). Let's now consider the more usual case where the concentrations of the metabolites (at least those of internal metabolites) are unknown. From (44), we obtain one or the other form for the thermodynamic constraint (existentially quantified on metabolite concentrations): 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: The solution space in R r of these constraints is thus . From the result of Lemma 2.1, it follows: . 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 2.4. The set Sol TC (resp., Sol 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 tsM 's. The same result holds for TC and TC b with M and O tsM .
It follows from (52), (53), (55), (56) that:  (18), (35) and (38), we get: 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., } 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.5. 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 [33]. (43). By applying the monotonic logarithm function and noting LM ∈ R m the vector whose components are given by LM j = ln(M j ), we obtain: ). 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. For any A ∈ R p×q and b ∈ R p , exactly one of the following statements holds: Applying this theorem in the formula above to can be rewritten as: is an open half-space with the hyperplane H lnKeq = {z ∈ R r | z T lnK eq = 0} as a frontier. Considering the decomposition of F C into FTs F C ≤η , this is equivalent to: lnKeq } for all η's maximal sign vectors in sign(F C). Now, note that, in this formula, F C ≤sign(v is the minimal face (for set inclusion) of F C ≤η containing v. We finally obtain finally: That is to say CF C TC is the union of all the maximal faces F C ≤η i of the FTs for F C such that A topological consequence of this characterization is that CF C TC \{0} is a connected set (actually arc-connected). Another consequence, regarding EFMs (or, equivalently, EFVs), is: 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 + lnKeq . We can equivalently write: where F C + lnKeq is the polyhedral cone obtained from F C merely by adding the homogeneous linear inequality v T lnK eq ≥ 0 [34]. Each flux vector of CF C TC is a conformal conical sum of these EFMs, but the converse is false as CF C TC is not convex. Thus the set of EFMs, considered globally, does not characterize the solution space CF C TC . More precisely, we have: i.e., CF 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 CF C TC , i.e., all vectors in this hull must satisfy the constraint TC: We could want to estimate the ratio of thermodynamically feasible EFMs on all EFMs, i.e., the ratio of EFMs(CF 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 lnKeq 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. So, 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 CF P TC = CF C TC ∩ P ILC , with P ILC = {v ∈ R r | Gv ≥ h} (24) and we obtain in the same way, for all η's maximal sign vectors in sign(F C): That is to say CF 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 + lnKeq 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 + lnKeq . Take care because the F P ≤η i 's involved are faces of F P ≤η included in H + lnKeq ∪ {0}, but not any face of F P ≤η included in H + lnKeq ∪ {0} is included in CF 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). CF P TC is no longer a connected set in general.
We can sum up these results as follows.
Proposition 2.6. 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 + lnKeq = {z | z T lnK 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 lnK eq ≥ 0 (and removing those EFMs that would belong to the frontier hyperplane of H + lnKeq ). 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 Gv ≥ 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 Gv ≥ h. All these results hold for constraint TC by just replacing theK i eq 's by the K i eq '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 + lnKeq . 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 From what precedes, we obtain: This can be rewritten as: is a face of the cone F C b ≤η+ (we note η+ the sign vector of dimension r+2m obtained by concatenation of η of dimension r and + of dimension 2m, thus F C b , 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 ≤sign(v) , this is the maximal one (without any constraint on the R 2m + factor). Thus finally: Note that any F C ≤η j given by (70) is a face of a certain F C ≤η i given by (65) (i.e., ∀η j ∃η i η j ≤ η i ). CF C TC b \{0} is no longer a connected set in general. A consequence, for EFMs (or, equivalently, EFVs), is: i.e., 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 CF 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 Each flux vector of CF C TC b is a conformal conical sum of these EFMs, but the converse is false as CF C TC b is not convex. More precisely, we have: i.e., CF 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 CF C TC b , i.e., all vectors in this hull must satisfy the constraint Such an E j is called a largest thermodynamically consistent (for bounded concentrations of internal metabolites) set (LTCbS) of EFMs [10] 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 (68). If inhomogeneous linear constraints ILC are added, we have, for all η's maximal sign vectors in sign(F C): That is to say CF 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 (70). We can sum up these results as follows.
Proposition 2.7. For F C a flux cone in R r , let's define its lift to R r+2m as the flux cone For any flux tope F C ≤η for F C and any face F C of F C ≤η , its lift to R r+2m 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+2m are entirely contained (except vectors from The thermodynamically feasible EFMs in F C ≤η are those EFMs in F C ≤η whose lifts (as rays) are contained (except vectors from {0} × R 2m ) in this half-space. They 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 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 Gv ≥ h, the space of flux vectors in FP satisfying TC b is a finite union of flux polyhedra, obtained as intersections of the flux cones above with the polyhedron defined by Gv ≥ h.

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: 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: We can also consider intermediate constraints, existentially quantified only on metabolite concentrations if enzyme concentrations are known: possibly with bounds on metabolite concentrations in the quantification: or only on enzyme concentrations if metabolite concentrations are known: possibly with enzymatic resource constraint and bounds on enzyme concentrations in the quantification: The solution space in R r of the constraint where sup (resp., inf ) applies to those M such that sign(κ(M)) = ts, which gives a polyhedron, but also, in the presence of an enzymatic resource constraint, by an algebraic, non-linear, 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 2.8 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 CF C KC = CF C TC , or CF C KC b = CF 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 CF C KC or CF 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 = r i=1 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 r i=1 (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 CF C KC or CF 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. So, 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 CF C KC or CF C KC b . In the case of CF 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 CF C KC , the maximum might not be attained at finite metabolite concentrations. This is the result already obtained in [31,47] and we followed a similar proof, but relying this time on a precise characterization of the solution space CF C KC or CF C KC b and of the EFMs given by the Proposition 2.8, 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 Bc = k, see next subsection), say {e l } (all components of which are fixed by e l k = 1), and, for each one, by computing the minimum (if it exists) of i∈supp(e l ) (w i e l i )/κ i (M) for M varying in its domain, such that sign(κ i (M)) = sign(e l i ) for all i ∈ supp(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 2.9. Given a kinetic metabolic network with possible bounds on metabolite concentrations, but not on enzyme concentrations, i.e., with solution space CF C KC or CF 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 CF C KC or CF C KC b and E T = r i=1 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 CF C KC b if the problem is solvable), then this solution is reached in particular at an EFM of CF C KC or CF C KC b .

Application to regulatory constraints
From Lemma 2.2, we obtain that CF C RC Bc (resp., CF P RC Bc ) is the disjoint union of F C ∩O's (resp., 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 itC (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 ofC (resp.,P ) in R r .C (resp.,P ) is defined as in (7) (resp., in (25)) 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 asF 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 orthantO ⊂ 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 Gv = h). In all, we obtain that CF C RC Bc ∩ O (resp., CF P RC Bc ∩ 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 F C ≤η for F C (resp., F P ≤η for F P ), CF C RC Bc ≤η (resp., CF P RC Bc ≤η ) is a disjoint union of open faces of F C ≤η (resp., F P ≤η ). , 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 onlyC, resp.,P ) of the open faces of C (resp., P ). We have:C ⊆ C • ⊆ C (resp.,P ⊆ P • ⊆ P ). From this geometrical structure of the solution space, we will now deduce what its EFMs and its EFVs (or EFPs) are. Let's 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 sign(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 CF C C (or constrained flux polyhedron subset CF 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 sign(F C) EFMs(F C ≤η ) (and idem for F P ) and remained true for CF 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. 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).
So, 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., CF C RC Bc ≤η (resp., CF P RC Bc ≤η ), thus we can limit ourselves to sign vectors η that are maximal in sign(CF C RC Bc ) (resp., sign(CF P RC Bc )). We saw that we could rewrite the Boolean constraint as a disjunction of two by two exclusive disjuncts, Bc = k D k , decomposing thus the solution space CF C RC Bc ≤η (resp., CF P RC Bc ≤η ) into the disjoint solution spaces for each disjunct, CF C RC D k ≤η (resp., CF 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 CF C RC D k ≤η and the reciprocal is also true: if a flux vector of the semi-open polyhedral cone CF 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 CF C RC D k ≤η of dimension at least two and is thus conformally decomposable in this face, i.e., is not elementary in CF C RC D k ≤η . It follows that the elementary flux vectors of CF C RC Bc ≤η are made up of all the elementary flux vectors of the CF C RC D k ≤η 's. We can sum up the results regarding EFVs as: This means that EFVs can be computed flux tope by flux tope and constraint-disjunct by constraintdisjunct. And the result holds also for EFPs (by considering vertices of F P ≤ η) and EFVs of CF P RC Bc .
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 CF C RC D k ≤η is not necessarily an EFM of CF C RC Bc ≤η .

Example 2.2. Consider the network comprising three irreversible reactions and one internal metabolite
A: R1 : → A, R2 : → A, R3 : A →, and the Boolean constraint Bc = ¬R1 ∨ R2, decomposed as Bc = D 1 ∨ D 2 , with D 1 = ¬R1 and D 2 = R1 ∧ R2. Take η = +. Then e 1 = (0, 1, 1) T , made up of R2, R3, which is the only EFM of CF C RC D 1 ≤η , is also the only EFM of CF C RC Bc ≤η . However any e 2 = (λ, 1, 1 + λ) T with λ > 0, made up of R1, R2, R3, is an EFM of CF C RC D 2 ≤η and is not an EFM of CF C RC Bc ≤η . Note that the way the constraint is decomposed matters. With the decomposition Bc = D 1 ∨ D 2 with D 1 = R2 and D 2 = ¬R1 ∨ ¬R2, the result for D 1 is identical to that for D 1 , but This means that, if it is natural to study each CF C RC D k ≤η (resp., CF 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: EFMs Proposition 2.11. Given an arbitrary regulatory constraint RC Bc with Bc = k D k , where the disjuncts D k are taken two by two inconsistent, and the associated constrained flux cone subset CF C RC Bc (resp., constrained flux polyhedron subset CF P RC Bc ), 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 CF C RC D k ≤η (resp., CF 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 Bc . This is not the case for EFMs. First, an EFM of CF C RC Bc ≤η is not necessarily an EFM of CF C RC Bc , but actually the biologically relevant minimality concept being signminimality and not support-minimality, we will consider EFMs for each flux tope CF C RC Bc ≤η . Second, an EFM of CF C RC D k ≤η is not necessarily an EFM of CF C RC Bc ≤η . EFMs of CF C RC Bc ≤η are actually obtained by collecting EFMs of CF 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 CF C RC D ≤η (resp., CF 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 ≤η ) ∩ Sol RC D those EFMs of F C ≤η that satisfy the constraint RC D and EFMs RC D = EFMs(CF 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 2.8) constraints, there is generally no longer identity between EFMs RC D and EFMs RC D [28]. Example 2.3. Consider the network of Example 2.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 EFMs RC D \EFMs RC D .
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 thatF ⊆ 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 ofF 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 = cone ⊕ ({e k }) (which is actually the same as cone({e k }) as we are in orthant O ≤η ) and thusF = { k∈K β k e k | β k > 0} = cone + ⊕ ({e k }) and the common minimal support of all vectors ofF is supp(F ) = k∈K supp(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 ofF 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). So, we consider a face F of F , of dimension at least two, such thatF ⊆ 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 supp(F ) = 1≤k≤N supp(e k ). Note that, for any k, supp(e k ) ∩ I = ∅ (and, as we have seen, I supp(e k )), because if for a certain k 0 we had supp(e k0 ) ∩ I = ∅, then e k0 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's note S 1 = supp(e 1 ) ∩ I and, for any k, 2 ≤ k ≤ N , S k = (supp(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 k0 = ∅, then the vectors of cone + ⊕ ({e k } k =k0 ) would verify the constraint RC D and have their supports included in supp(F ), thus equal to it as it is minimal for vectors in CF 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 supp(e k ) ⊇ S k by construction and supp(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 supp(e k ), where {e k } 1≤k≤N , N ≥ 2, are vectors in EFMs(F ) verifying supp(e k ) ⊇ S k and supp(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 supp(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 . So, 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 supp(EFMs RC D ). We have thus proved the following result. for all i ∈ I. Let EFMs RC D = EFMs(F C ≤η ) ∩ Sol RC D be the EFMs of F C ≤η that satisfy RC D and EFMs RC D = EFMs(CF C RC D ≤η ) be the EFMs of CF C RC D ≤η . We get EFMs RC D ⊆ EFMs RC D but, unlike the case of thermodynamic and kinetic (as in Proposition 2.8) constraints, there is generally no longer identity between EFMs RC D and EFMs RC D . EFMs RC D correspond to the extreme rays (faces of dimension one) of F • , i.e., the edges of F whose support contains I. EFMs RC D \EFMs RC D are the vectors of the F 's for all F faces of F of dimension at least two, such thatF ⊆ 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 EFMs RC D \EFMs RC D are obtained as 1≤k≤N supp(e k ), where {e k } 1≤k≤N , N ≥ 2, are vectors in EFMs(F ) verifying supp(e k ) ⊇ S k and supp(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 EFMs 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 Gv ≥ h, may not respect the structure of CF C RC Bc 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 {v i = 0}, in order to determine EFMs RC D \EFMs RC D . And in practice, ILC is generally only used to bound (below and/or above) fluxes, in which case each extremal ray R + 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) nondecomposable 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 2.3 can be decomposed as (0.7e 2 +0.3e 3 )+(0.3e 2 +0.7e 3 )). More straightforwardly, we can notice that EFMs RC D is equal to EFVs(CF 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 supp(e 2 ) ∪ supp(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 thatF ⊆ F • , not covered by Proposition 2.12, i.e., owning at least one facet F with interior included in F • , it so happens that all vectors ofF are actually support-wise decomposable, as each such vector can always be decomposed into a vector ofF 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: 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 nonstrictly-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: (84) 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 2.12, let swNSDFVs RC D = swNSDFVs(F C ≤η ) ∩ Sol RC D and swNSDFVs RC D = swNSDFVs(CF 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 2.8) 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 2.5. Continuing with the network of Example 2.3, e 1 + e 3 ∈ swNSDFVs RC D \EFMs RC D . More precisely, we have (by considering only a representative for each ray) So, we extend the results of Proposition 2.12 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 = cone ⊕ ({e m } m∈M ). We get the following structure for F • regarding EFMs and swNSDFVs, given in the form of an algorithm: 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 sign(CF C RC D ), this case corresponds to I = ∅). We consider the case R ⊂ M here below.
• Let's consider successively all faces F of F of dimension at least two, such thatF ⊆ 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 . ThusF = cone + ⊕ ({e k } k∈K ) and I ⊆ k∈K supp(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 supp(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, thenF ⊆ swNSDFVs RC D \EFMs RC D and k∈K supp(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 supportwise non-strictly-decomposable, independently of the choice of the nonnegative 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 inF 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 inF .
• 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 cone + ⊕ ({e l } l∈L ), is not empty inF (as there are at least two such facets) and is made up of the support-wise strictlydecomposable vectors ofF (by construction): cone + ⊕ ({e l } l∈L ) =F \swNSDFVs RC D . Two subcases must therfeore be distinguished.
• If L = K, i.e., the strict conical sum of the interiors of these facets is equal toF , thenF ⊆ F • \swNSDFVs RC D and k∈K supp(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 strictlydecomposable, independently of the choice of the nonnegative coefficients fixing the contribution of each e k in the distribution of the fluxes).
• If L ⊂ K, thenF is split into two nonempty subsets: is made up of the vectors of cone + ⊕ ({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 ofF are support-wise non-strictlydecomposable and part are support-wise strictly-decomposable, while having the same non-minimal support k∈K supp(e k ) (still equal to l∈L supp(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-strictlydecomposable support. Note that, in this subcase, the support-wise non-strictly-decomposable vectors ofF constitute the complementary, in the open cone cone + ⊕ ({e k } k∈K ), of the open subcone cone + ⊕ ({e l } l∈L ) which is thus a finite disjoint union of semi-open cones, each of which is conically generated (with positive or nonnegative 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 2.12).
Proposition 2.13. Let swNSDFVs RC D = swNSDFVs(CF C RC D ≤η ) be the support-wise non-strictlydecomposable vectors of CF C RC D ≤η . We obtain EFMs RC D ⊆ swNSDFVs RC D but, unlike the case of thermodynamic and kinetic (as in Proposition 2.8) constraints, there is in general no longer identity between EFMs RC D and swNSDFVs RC D . Consider all faces F of F of dimension at least two, such that {(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. And 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 +ye 3 ), with arbitrary k, 0 < k < 1, i.e., into two support-wise non-strictly-decomposable vectors respectively in {(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 xe 2 + ze 3 + (y − z)e 4 and belongs to cone + ⊕ ({e 2 , e 3 , e 4 }) if z < y and to cone + ⊕ ({e 2 , e 3 }) if z = y.
However take care, in the case of a general constraint Bc = k D k , that, if the support-wise strictlydecomposable vectors of each CF C RC D k ≤η are support-wise strictly-decomposable in CF C RC Bc ≤η , it is not true for support-wise non-scrictly-decomposable vectors and some of them for CF C RC D k ≤η may be decomposable in CF C RC Bc ≤η .
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.

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 Bc ), 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: . Now, from propositions 2.5 and 2.8, CF C C1 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: From proposition 2.10, each CF C i C2 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 2.11, we can also conclude that the EFVs of CF C C1∧C2 are exactly the EFVs of F C that satisfy both constraints C 1 and C 2 .
Proposition 2.14. 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 Bc 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 2.13 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.

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 CF C C or CF P C . More precisely, in identifying certain general features regarding constraints C (some of which are present in particular in biological constraints (41), (46),(49) 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.

Sign-invariant constraints
A constraint C x (v) (resp., ∃xC x (v)), is said to be sign-invariant if it depends only on sign(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): It follows from these definitions that, if the C x (v)'s are sign-invariant for all x, then ∃xC x (v) is signinvariant. Note that the property for a constraint to be support-invariant (i.e., to depend only on supp(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. On the other hand, the kinetic constraint KC E,M (v) (46) is not sign-invariant, and this is also the case of KC E (v) (77), 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) (79) 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 sign(v ) = sign(v), then, for all i ∈ supp(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  The structure of the solution space of a sign-invariant constraint follows directly from its definition. Actually, if v satisfies a given sign- union of open orthants. This result can be compared to Lemma 2.2. More precisely, if we consider a support-invariant constraint C, the same reasoning applies and gives Sol C = {v|C(v)} S supp(v) , where S ρ = {x ∈ R r | supp(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 Bc = i D ρ i , thus to the regulatory constraint RC Bc . 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 |supp(ρ)| sign vectors η with support ρ, given by η i = + or − if ρ i = 1 and 0 else, and we obtain S ρ = {η i |supp(η i )=ρ}O η i and Sol C = {η i |supp(η i )=supp(v)∧C(v)}O η i , i.e., that Sol C is a disjoint union of open orthants, which is precisely the statement of Lemma 2.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 Sol C can be reconstituted from Sol C ∩ O η . And, 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.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 3.2. 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 Sol C of vectors in R r satisfying the sign-invariant constraint C = C x (v) (resp., C = ∃xC x (v)) is a disjoint union of open orthants, which can be grouped together according to Lemma 2.2 to provide a disjoint union of semi-open orthants without any possible further merging between any two of them. Theorem 3.3. Let C = C x (v) (resp., C = ∃xC x (v)) be a sign-invariant constraint and CF C C (resp., CF P C ) be the associated constrained flux cone subset (resp., the associated constrained flux polyhedron subset). Then CF C C (resp., CF 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 sign(F C) (resp., sign(F P )). They can be grouped together according to Lemma 2.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): , 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 3.3 applies in particular to all thermodynamic constraints, to regulatory constraints and to those kinetic constraints described in Lemma 3.1. Especially, we directly obtain the results of Propositions 2.10 and 2.11 for regulatory constraints.
It is important to note that our only knowledge of the EFVs for CF C C (resp., the EFPs and EFVs for CF P C ) is really a long way from characterizing the solution space CF C C (resp., CF P C ): actually they are just the ExVs of each tope CF C C ≤η (resp., the ExPs and ExVs of each CF 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 CF C C (resp., CF P C ).
As in each closed r-orthant, i.e., in each flux tope, a sign-invariant constraint identifies with a regulatory constraint, Propositions 2.12 and 2.13 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 CF 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 CF C RC D ≤η .
Let us adopt notations similar to those used in Propositions 2.12 and 2.13, i.e., let EFMs C = EFMs(F C ≤η )∩ Sol C be the elementary flux modes of F C ≤η that satisfy C (i.e., the elementary flux vectors of CF C C ≤η ), EFMs C = EFMs(CF C C ≤η ) be the elementary flux modes of CF C C ≤η and swNSDFVs C = swNSDFVs (CF C C ≤η ) be the support-wise non-strictly-decomposable vectors of CF C C ≤η . We have thus EFMs C ⊆ EFMs C ⊆ swNSDFVs C . The proof of Proposition 2.12 adapts straightforwardly: EFMs C correspond to the open faces of dimension one (extreme rays) of F C ≤η in CF C C ≤η and EFMs C \EFMs C are made up of all the open facesF in CF 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 thatF is in CF C C ≤η . The result is based on the properties that all vectors of an open faceF in CF C C ≤η have the same support, and the common support of vectors ofF , where F is a face of F , is strictly included in that of vectors ofF . In the same way, the proof of Proposition 2.13 is easily adapted for what concerns swNSDFVs C \EFMs C . For this we consider successively all the open faces F in CF C C ≤η where F is a face of dimension at least two of F C ≤η which owns at least one proper face F such thatF is in CF C C ≤η . For each such givenF , 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 cone + ⊕ ({e k } k∈K ) =F and cone + ⊕ ({e l } l∈L ) = cone + ⊕ ({F j } j∈J ) and we obtainF ∩ cone + ⊕ ({e l } l∈L ) =F \swNSDFVs C , as by construction those vectors ofF which belong to cone + ⊕ ({F j } j∈J ) are precisely those inF which are support-wise strictly-decomposable in CF C C ≤η , the decomposition being achieved along vectors in theF j 's, whose support is strictly included in the common support of vectors inF . So three different exclusive cases can be distinguished forF : , 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 supp(e l ) ⊂ k∈K supp(e k ) (this is for example the case if |J| = 1). ThenF ⊆ swNSDFVs C \ EFMs C is entirely made up of support-wise non-strictly-decomposable vectors and k∈K supp(e k ) is a non-strictly-decomposable non-minimal support for vectors in CF C C ≤η .
• cone + ⊕ ({e l } l∈L ) =F , i.e., L = K. ThenF ⊆ CF C C ≤η \swNSDFVs C is entirely made up of support-wise strictly-decomposable vectors and k∈K supp(e k ) is a strictly-decomposable support for vectors in CF C C ≤η .
• cone + ⊕ ({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 l i = 0, which is still equivalent to l∈L supp(e l ) = k∈K supp(e k ). ThenF is split into two nonempty subsets: cone + ⊕ ({e l } l∈L ) =F \swNSDFVs C , made up of support-wise strictly-decomposable vectors, andF \cone + ⊕ ({e l } l∈L ) ⊆ swNSDFVs C \EFMs C , made up of support-wise non-strictly-decomposable vectors (note thatF \cone + ⊕ ({e l } l∈L ) is made up of those vectors ofF the conical decomposition of which on the e k 's requires at least one e k with k ∈ K\L). As all vectors ofF have the same non-minimal support k∈K supp(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-strictlydecomposable support. Note that the support-wise non-strictly-decomposable vectors ofF are obtained as the complementary, in this open cone, of the open sub-cone cone + ⊕ ({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 nonnegative 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 cone ⊕ ({e l } l∈L )).
We have finally proved the following result (while keeping the notations of Theorem 3.3).
Theorem 3.4. For C a sign-invariant constraint, F C ≤η a flux tope and CF C C ≤η the associated subset of the solution space (disjoint union of open faces of F C ≤η from Theorem 3.3), let EFMs C = EFMs(F C ≤η )∩ Sol C be the elementary flux modes of F C ≤η that satisfy C (equal to EFVs(CF C C ≤η )), EFMs C = EFMs(CF C C ≤η ) be the elementary flux modes of CF C C ≤η , and swNSDFVs C = swNSDFVs (CF C C ≤η ) be the support-wise non-strictly-decomposable vectors of CF C C ≤η , with EFMs C ⊆ EFMs C ⊆ swNSDFVs C . EFMs C correspond to the open faces of dimension one (extreme rays) of F C ≤η belonging to CF C C ≤η . Consider now successively all the open facesF in CF 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 thatF belongs to CF C C ≤η . Thus cone + ⊕ ({e k } k∈K ) =F and let D(F ) = cone + ⊕ ({e l } l∈L ). Consequently: if L = ∅, thenF ⊆ EFMs C \EFMs C ; if L = ∅ and l∈L supp(e l ) ⊂ k∈K supp(e k ), thenF ⊆ swNSDFVs C \ EFMs C and k∈K supp(e k ) is a non-strictly-decomposable non-minimal support for vectors in CF C C ≤η ; if L = K, thenF ⊆ CF C C ≤η \swNSDFVs C and k∈K supp(e k ) is a strictly-decomposable support for vectors in CF C C ≤η ; if L ⊂ K and l∈L supp(e l ) = k∈K supp(e k ), then D(F ) =F \swNSDFVs C , a non-empty open cone, constitutes the support-wise strictly-decomposable vectors inF andF \D(F ) ⊆ swNSDFVs C \EFMs C , a non-empty finite disjoint union of semi-open cones, constitutes the support-wise non-strictly-decomposable vectors inF , 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 EFMs C \EFMs C and swNSDFVs C \EFMs C by collecting these vectors for eachF .
Note that even the knowledge of swNSDFVs C is not enough to completely reconstruct the tope CF 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., ∃xC 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: Obviously a sign-invariant constraint is contracting-sign-invariant and, if the C x (v)'s are contractingsign-invariant for all x, then ∃xC x (v) is contracting-sign-invariant. M (v) is satisfied (for a certain E verifying c T E ≤ W and E ≤ E + ) and sign(v ) = sign(v) with |v i | ≤ |v i | for all i, we saw that, by taking We will call 0-star domain a subset SD of R r which has the property that the whole open segment joining 0 to any element of SD is included in SD: v ∈ SD, 0 < λ < 1 ⇒ λv ∈ SD. Any cone is a 0-star domain.
Theorem 3.6. If C = C x (v) (or C = ∃xC x (v)) is contracting-sign-invariant, then CF C C is a 0-star domain and there exists a neighborhood N = ]−δ, +δ[ r of 0 for a certain δ > 0 such that CF 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 (87) regarding EFVs and EFMs hold in N , i.e., for sufficiently small fluxes. The same holds for CF P C with the flux topes F P ≤η if 0 is an interior point of F P , and in this case results (88) regarding EFVs and EFPs hold in N .
Proof. The property of being a 0-star domain is a direct consequence of the definition: if a contractingsign-invariant constraint is satisfied by a vector v, it is satisfied by λv for all 0 < λ < 1. From the proof of Theorem 3.3, 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 = min i∈supp(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 Gv ≥ h, verify h < 0.
Theorem 3.6 tells us that the result of Theorem 3.3 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 3.5.

Sign-monotone constraints
We now consider a property of compatibility of a constraint with signs that is stronger than signinvariance. A constraint C x (v) (resp., ∃xC 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): It follows from these definitions that, if the C x (v)'s are sign-monotone for all x, then ∃xC 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, sign(v ) ≤ sign(v) is equivalent to supp(v ) ⊆ supp(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.  49) is not generally sign-monotone: consider for example Bc reduced to a positive literal. Actually, RC Bc (v) is sign-monotone if and only if Bc 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 signmonotone).
Lemma 3.7. 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 Bc is not sign-monotone, except if Bc 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 3.8. The set Sol C of vectors in R r satisfying the sign-monotone constraint C = C x (v) (resp., C = ∃xC x (v)) is a union of closed orthants.
Obviously we can keep only those v's for which sign(v) is maximal, in order to avoid any inclusion between the closed orthants. This result can be compared to Lemmas 2.1, 2.4 and first part of Proposition 2.8, 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 3.9. Consider indifferently a sign-monotone constraint C x (v) or ∃xC 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 CF C C (resp., the associated constrained flux polyhedron subset CF P C ). Then CF C C (resp., CF 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 sign(F C) (resp., sign(F P )) and we obtain: Proof. The (non-disjoint) decomposition of the solution space into faces of the flux topes F C ≤η (resp., F P ≤η ) follows from Lemma 3.8 or directly from the fact that all vectors of the conical hull cone(v 1 , . . . , v n ) = {β 1 v 1 +. . .+β n v n | β 1 , . . . , β n ≥ 0} and of the convex hull conv(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 supp(v 1 ) ∪ . . . ∪ supp(v n ), which is the support of any vector v in cone + (v 1 , . . . , v n ) or conv + (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 3.3 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 supp(v ) ⊂ supp(v). And we could choose λ ∈ R such that v = v − λv is nonzero, belongs to F C and verifies supp(v ) ⊂ supp(v) and sign(v ) ≤ sign(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. So, v is an elementary flux mode in F C and we get the result for EFMs.
Of course, in the decomposition of CF C C (x) or CF C C (resp., CF P C (x) or CF 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 3.9 applies in particular to all thermodynamic constraints, to those kinetic constraints and to those very few regulatory constraints described in Lemma 3.7. In particular, we directly obtain (except of course the reference to ts M ) Proposition 2.5 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 2.8, 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 CF C C (resp., and of EFPs of CF P C ) is not good enough to reconstruct the structure of the solution space CF C C (resp., CF 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: cone ⊕ (E i ) ⊆ CF C C (resp., and conv ⊕ (E i ) ⊆ CF 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., ∃xC 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: Obviously a sign-monotone constraint is contracting-sign-monotone and, if the C x (v)'s are contractingsign-monotone for all x, then ∃xC x (v) is contracting-sign-monotone. Also, any contracting-sign-monotone constraint is contracting-sign-invariant.
Example 3.4. The argument given in Example 3.2 to establish that the kinetic constraints KC b M (v) (80) and KC b (v) (76) 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.
is contracting-sign-monotone, then CF C C is a 0-star domain and there exists a neighborhood N = ]−δ, +δ[ r of 0 for a certain δ > 0 such that CF 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 (93) regarding EFVs and EFMs hold in N , i.e., for sufficiently small fluxes. The same holds for CF P C with flux topes F P ≤η if 0 is an interior point of F P , and in this case results (94) regarding EFVs and EFPs hold in N .
Proof. The proof becomes identical to that of Theorem 3.6 by using the proof of Theorem 3.9.
Theorem 3.11 tells us that the result of Theorem 3.9 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 3.10, providing the result for the structure of CF C KC b M and CF C KC b quoted in the last part of Proposition 2.8.
Finally, proposition 2.14 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 signinvariant the conjunction of the two constraints is sign-invariant and Theorem 3.3 applies.
Theorem 3.12. 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 3.4 applies also to determine elementary flux modes and support-wise non-strictly-decomposable vectors.

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's briefly remember the principle of the DD method [29]. 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. And, 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 (6), what has been shown as the most efficient initialization is to start from a basis of the nullspace of S ( [45] 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 3.9, the EFVs (identical to EFMs) of the constrained flux cone subset CF 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 CF P C ). So, 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 3.13. Given a sign-monotone constraint C, the EFVs (or EFMs) of the constrained flux cone subset CF 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 [11,13] 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 2.8), Proposition 2.6 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 CF 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 CF 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, 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 CF C C and the support-wise non-strictly-decomposable vectors in CF C C and, as only the last two are biologically relevant, computing only the first ones would be of limited interest.

Conclusion
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 or elementary flux vectors, which 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, 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 and, as there are considerably fewer of them, improving the scalability of this computation, we took in consideration in this paper one 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 faces of the flux topes of F C) and that the EFMs, which 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 computable simply by adding one homogeneous linear inequality to the initial (stoichiometric and reactions irreversibility) ones. The situation is more complex for constraints which are only sign-invariant (i.e., once satisfied by a certain vector are satisfied by any vector with the same sign). We demonstrated that in fact support-invariant constraints coincide with regulatory (Boolean) constraints, so, in this sense, this is not a generalization. For sign-invariant constraints (i.e., support-invariant in each closed r-orthant separately), we demonstrated that the solution space is a finite disjoint union of open (i.e., without their facets) convex polyhedral cones 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 to EFMs and that EFMs identify to support-wise non-strictly-decomposable vectors: there are in general 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. Firstly 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. And this generalizes to a metabolic network modeled by a flux polyhedron F P by using EFPs. This gives an indirect method for computing MCSs [15,18] 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 [4,44], 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). Secondly, 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. And 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 2.7 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 open polyhedral cones (or open polyhedra) is more complex and other 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 [28,26] will be investigated.