### Computational Procedure

Our algorithm needs three input objects: the stoichiometric matrix *S* ∈ ℝ ^{
m×n
} of the network *is*
, the set of irreversible reactions *Irr* ⊆ {1, . . ., *n*}, and the set of reactions ∑ ⊆ {1, . . ., *n*} in the subnetwork of interest, while as an output it will return the complete set of ProCEMs. The computation of ProCEMs is achieved in three main consecutive steps.

**Step 1 - Preprocessing:** The aim of this step is to remove inconsistencies from the metabolic network and to transform it into a form suitable for the projection in Step 2. First, based on ∑ we sort the columns of

*S* in the form:

$\stackrel{\u0304}{S}=\left(\u0100\phantom{\rule{0.3em}{0ex}}\stackrel{\u0304}{B}\right)$

(3)

where the reaction corresponding to the

*i*-th column belongs to ∑ iff the

*i*-th column is in

*Ā*. Next, the blocked reactions [

37] are removed. Finally, each of the reversible reactions is split into two irreversible "forward" and "backward" reactions. The final stoichiometric matrix will be in the form:

${S}^{\prime}=\left(A\phantom{\rule{0.3em}{0ex}}B\right)$

(4)

where the columns of

*A* represent the "interesting" reactions after splitting reversible reactions and removing the blocked reactions. In the following, we assume that

*A* (resp.

*B*) has

*p* (resp.

*q*) columns. Given

*S*', the steady-state flux cone in canonical form will look as follows

$C=\left\{\left(x,y\right)\in {\mathbb{R}}^{p+q}|G\cdot x+H\cdot y\le 0\right\},$

(5)

where matrix

*G* (resp.

*H*) represent the columns to be kept (resp. eliminated):

$G=\left(\begin{array}{c}\hfill -A\hfill \\ \hfill A\hfill \\ \hfill -{I}_{p}\hfill \\ \hfill {0}_{q,p}\hfill \end{array}\right),\phantom{\rule{0.3em}{0ex}}H=\left(\begin{array}{c}\hfill -B\hfill \\ \hfill B\hfill \\ \hfill {0}_{p,q}\hfill \\ \hfill -{I}_{q}\hfill \end{array}\right)$

(6)

Here *I*_{
p
} denotes the *p* × *p* identity matrix, and 0 _{
p,q
} the *p* × *q* zero matrix.

**Step 2 - Cone Projection:** In this step, the flux cone is projected, eliminating the reactions corresponding to columns in

*H*. Several methods have been proposed in the literature for the projection of polyhedra [

48]. For our purpose we chose the

*block elimination method*[

49]. This method allows us to find an inequality description of the projected cone by enumerating the extreme rays of an intermediary cone called the

*projection cone*. In our case, the projection cone is defined as

$W=\left\{w\in {\mathbb{R}}^{2m+p+q}|{H}^{T}\cdot w=0,w\ge 0\right\},$

(7)

where *H*^{
T
} denotes the transpose of *H*.

We enumerate the extreme rays {

*r*^{1},

*r*^{2}, . . .,

*r*^{
k
} } of

*W* using the double description method [

50]. The projected cone is given by

${\mathcal{P}}_{X}\left(C\right)=\left\{x\in {\mathbb{R}}^{p}|R\cdot G\cdot x\le 0\right\},$

(8)

where

$R={\left({r}^{1}...{r}^{k}\right)}^{T}.$

(9)

This representation of the projected cone contains as many inequalities as there are extreme rays in *W*, thus a large number of them might be redundant [48]. These redundant inequalities are removed next (see below).

**Step 3 - Finding ProCEMs:** In the final step, the extreme rays of the projected cone, i.e., the ProCEMs, are enumerated. Similarly as in Step 2, the double description method is employed to enumerate the extreme rays of ${\mathcal{P}}_{X}\left(C\right)$.

With the block elimination algorithm, it is also possible to perform the projection in an iterative manner. This means that rather than eliminating all the "uninteresting" reactions in one step, we can partition these in *t* subsets and then iteratively execute Step 2, eliminating every subset of reactions one by one. By proceeding in this fashion, the intermediary projection cones, *W*^{1}, *W*^{2}, . . ., *W*^{
t
} get typically smaller, thus enumerating their extreme rays requires less memory. On the other side, the more sets we partition into, the slower the projection algorithm usually gets.