A Graph-Based Approach to Analyze Flux-Balanced Pathways in Metabolic Networks

TL;DR

This paper introduces a graph-based method for identifying flux-balanced pathways in metabolic networks, simplifying the analysis of elementary flux modes by leveraging graph models for reaction categorization.

Contribution

The paper proposes a novel graph data model approach to efficiently find a subset of EFMs, improving pathway analysis in metabolic networks.

Findings

Enables categorization of pathways based on metabolites and reactions

Reduces computational complexity in EFM analysis

Provides insights into metabolic network structure

Abstract

An Elementary Flux Mode (EFM) is a pathway with minimum set of reactions that are functional in steady-state constrained space. Due to the high computational complexity of calculating EFMs, different approaches have been proposed to find these flux-balanced pathways. In this paper, an approach to find a subset of EFMs is proposed based on a graph data model. The given metabolic network is mapped to the graph model and decisions for reaction inclusion can be made based on metabolites and their associated reactions. This notion makes the approach more convenient to categorize the output pathways. Implications of the proposed method on metabolic networks are discussed.

Tables

Table 1. Table 2: Number of EFMs starting from Glucose in some networks.

Network	Size ($m\ast n$)	#EFMs
Glycolysis and Pentose Phosphate pathway in $T.$ $brucei$	26$\ast$35	4
Human red blood cell	20$\ast$50	20
$E.$ $coli$ core (Escherichia coli str. K-12 substr. MG1655)	53$\ast$94	47

Equations

$$r_j:m_1+m_2+...+m_i+...\rightleftharpoons m_1^{\prime }+m_2^{\prime }+...+m_i^{\prime }+...$$

$$\begin{array}[]{l}O_{f_{ijp}}={\rm}\frac{{I_{f_{ijp}}I_{c_{ij}}}}{{O_{c_{ij}}}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,Forward\\ \\ I_{f_{ijp}}={\rm}\frac{{O_{f_{ijp}}O_{c_{ij}}}}{{I_{c_{ij}}}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,Backward\\ \end{array}$$

$$\begin{array}[]{l}I_{f(n)_{ijp}}=\frac{{I_{f(o)_{ijp}}I_{c_{ij}}-I_{f_{iep}}I_{c_{ie}}}}{{I_{c_{ij}}}},\\ \\ I_{f(n)_{ijp}}=\frac{{I_{f(o)_{ijp}}I_{c_{ij}}+O_{f_{iep}}I_{c_{ie}}}}{{I_{c_{ij}}}},\\ \\ O_{f(n)_{ijp}}=\frac{{O_{f(o)_{ijp}}O_{c_{ij}}+O_{f_{iep}}I_{c_{ie}}}}{{O_{c_{ij}}}},\\ \\ O_{f(n)_{ijp}}=\frac{{O_{f(o)_{ijp}}O_{c_{ij}}-I_{f_{iep}}I_{c_{ie}}}}{{O_{c_{ij}}}}\\ \end{array}$$

Full text

Keywords: Elementary Flux Mode (EFM); Graph Data Model; Metabolic Network

1 Introduction

Metabolic network models are among the well-studied models in biotechnology. The reconstruction of these networks is possible by collecting the gene-protein-reaction information from related genomic data and literature [1]. It is important to explore biologically relevant pathways in metabolic networks. Forcing constraints to a reconstructed biochemical network results in the definition of achievable cellular functions [2]. Mathematical representation of constraints are as flux-balance constraints (e.g., conservation of flux) which means the network should be at the steady-state condition, and flux bounds which limit numerical ranges of network parameters and coefficients such as the minimum and maximum range of fluxes for each reaction. Flux Balance Analysis (FBA), a method to predict the optimal growth rate of a a certain species when grown on particular set of metabolites [3], and Elementary Flux Mode (EFM) analysis, an approach to decompose a network to minimal functional pathways [4], are among the approaches used in constraint-based analysis of metabolic networks.

EFMs [5, 6] have been used in several biological applications such as bioengineering [7], phenotypic characterization [8], drug target prediction [9] and strain design [10]. The incorporation of kinetic analysis into EFMs enables a more complete description of cellular functions for which kinetics play a dominant role [11].

Several methods were introduced for finding EFMs based on double-description method [5]. Double-description is a technique to enumerate all extreme rays of a polyhedral cone. An improved approach to the primary method was introduced in which the null-space of the stoichiometric matrix is used instead of the matrix itself to generate EFM candidates [12, 13, 14]. Various effective computational approaches have been proposed to speed-up previous methods for computing EFMs [15, 16]. Some of these approaches led to the development of computational tools such as Metatool [17] and EFMtool [16]. Besides, methods based on linear programming have been proposed that explore a set of EFMs with specific properties, such as $K$-shortest EFMs [18], or EFMs with a given set of target reactions [19]. A new set of methods based on graph-theory, [20, 21] tries to overcome the scalability problem of the double-description-based techniques.

According to the complexity of computing EFMs as discussed in [22, 23], an additional challenge of extracting biological properties from large set of EFMs [11] has been raised. As stated in [23], the complexity of enumerating all EFMs still remains unknown and the computational complexity of enumerating EFMs containing a specific reaction is hard. Even for networks with the same number of reactions and metabolites, knowing the number of EFMs in one network cannot necessarily help in finding the number of EFMs in the other one, as they have different connectivities and topologies. This fact emphasizes the inherent structural information reflected by EFMs [24] and is a motivation to extract a set of biologically meaningful EFMs according to a biological aspect such as motifs [25] and thermodynamics [26].

The main focus of this paper is to introduce a data model based on the AND/OR graph and to propose an approach to find flux-balanced pathways according to pathway topology and reaction stoichiometries. It is shown that the computed flux-balanced pathways are a biologically relevant subset of EFMs which include external input and output metabolites. In other words, the subset comprises all EFMs connecting input and output external metabolites. Besides, based on the introduced graph data model, an upper-bound for the complexity of exploring EFMs containing external metabolites is calculated. Using the topology of the network gives us the opportunity to make decisions according to the metabolite/reaction aspects to preserve or eliminate a particular reaction or metabolite in a certain pathway. Introducing a model to consider both topology and stoichiometry of a metabolic network for network analysis may lead to better biological decisions and output categorization. The graph structure makes it easier to set-up rules for pathways and to explore the intended solution space via rules. Finally, our method can potentially be implemented on hardware (through a system design approach) more conveniently.

The rest of the text is organized as follows. Some required concepts and preliminaries are provided in Section 2. Before explaining the proposed approach to find the elementary flux modes in a given metabolic network (Section 3.2), we first describe the modified AND/OR graph model and its properties in Section 3.1. Finally, Section 4 is devoted to results and discussions and Section 5 concludes the paper.

2 Preliminaries

In this section, some basic concepts along with the formal definition of elementary flux modes are provided.

2.1 Metabolic Networks

Metabolic networks model the metabolism of living cells in terms of a set of biochemical reactions. Definitions in this section are derived from [27]. The biochemical reactions can be irreversible, i.e., the reaction can be active only in one direction, or reversible, i.e., the reaction can be active in both directions. The contributing metabolites in a reaction can be either substrates or products. Substrates are consumed and products are produced during the operation of a reaction. The topology of a metabolic network is characterized by its $m\times n$ stoichiometric matrix, $S$, where $m$ and $n$ correspond to the number of metabolites and reactions, respectively. In this paper, external metabolites are not included in $S$. The value $S_{ij}$ represents the stoichiometric coefficient of the metabolite $i$ in the reaction $j$. $S_{ij}$ is positive/negative if the metabolite $i$ is produced/consumed. If this coefficient is zero it means that the metabolite $i$ does not contribute to the reaction $j$. The network is considered in the steady-state if for each internal metabolite, the rates of consumption and production are equal. The reactions connected to the external metabolites are called Boundary reactions.

2.2 Elementary Flux Mode Definition

A flux vector $\textbf{v}\neq\textbf{0}$ and $\textbf{v}\in{R}^{n}$ is considered an EFM if it meets the following conditions:

• •

$v_{i}\geq 0$ for all $i\in$ irreversible reactions (thermodynamic constraint).

• •

$\textbf{S}.\textbf{v}=\textbf{0}$ (steady-state condition).

• •

There is no $\textbf{v}^{\prime}\in{R}^{n}$ with $\textbf{supp}(\textbf{v}^{\prime})\subset\textbf{supp}(\textbf{v})$, where support of a mode is defined as $\textbf{supp}(\textbf{v})=\{i|v_{i}\neq 0\}$, (minimality limitation).

3 Methods

In this section, the proposed graph-based elementary flux mode analysis is discussed in detail. To do so, the data structure used in our approach is defined, and then, the proposed approach is introduced accordingly.

3.1 Data Model

The proposed data model is based on the conventional AND/OR graph in computer science. However, we provide a different definition with additional features to make the model appropriate for our proposed algorithm. In our model, the coefficients of the metabolites in reactions are embedded in the graph structure as attributes of each node (explained in Definition 3, see below). A $pathway$ is defined as a set of reactions with their associated fluxes. The flux of each reaction in a certain pathway is also embedded in the graph.

Definition 1. The Modified Graph for representing a metabolic network, denoted as $\mathbf{MG}$, is defined as a set of Nodes, $\mathcal{N}_{i}$, i.e., $\mathbf{MG}$ = $\{\mathcal{N}_{i}|0\leq i\leq{M-1}\}$, where $M$ is the number of internal metabolites in the metabolic network.

Definition 2. Each node $\mathcal{N}_{i}$ in $\mathbf{MG}$ is a 3-tuple $\mathcal{N}$ = ($i$, $I$, $O$) where

• •

$i$ is the tag of a metabolite,

• •

$I$ is an array of input reactions that produce the metabolite and

• •

$O$ is an array of output reactions that consume the metabolite.

Definition 3. Each $I$/$O$ in Definition 2 contains the following data:

• •

The reaction $j$, $0\leq j\leq r-1$, where $r$ is the number of reactions consuming/producing the metabolite $i$,

• •

$I_{M_{j}}$/$O_{M_{j}}$, an array of the metabolites consumed/produced by reaction $j$. In other words, $I_{M_{j}}$/$O_{M_{j}}$=$\{m_{kj}|0\leq k\leq m-1\}$, where $m$ is the number of consumed/produced metabolites by the reaction $j$,

• •

I_{{\mathord{\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\frown}\over{M}}}_{j}}/O_{{\mathord{\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\frown}\over{M}}}_{j}}, an array of the metabolites produced/consumed by reaction $j$. In other words, I_{{\mathord{\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\frown}\over{M}}}_{j}}/O_{{\mathord{\buildrel\lower 3.0pt\hbox{\scriptscriptstyle\frown}\over{M}}}_{j}}=$\{m_{kj}|0\leq k\leq m-1\}$, where $m$ is the number of produced/consumed metabolites by the reaction $j$ excluding the metabolite $i$ itself,

• •

The direction of the reaction for reversible reactions,

• •

$I_{c_{ij}}$/$O_{c_{ij}}$, the coefficient of the reaction $j$ for the produced/consumed metabolite $i$ in the stoichiometric matrix $S$ and

• •

$I_{f_{ijp}}$/$O_{f_{ijp}}$, the flux of the input/output reaction $j$ for a certain pathway $p$.

A general configuration of a node $\mathcal{N}_{i}$ in $\mathbf{MG}$ is shown in Figure 1. The incoming arcs (i.e., $I$ in Definition 2) in each $\mathcal{N}_{i}$ produce the metabolite $i$ and the outgoing arcs (i.e., $O$ in Definition 2) consume it. Therefore, consumed metabolites, $m_{i}$, and produced metabolites, $m^{\prime}_{i}$, contributing to reaction $j$ as shown in Equation 1, are as arcs between $\mathcal{N}_{i}$ nodes each associated to one metabolite. $I$ arcs and $O$ arcs in $\mathcal{N}_{i}$ nodes are related to each other by reaction tags in Definition 3.

[TABLE]

3.2 GB-EFM: Graph-Based EFM Analysis Algorithm

The proposed algorithm for finding EFMs of a metabolic network is based on the following facts:

Starting from boundary metabolites, each incoming flux in a pathway produce some metabolite(s), and each produced metabolite should be consumed and the flux of the contributed reactions can be obtained according to the stoichiometric coefficients in order to make sure that the constructed pathway is in steady-state. 2. 2.

To guarantee the elementarity of the produced pathways, two rules are followed: (1) Only one output from the node $\mathcal{N}_{i}$ should be considered when constructing a pathway and (2) multi-path condition should be checked in a node which has more than one input and one output. The multi-path condition occurs in a node when independent pathways share part of a pathway while traversing the graph. 3. 3.

When one metabolite is consumed by a reaction, the presence of its AND-related metabolites are required. For example, when metabolites $A$ and $B$ are AND-related, if $A$ is consumed during the reaction $A$+B$$\rightarrow$$C+$D$, then $B$ should also be present which means $A$ and $B$ should be consumed simultaneously. Therefore, the algorithm is designed to have a forward/backward flow for the produced/consumed metabolites. In this example, the metabolites $A$ and $B$ are consumed while $C$ and $D$ are produced.

The proposed graph-based EFM analysis algorithm is composed of the following five steps:

Step 0. Construction of $\mathbf{MG}$. The graph $\mathbf{MG}$ is constructed based on the stoichiometric matrix as stated in Algorithm 1.

Step 1. Finding independent pathways with dependent nodes. In this step, starting from the metabolites connected to the boundary reactions, all possible pathways are constructed by traversing the internal metabolites. Proposition 1 shows that the pathways are constructed with minimum number of possible reactions using Definition 4 and Definition 5. The semi-minimality occurs since different pathways are constructed from different outputs and each node is visited multiple times only if required.

Definition 4. Forward/Backward Flow. The forward flow is applied for the output metabolites of the selected reaction and the backward flow is applied for the and-related metabolites of that reaction.

Definition 5. Primary and Secondary Reactions. On each pathway, each node has a primary input and a primary output which directs the flow of the pathway and tags the metabolite as visited. When an edge (i.e., a reaction) enters an already visited node in that pathway in a forward/backward flow as an input/output, that reaction is tagged as a secondary input/output of that node.

Proposition 1**.**

Starting from metabolites connected to the boundary reactions, $r$ pathways are constructed from each node where $r$ is the number of output/input reactions of that node in the forward/backward flow. The pathways are semi-minimal based on the following statements:

• •

each output/input reaction in forward/backward flow is traversed once and labeled as primary, and

• •

when an edge enters an already visited node, the pathway is considered as closed in the node and the edge is labeled as secondary.

In the forward flow of a node, for each output reaction, an independent pathway is constructed. Each pathway is traversed independently of the other pathways. In the backward flow of a node, for each input reaction an independent pathway is constructed. In Figure 2 the forward/backward flow of each node $\mathcal{N}_{i}$ is shown in an illustrative example. The information of each pathway is saved independently.

Step 2. Eliminate multi-way pathways. By using Proposition 2 and eliminating multi-way pathways (defined in Definition 7) which include nodes with multi-path condition (defined in Definition 6), the set of semi-minimal paths is reduced to the set of minimal pathways. An example of this is shown in Figure 3.

Definition 6. Multi-Path Condition. This is a condition that occurs in a node with more than one input and more than one output and two independent pathways from different origins share this node.

Definition 7. Multi-Way Pathways. Semi-minimal pathways that contain reactions from two independent pathways are referred to as multi-way pathways.

Proposition 2**.**

Since different pathways are constructed from different outputs of a node and there are both forward and backward flows, depending on the order of traversing the nodes, there exist multi-way pathways containing reactions from two independent pathways. The pathway is minimal if it is not a multi-way pathway. In other words, it has no node with multi-path condition applied to it.

To solve the multi-path condition, for each produced pathway in Step 1, we traverse the subgraph of that pathway and tag each reaction with its associated metabolites. Since the structure of the pathway is known, one can label the outgoing/ingoing reactions in forward/backward flow of a pathway.

Then, in a pathway, for each node with more than one input and more than one output reaction, if one can find two routes crossing the node and the origins of the routes are different, the pathway should be discarded. This is checked according to the labels, i.e., $[node.pathNum]$, which has been assigned to each reaction. $pathNum$ is a variable which counts the number of pathways generated by a node. In other words, this pathway can be functional by either of the found routes in the node.

Step 3. Adjusting reaction fluxes. In this step, every node $i$ in a pathway $p$ is examined and the fluxes of the primary input and the primary output of that node are calculated using the following equations:

[TABLE]

In Equation 2, $j$ refers to the primary input/output reaction. Proposition 3 explains how all reactions in a pathway $p$ get a flux.

Proposition 3**.**

Starting from a boundary reaction with a certain flux, all nodes are being traversed and all reactions in a pathway get a flux using Equation 2, as each node has a primary input and a primary output to be used in this equation. Equation 2 states that the amount of the incoming and outgoing fluxes in $\mathcal{N}_{i}$ should be equal.

Step 4. Balancing flux. After Step 3, the consumption/production balance of the nodes with secondary reactions may be disturbed. To fix this, we go through these nodes and find a way forward/backward to output/input reactions to recalculate the extra production/consumption of the unbalanced nodes. Using Proposition 4, the pathways with no balance rate are discarded. If $\sum{I_{c_{ij}}I_{f_{ijp}}}-\sum{O_{c_{ij}}O_{f_{ijp}}}=0$, the balance condition is accepted for this node. Otherwise, on this pathway, the extra flux is added to the next reactions on the pathway, towards the boundary reactions using Equation 3. In this case, the difference between the input rates and the output rates is equal to $I_{c_{ie}}I_{f_{iep}}$ or $O_{c_{ie}}O_{f_{iep}}$ where $e$ implies the reaction which causes the extra production/consumption of the $\mathcal{N}_{i}$. The old/new value of the flux is shown by $o$/$n$ labels in Equation 3 and $j$ refers to the primary input/output reaction.

[TABLE]

Proposition 4**.**

Flux dependencies can prevent the rate of the production/consumption of an internal metabolite in a given topology of a pathway from being zero. Equation 3 is used to calculate an update for the consumption/production rate of the node. In the case of

• •

$I_{f(n)_{ijp}}=0$*             or             $O_{f(n)_{ijp}}=0$ or*

• •

a repetitive loop over a node (i.e., getting back to a node from the same reaction multiple times), while trying to find a way out in a subgraph of a given pathway,

the pathway is discarded.

An illustration of Steps 3 and 4 on a sample graph is shown in Figure 4.

The proposed GB-EFM analysis algorithm is presented in Algorithm 2 and is summarized in Theorem 1.

Theorem 1**.**

Starting from boundary nodes in a metabolic network, GB-EFM analysis algorithm produces all minimal flux-balanced pathways connecting input and output metabolites with the following properties:

• •

${R_{k}}\not\subset{R_{l}}$*,    for all    $\mathcal{P}_{k}$ and $\mathcal{P}_{l}$,       $k,l\in[0,P-1]$      $k\neq l$. *

• •

$\sum\limits_{j}{I_{f_{ijk}}I_{c_{ij}}}=\sum\limits_{j}{O_{f_{ijk}}O_{c_{ij}}}$,    for all nodes     $\mathcal{N}_{i}$,      $0\leq i\leq M-1$.

where $\mathcal{P}_{k}$ refers to a pathway $0\leq k\leq P-1$ and $R_{k}$ is the set of contributing reactions in the pathway. The above properties indicate that the produced pathways are elementary flux modes.

Proof.

Step 0 provides the required interface from the stoichiometric matrix. By using Propositions 1 and 2, semi-minimal and then minimal pathways are produced, respectively. According to Propositions 3 and 4, steady-state condition is checked for the given pathways and fluxes are calculated. Non-qualified pathways are discarded in Steps 2 and 4 based on Propositions 2 and 4. The remaining pathways have minimum functional reactions with balanced fluxes as elementary flux modes. Therefore, the produced pathways are a subset of EFMs connecting input and output external metabolites. ∎

A complete framework of Algorithm 2 is illustrated in Figure 5.

4 Results and Discussion

4.1 Computational Complexity

To calculate the computational complexity of Algorithm 2, the complexity of each step is analyzed below:

• •

In Step 1, the maximum number of created pathways is calculated as $P_{\rm{max}}=M_{B}\times(r_{\rm{max}})^{M_{\backslash B}}$, in which $M_{B}$ is the number of boundary metabolites, $M_{\backslash B}$ is the number of all other internal metabolites except the boundary ones and $r_{\rm{max}}$ is the maximum number of input/output reactions for all nodes.

• •

In Step 2, for each pathway $\mathcal{P}_{i}$, there are two traversals with the complexity of (1) $M_{\mathcal{P}_{i}}$, the number of all nodes in a pathway, and (2) $M_{t\mathcal{P}_{i}}$, the number of all tagged nodes with more than one input and one output in that pathway. Therefore, the order of maximum traversals is $P\times[M_{x\mathcal{P}_{i}}\times M_{\mathcal{P}_{i}}]$ or $P\times M^{2}$ where $M_{\mathcal{P}_{i}}$ and $M_{t\mathcal{P}_{i}}$ $\leq$ $M$.

• •

In Step 3, all nodes in $\mathcal{P}_{i}$ should be traversed which requires ${P^{\prime}}\times M_{\mathcal{P}_{i}}$ iterations. ${P^{\prime}}\leq P$ because some pathways may be discarded in Step 2. Therefore, the worst-case computational complexity for this step would be $P\times M$.

• •

In Step 4, for each $\mathcal{P}_{i}$ and for all nodes with secondary reactions, (i.e., $M_{v\mathcal{P}_{i}}$ nodes), at least half of the nodes should be traversed to lead the extra production/consumtion of metabolites to input/output. Therefore, ${P^{\prime}}\times[M_{v\mathcal{P}_{i}}\times M_{\mathcal{P}_{i}/2}]$ is the computational complexity of this step. Since $P^{\prime}\leq P$ and $M_{v\mathcal{P}_{i}},M_{\mathcal{P}_{i}/2}\leq M$, $P\times[M^{2}/2]$ is the computational complexity in the worst-case.

Considering all items above, $P\times O(M^{2})$, or as a result ${r_{\rm{max}}}^{M}\times O(M^{2})$ is the computational complexity of the algorithm if we assume $M_{B}\ll M$. The computational complexity of the approach can be considered as an upper-bound for the number of EFMs including external reactions.

4.2 EFM Topology Description

Three different topologies for an EFM can be observed in metabolic networks:

• •

Internal pathways with no boundary reaction(s),

• •

pathways with only input or only output reaction(s) consisting of an internal loop,

• •

pathways from input reaction(s) to output reaction(s).

These topologies are illustrated in Figure 6. The main focus of Algorithm 2 is the third topology. However, in order to expand this focus to all EFMs, the algorithm can be easily modified to keep pathways without the contribution of output reactions (in Case 2) rather than discarding them, or starting from internal reactions (in Case 1).

However, from the biological point of view, internal pathways with no boundary reactions, (i.e., loops) are called “futile cycles” and are not biotechnologically relevant [28], and the pathways with boundary inputs and no output implicate inconsistency in the network [29]. Therefore, a meaningful subset of EFMs are targeted here to keep more relevant EFMs in the output.

4.3 Application to Metabolic Network Models

To demonstrate the functionality of the model, the GB-EFM method has been implemented in C++ and tested on an Intel Core-i5 CPU with 4 GB RAM. The source code and results are provided on https://github.com/marabzadeh/GB-EFM. A small metabolic network comprising tricarboxylic acid cycle, glyoxylate shunt and adjacent reactions of amino acid synthesis in E. coli [30] is considered to describe the steps of the algorithm and proof the concept. The network has one boundary input reaction and four boundary output reactions among all 25 reactions and 18 internal metabolites as shown in Table 1. Reactions $R_{7}$, $R_{12}$, $R_{13}$ and $R_{24}$ are boundary output reactions. The elementary flux modes were obtained using EFMtool [16] as well as the proposed method. In the case of considering all reactions as irreversible, EFMtool produces 12 EFMs in which 3 are cycles, 2 are inconsistent routes with only an input reaction and no output (mentioned in [29] as elementary leakage modes) and the other 7 are consistent, non-cyclic pathways. Considering the reversible reactions in the table as reversible, 16 EFMs are produced, in which 3 and 2 are cycles and inconsistent routes belonging to the first and the second category as defined in Section 4.2, respectively, and 11 are pathways from input to an output as produced by GB-EFM. All information are provided in Supplementary file 1. The information of inconsistent routes with further description are provided in Supplementary file 2.

The number of pathways produced by the proposed algorithm after Step 1 is 40. After Step 2, this is reduced to 30 and after applying Steps 3-4, 5 pathways (final EFMs) remain. These pathways are the same as the eleven “acyclic” EFMs obtained by EFMtool. The resulting EFMs are illustrated in Figure 7.

Besides, a set of moderate-sized metabolic networks are considered as further test cases. The first network is the pentose phosphate pathway in $trypanosoma$ $brucei$ [32], the second one is the simple network of human red blood cell metabolism [33] and the third one is an $E.$ $coli$ core model obtained from [34]. GB-EFM found all EFMs from Glucose as reported in Table 2. Computational time for all the networks was less than one second. The currency metabolites such as ADP, ATP, AMP, CO2, H2O, O2, H2, NH4, Pi, NAD, NADH, NADP and NADPH were removed from the input stoichiometry matrix of all samples. As stated in the literature [35, 18], some simplifications to the system's model can be performed in order to reduce the complexity of the problem, such as setting currency metabolites like cofactors and energy metabolites to external. According to  [35], for energy currency metabolites like ATP, NADH and FADH, since their concentration is assumed to be constant they are not required to be balanced by an EFM. PG is an input metabolite and Asp$\_$ex, Glu$\_$ex, Ala$\_$ex and Suc$\_$ex are output metabolites.

4.4 Discussion

According to the graph data model, the order of traversing the nodes and the relations between reactions and metabolites can be maintained for each pathway. Besides, decisions to select or eliminate a particular reaction or metabolite can be made during the algorithm. Since in this approach the traversal is on the graph, many pathways are produced which may be discarded afterwards.

The proposed data model gives the opportunity of different levels of parallelism to trace pathways. Therefore, the complexity of implementation can be reduced by optimizing the method to only focus on pathways that are desirable in EFM analysis and by exploiting advantages of the possible parallelism in the method. The graph structure makes a good track of unbalanced metabolites. Our approach tries to decide if a certain “pathway” is EFM or not, merely based on topology rules (instead of recognizing the elementarity through rank test or comparing new candidates with produced ones).

Our strategy brings up the opportunity of object-oriented programming by considering a metabolite as an object. The method at least in its present form, might not be applicable to large-scale networks. However, the graph structure makes it easier to set-up rules for pathways and explore the intended solution space via rules. The proposed method can potentially be implemented through a system design approach on hardware more conveniently.

A hypergraph is a simplified type of the modified AND/OR graph. The application of hypergraphs in biological networks has been reviewed in  [36]. In the introduced model, each arc has two different input and output coefficients. Besides, for each arc, a dynamic label for each pathway is defined, that is, the flux of the reaction represented by that arc. The proposed framework can be considered as an extension to the well-known flow network problem for hypergraphs. The edges are directed and each has been associated with a weight. However, the weight of an output arc of a node is different from the weight when the arc enters another node as an input. Besides, the hypergraph structure complicates pathways topology [37]. The complication arises since fluxes of contributing input and output reactions over a node may be related. This results to the non-linear property of fluxes. The non-linearity has been considered in our flux calculation procedure (see Eq. 2 and Eq. 3).

The first introduced double-description method for finding EFMs explores the whole set of pathways to find pathways that are in the steady-state solution space and then selects EFMs by comparing the reaction subset of pathways [5]. In the improved double-description versions of the method which use null-space of the stoichiometry matrix as an input, different combinations of pathways in the null-space are calculated and then the reaction subset of pathways are compared [12, 13, 14]. The main difference of the GB-EFM method with double-description-based methods is that GB-EFM first calculates the shortest pathways in the solution space and then checks to see if these pathways can be in the steady-state condition by using some rules on the topology of the pathway.

The method of [20] explores the shortest pathways in the AND/OR graph of the network. The authors did not take into account the stoichiometry coefficients. Thus, the resulting pathways may or may not be EFM. In the method of [21], all reaction combinations to balance a metabolite are calculated and all EFMs are computed accordingly. The order of traversing unbalanced metabolites is based on the graph model. The proposed model in this paper traverses the graph and keeps both the stoichiometry information and the minimality of the pathway with the opportunity of not exploring the whole solution space. The methods based on linear programming usually restrict their method to integer variables [18] which is limiting. In the case of the proposed method this is not a serious concern.

The presented approach, as discussed earlier, aims at finding EFMs by specific characteristics, i.e., the third category is Section 4.2. This brings up the fact that we are looking for new ways to reduce the solution space. EFMtool and the method of [21] target the whole set of EFMs, which means either they find the complete set of EFMs of a metabolic network (which is a lot of non-categorized output information) or they fail to find any EFM. In comparison, our method is looking for a meaningful subset of EFMs which leads to a categorized output set with the cost of losing speed.

In several approaches such as [16, 18], reversible reactions are considered as two irreversible reactions and the two-cycle EFMs are being removed during the EFM calculation. As it has been studied in [38], by splitting reversible reactions, the size of the flux cone description increases by one. However, in the proposed method, while reversible reactions are added to both input and output of a metabolite, they are treated as reversible routes. When a reversible route is used in one direction, the other direction is disabled. Consequently, the proposed method preserves the dimension of the flux cone.

5 Conclusion and Future Work

In order to calculate the elementary flux modes of a given metabolic network, an algorithm was proposed based on a modified AND/OR graph. It was shown in the paper that the steps in the algorithm leads to exact pathways according to EFM definition. The worst-case computational complexity of the algorithm was calculated as $O({r_{\rm{max}}}^{M}\times M^{2})$ which is a newly reported upper-bound for the number of EFMs with external metabolites. Additionally, a set of test cases was provided to prove the concept of the algorithm.

The main focus of this paper was to introduce a model to find minimal flux-balanced pathways in metabolic networks. The model can also be applicable to several constraint-based pathway analysis approaches. Using the potential parallelism in the method to speed-up the algorithm on a hardware structure and proposing different categorization of EFMs based on reactions/metabolites characteristics, as is possible according to the graph data model, are considered as our future research.

Acknowledgements

Authors would like to thank Dr. Nathan Lewis (UCSD) for helpful discussions on the concept.