Aligning Metabolic Pathways Exploiting Binary Relation of Reactions

Metabolic pathway alignment has been widely used to find one-to-one and/or one-to-many reaction mappings to identify the alternative pathways that have similar functions through different sets of reactions, which has important applications in reconstructing phylogeny and understanding metabolic functions. The existing alignment methods exhaustively search reaction sets, which may become infeasible for large pathways. To address this problem, we present an effective alignment method for accurately extracting reaction mappings between two metabolic pathways. We show that connected relation between reactions can be formalized as binary relation of reactions in metabolic pathways, and the multiplications of zero-one matrices for binary relations of reactions can be accomplished in finite steps. By utilizing the multiplications of zero-one matrices for binary relation of reactions, we efficiently obtain reaction sets in a small number of steps without exhaustive search, and accurately uncover biologically relevant reaction mappings. Furthermore, we introduce a measure of topological similarity of nodes (reactions) by comparing the structural similarity of the k-neighborhood subgraphs of the nodes in aligning metabolic pathways. We employ this similarity metric to improve the accuracy of the alignments. The experimental results on the KEGG database show that when compared with other state-of-the-art methods, in most cases, our method obtains better performance in the node correctness and edge correctness, and the number of the edges of the largest common connected subgraph for one-to-one reaction mappings, and the number of correct one-to-many reaction mappings. Our method is scalable in finding more reaction mappings with better biological relevance in large metabolic pathways.


Introduction
In the last few decades, the quantity and quality of metabolic data in biological databases such as KEGG (Kyoto Encyclopedia of Genes and Genomes) [1] and Metacyc [2] are greatly increased. The comparative analysis of this vast quantity of metabolic data provides insights into biology and life science applications [3]. An effective way of such analysis is to find the similarity between metabolic pathways by aligning them. The similarity between two pathways is often modeled as a function of the similarity between the aligned nodes or matching edges [4]. By comparing the similarity between metabolic pathways, we can reconstruct phylogeny and infer unknown function or evolution of pathways [5], reveal similar connecting pattern of metabolic pathways [6], and study the structural and functional relevance among species [7,8]. The complexity of the pathway alignment problem stems from its close relationship with graph and subgraph isomorphism problems, which are GI (Graph isomorphism)-Complete and NP-Complete respectively [4]. Thus, it may become impractical to find an accurate solution for this problem as the size of the pathways grows. Due to both computational hardness of pathway alignment and the increasing amount of available metabolic data, obtaining topologically and biologically accurate alignments is a challenging task [9].
In the metabolic pathway alignment problem, metabolic pathways are usually represented as directed graphs, where a node denotes a molecule which can be specified as reaction, enzyme, or compound and an edge represents the interactions between molecules. A one-toone mapping between nodes in two metabolic pathways maps a node from one pathway to a node in the other. A one-to-many mapping between the nodes in two metabolic pathways maps a node from one pathway to a connected subgraph of many nodes in the other [10]. The size of a one-to-many mapping is determined by the number of nodes in this mapping. Performing an alignment is often considered as finding one-to-one mappings or one-to-many mappings between molecules in metabolic pathways [10].
The graph isomorphism problem asks to decide whether two given graphs are isomorphic, and the subgraph isomorphism problem asks to decide whether one graph is isomorphic to a subgraph of another [11]. A straightforward method for identifying the similarity between metabolic pathways is to transform metabolic pathway alignment problem into graph-based isomorphism problem. Considerable efforts were devoted to aligning metabolic pathways in this way [3,[11][12][13][14][15][16][17]. For example, Pinter et al. [12] used enzyme graph to describe metabolic pathway and proposed a tree-based pathway search method called MetaPathwayHunter to align the enzyme graphs by using a graph theoretic approach. Although MetaPathwayHunter obtains a high efficiency in the alignments, the pathways are restricted to trees. To alleviate this restriction, Wernicke and Rasche [13] reduced the pathway alignment problem to subgraph homeomorphism problem and presented an alignment tool METAPAT. METAPAT does not restrict the topology of the metabolic networks in the alignments. Given two metabolic networks G P and G H where G P is represented as the pattern network and G H is represented as the host network, METAPAT determines whether G H contains a subgraph that is isomorphic to G P [13]. Owing to the fact that subgraph homeomorphism problem is NP-complete, METAPAT could be computationally hard with the increasing size of the networks. Meanwhile, Yang and Sze [14] proposed two metabolic pathway matching methods Path-Match and GraphMatch. PathMatch reduces the path matching problem to finding the longest weighted path in a directed acyclic graph while GraphMatch reduces the graph matching problem to finding the highest scoring subgraphs in a graph. Both PathMatch and Graph-Match can effectively and accurately extract biologically meaningful pathways, but finding the matching is time consuming owing to the exhaustive search of subgraphs. Although graphbased isomorphism is the most straightforward idea for aligning pathways, the computational complexity of the graph-based isomorphism problem prohibits its practical application because implementation requires tremendous computing resources as the size of the pathways grows.
In addition, some other methods align metabolic pathways by employing dynamic programming [18][19][20]. In such alignment methods, the similarity between two pathways is defined by the sum of both node and edge matching scores in the similarity objective function. Then, the alignment of pathways is solved by maximizing the similarity objective function between two pathways over all feasible combinations. MNAligner [18] is one example of such methods. MNAligner uses the integer quadratic programming to formulate the alignment of two pathways and find conserved patterns between pathways. To align both two and multiple pathways, Tohsato et al. [19] exploited the global alignment algorithm using dynamic programming to find common pattern from pairwise alignment and then extend pairwise alignment to multiple alignment. Tohsato et al.'s methods were successfully applied to pathway analyses of sugar, DNA and amino acid metabolisms. However, dynamic programming methods do not work well for the large pathway alignment problem since solving the large-scale dynamic programming is time consuming.
Although the above-mentioned methods have achieved considerable progress, there still remains a big challenge. Ay et al. [10] reported that the methods which only search for one-toone mappings between molecules could not identify biologically relevant mappings when different organisms perform the same or similar function through a varying number of steps. An example is shown in Fig 1, where both paths transform LL-2,6-diaminopimelate into 2,3,4,5-tetrahydrodipicolinate. The upper path denotes the shortcut used by plants to synthesize L-lysine. Due to the lack of the gene encoding LL-DAP aminotransferase (2.6.1.83) catalyzing reaction R07613, H. sapiens has to employ a three-step process, as shown with the lower path in Fig 1, to accomplish this transformation. The upper and lower paths should be mapped together in a meaningful alignment when the lysine biosynthesis pathways of human and a plant are aligned. However, due to the different number of reactions in these two paths, traditional methods that are restricted to finding one-to-one mappings fail to uncover the mapping in Fig 1. Motivated by this challenge, researchers develop the other type of alignment methods that allows not only one-to-one mappings but also one-to-many mappings between reactions of two metabolic pathways to tackle this problem. Ay et al. [10] proposed for the first time a one-to-many alignment model and an alignment method called SubMAP which searches oneto-many mappings between reactions of two metabolic pathways. SubMAP successfully identifies biologically relevant mappings of alternative subnetworks, and is scalable for metabolic pathways of arbitrary topology. To improve the quality of one-to-many alignments of metabolic pathways, Abaka et al. [21] presented a constrained alignment method CAMPways where its goal is to maximize the topological similarity while satisfying some constraints on homological similarity. However, due to the cost in exhaustive search of reaction sets, these two methods do not work well for finding reaction mappings in large-scale metabolic pathways.
In this work, we study the problem of aligning two metabolic pathways, which is briefly described as follows. To align two given metabolic pathways, we want to find a set of one-toone, one-to-many or many-to-many mappings between reactions, and maximize the sum of the similarity scores of these mappings. The similarity score of such mapping is evaluated as a function of the similarity between the aligned reactions in the mapping (see Section 'Third Stage' for details). Recall that one-to-many or many-to-many mappings between reactions are used to identify the mappings of alternative pathways that have similar or the same functions through different sets of reactions [10]. High similarity score indicates that the corresponding alternative pathways perform similar or the same functions with high probability.
Our work is based on the observation that connected relation between reactions can be formalized as binary relation of reactions in the metabolic pathway. Motivated by this observation, we propose an alignment method called MPBR for aligning a pair of metabolic pathways exploiting binary relation of reactions. We formalize connected relation between reactions as binary relation of reactions in metabolic pathway. We exploit for the first time the multiplications of zero-one matrices of binary relation of reactions in finding reaction sets. We show that the multiplications of zero-one matrices of binary relation of reactions can be completed in finite steps. As a consequence, we efficiently obtain such reaction sets in a small number of steps without the need of exhaustive search. Furthermore, distinguishing from measuring the topological similarity of reactions based on the direct neighbors of the reactions [10] or the conserved edges induced by the pairs of reaction mappings in the alignment [21], we measure the topological similarity of nodes (reactions) by comparing the structural similarity of the kneighborhood subgraphs of the nodes, which helps to improve the accuracy of the alignments due to the use of more topological information of the neighbors of the reactions. Our experimental results on the KEGG database show that when compared with other state-of-the-art methods, in most cases, MPBR obtains better topological and biological quality of the alignments than CAMPways and SubMAP, and accurately returns more biologically relevant reaction mappings.
The rest of the paper is organized as follows. Section 'Method' presents our method MPBR. Section 'Results' shows experimental results. Section 'Conclusions' concludes the paper.

Preliminaries
To start with, we introduce some definitions and notations. A directed graph G p = (V p ,E p ) is used to denote metabolic pathway P. V p = {r 1 ,r 2 ,. . .,r i ,. . .,r k } is the node set of G p and each node r i represents a reaction in P, i = 1,2,. . .,k. E p is the set of directed edges of G p . There is a directed edge (r i , r j )2E p from r i to r j if and only if at least one output compound of r i is an input compound of r j , i = 1,2,. . .,k and j = 1,2,. . .,k. If both r i and r j are reversible, there is also a directed edge (r j , r i )2E p from r j to r i . Similarly, a directed graph G p 0 = (V p 0 ,E p 0 ) is used to denote metabolic pathway P 0 . Fig 2(A) shows a directed graph for the metabolic pathway of lysine biosynthesis.
For G p = (V p ,E p ), let reaction set RS = {r 1 ,r 2 ,. . .,r n } of size n be a subset of V p such that the induced subgraph of the reactions in RS is linearly connected in the underlying graph, n = 1,2,3,. Problem Statement: Given two pathways G p and G p 0 , we aim to find a set of mappings (RS i , RS j 0 ) between RS n and RS m0 in the alignment of G p and G p 0 such that the sum of the similarity scores of the mappings is maximized, i = 1,2,. . .,N and j = 1,2,. . .,M.
In the following, we introduce how to formalize connected relation between reactions as binary relation of reactions in metabolic pathway. A relation between two related elements of two sets is called binary relation [22]. Accordingly, we formalize the binary relation between reactions A and B as the relation that A is connected with B in a metabolic pathway. In the following section, we present our method MPBR.

MPBR method
For a pair of metabolic pathways G p = (V p , E p ) and G p 0 = (V p 0 , E p 0 ), the goal of MPBR is to find the reaction mappings between G p and G p 0 . Without loss of generality, we assume that |V p | | V p 0 |, reaction sets RSV p and RS 0 V p 0 . MPBR consists of three main stages (as shown in Fig  3): (1) Find all reaction set RS of size n for G p , and find all reaction sets RS 0 of size m for G p 0 (as detailed in Subsection 'First Stage'); (2) Construct a similarity matrix B M by computing the similarity between the reactions in G p and G p 0 (as detailed in Subsection 'Second Stage'); (3) Find mapping (RS, RS 0 ) such that the similarity score of mapping (RS, RS 0 ) is maximized (as detailed in Subsection 'Third Stage'). A set RS map of mappings (RS, RS 0 ) is the result for aligning G p and G p 0 . Fig 3 shows an example illustrating the process of aligning a pair of sample pathways.
First Stage: Finding the candidate reaction sets. In this subsection, we discuss how to exploit the multiplications of zero-one matrices for binary relation of reactions to create the set RS n = {RS 1 , RS 2 ,. . ., RS N } in G p , and the set RS m0 = {RS 1 0 , RS 2 0 ,. . ., RS M 0 } in G p 0 respectively. For metabolic pathway G p , there is a path from r 1 to r n if there is a sequence of reactions r 1 , r 2 ,. . .,r n with edges (r 1 , r 2 ), (r 2 , r 3 ),. . ., and (r n-1 , r n ) in G p . Accordingly, we derive theorem 1.
Theorem 1: For reactions r i and r j , there is a path of length n from r i to r j in G p if and only if M p n [i,j] = 1, where n is a positive integer, i = 1,2,. . .,n and j = 1,2,. . .,n. Proof: We will use mathematical induction. There is a path of length one from r i to r j in G p if and only if M p [i,j] = 1, so the theorem is true when n = 1.
Assume that the theorem is true for a positive integer n. There is a path of length n+1 from r i to r j if and only if there is a reaction r k in G p such that there is a path of length one from r i to  [i, j] = 1. That is, there are n reactions in this path. Thus, we can construct reaction set RS of size n containing these n reactions. Finally, we search every path of length n between two reactions in each reaction pair of NS in G p , and create each reaction set RS of size n containing the reactions in each path to construct the set RS n = {RS 1 , RS 2 ,. . ., RS N } in G p . Similarly, we can find each reaction set RS 0 of size m and construct the set Fig 3(B) shows an example of reaction sets of size 3. In this example, we first obtain the reaction pairs (r i , r j ) with M p 2 [i, j] = 1, i = 1,2,3 and j = 1,2,3,4,5. These reaction pairs are (r 1 , r 3 ), (r 2 ,r 1 ), (r 3 ,r 2 ), (r 2 ,r 4 ) and (r 3 ,r 5 ). Next we create a reaction pair set NS = {(r 1 ,r 3 ), (r 2 ,r 1 ), (r 3 , r 2 ), (r 2 ,r 4 ), (r 3 ,r 5 )} by these reaction pairs. Then, we search the paths of length 2 between two reactions in each reaction pair of NS. These paths are Path 1(r 1 !r 2 !r 3 ), Path 2 (r 2 !r 3 !r 1 ), Path 3 (r 3 !r 1 !r 2 ), Path 4 (r 2 !r 3 !r 4 ) and Path 5 (r 3 !r 4 !r 5 ) respectively. Finally we create reaction sets RS 1 = {r 1 ,r 2 ,r 3 }, RS 2 = {r 2 ,r 3 ,r 1 }, RS 3 = {r 3 , r 1 , r 2 }, RS 4 = {r 2 , r 3 , r 4 } and RS 5 = {r 3 , r 4 , r 5 } with the reactions in Path 1, Path 2, Path 3, Path 4, and Path 5 respectively.
Based on theorem 1 and the above searching procedure of reaction sets in G p , we drive the following property. }. In other words, when n- 2S. According to property 1, we need to perform M p n-1 to find reaction set RS of size n in Thus, in this way, we avoid the exhaustive search of reaction sets, and produce candidate reaction sets in finite steps.
Second Stage: Calculating the similarity of reactions. The second stage aims to compute the similarity values between any two reactions in G p and G p 0 , which combines topological and homological similarities of reactions, and construct a |V p |×|V p 0 | similarity matrix B M , where B M [u,v] is the similarity value between nodes (reactions) u and v, u2V p , v2V p 0 . We first introduce how to compute topological similarity of nodes (reactions) in metabolic pathways. Our computation of topological similarity of nodes is based on the following observation. If node u is mapped to node v, then their neighbors in the corresponding graphs should also be similar. From this observation, we measure topological similarity of nodes by comparing the structural similarity of the k-neighborhood subgraphs of nodes. Next, we discuss how to compare the k-neighborhood subgraphs of nodes to compute topological similarity of nodes. N k (u) is defined as the k-neighborhood node set of u in G p and N k (u) is a subset of V p , u= 2N k (u), where k is an integer and k!0. The shortest distance between u and x2N k (u) is defined as the number of edges of the shortest path between u and x, which is not exceeding k.
. .,and dðx jN k ðuÞj Þ sorted in a nonincreasing order, and the degree sequence of the nodes in N k (v) is d(y 1 ), d(y 2 ),. . .,d(y j ),. . .,and dðx jN k ðvÞj Þ sorted in non-increasing order. We compare the k-neighborhood subgraph of u with the k-neighborhood subgraph of v to compute the topological similarity T(u,v) between u and v: where K m is the maximal value of k and is determined by the user, dðy i Þ are the sum of degrees of nodes in N k (u) and N k (v) respectively. The impact of edges on the topological similarity T(u,v) is evaluated by In the following, we discuss how to compute homological similarity between reactions. Since reactions consist of the input and output compounds and enzymes, we measure the homological similarity between reactions by the similarities of these components. Thus the homological similarity Bsim(u,v) between u and v is computed by the following equation: where u e is the enzyme catalyzing reaction u, v e is the enzyme catalyzing reaction v, Esim(u e ,v e ) is the similarity score between enzyme u e and enzyme v e . We employ the enzyme similarity score defined in [16]   Both homological and topological similarities of reactions provide significant information for the alignment of metabolic pathways. We are now ready to define our similarity S(u,v) between u and v, which is computed by the following equation: where σ is a balancing parameter between the weight of T(u,v) and the weight of Bsim(u,v), 0 σ 1.
In the second stage, we use Eq (3) to calculate the similarity values between any two reactions in two pathways, and construct a similarity matrix B M for the reactions using these similarity values. For example, when k = 2, σ = 0.5, for simplicity, we assume that the homological similarities between any two reactions in sample pathways G p and G p 0 are 0.5, a similarity matrix B M for the reactions in G p and G p 0 is shown in Fig 3(C). In summary, we first utilize the multiplications of zero-one matrices for binary relation of reactions to find reaction set RS of size n for G p and reaction set RS 0 of size m for G p 0 . Second, in order to improve the topological and biological accuracy of the alignments for metabolic pathways, we propose a measure of topological similarity of nodes (reactions), which compares the structural similarity of the k-neighborhood subgraphs of the nodes. Then, we measure the similarity between reactions by combining topological and homological similarities of reactions and build a similarity matrix B M between all reactions in two pathways. Finally, we employ a greedy search to find a set of reaction mappings (RS, RS 0 ) where the sum of the similarity scores of each mapping is maximized.

Results
MPBR is implemented in Java, the data and program are available at http://210.36.16. 170:8080/MPBR/MPBR.zip. Currently, CAMPways and SubMAP are the two available alignment softwares that allow one-to-many reaction mappings in the alignment of metabolic pathways. We downloaded CAMPways and SubMAP from http://code.google.com/p/campways/ and http://bioinformatics.cise.ufl.edu/SubMAP.html respectively. In this section, we experimentally compared the performance of MPBR with that of CAMPways and SubMAP, and discussed three sample alignments.
From the metabolic pathways of the KEGG database retrieved and reformatted by Ay et al. [10], Abaka et al. [21] provided a dataset of 11 metabolic pathways to evaluate alignment quality. Each pathway in this dataset corresponds to one of the above metabolisms. Following the state-of-the-art method CAMPways [21], we also evaluate alignment quality using this dataset. The experimental evaluations are divided into the pathway alignments between species within the same domain and the pathway alignments between species from different domains. Similar to CAMPways, Homo sapiens (hsa) and Mus musculus (mmu) are selected as two representative species from the eukaryota domain, while Escherichia coli (eco) and Agrobacterium tumefaciens (atc) are selected as two representative species from the bacteria domain.
By using the method proposed by Abaka et al. [21], we merge all pathways of the above metabolisms into a large pathway for each of these four species. Thus, we totally obtain four large merged pathways, namely hsa-1.12 with 1520 nodes, mmu-1.12 with 1466 nodes, eco-1.12 with 1104 nodes and atc-1.12 with 1127 nodes. We also use these four large merged pathways to evaluate the performance of the alignment methods. Moreover, we use eight real metabolic pathways eco00230, eco00240, hsa00230, hsa00240, atc00230, atc00240, mmu00230 and mmu00240 from these four species as test pathways. These eight metabolic pathways are obtained from the literature [10] and they are represented by eco-1.13, eco-1.14, hsa-1.13, hsa-1.14, atc-1.13, atc-1.14, mmu-1.13 and mmu-1.14 respectively in this paper. S1 Table presents the number of nodes and the number of edges of the pathways used in the experiments.
The experimental comparisons are conducted based on six criteria. Next, we introduce these criteria in detail [ where u and v are the nodes in the first pathway, (u,v) is an edge in the first pathway, E is the edge set of the first pathway, E 1 is the matched edge set of the first and second pathways, g(u) and g(v) are the mapping nodes of u and v in the second pathway respectively, and E 2 is the edge set of the second pathway.
2. Node Correctness (NC) is the percentage of nodes of the first pathway that are aligned to the correct nodes of the second pathway. NC is calculated by the following equation [24]: where u is a node in the first pathway, V 1 is the node set of the first pathway, f is the correct node mapping, and g is the alignment mapping. For the correct node mapping f, we use measurement FGC (functional group conversion category), which was previously used to define the correct mapping between pathways in [21], to judge whether the node mapping is correct. Specifically, FGC category is a part of the RCLASS database [27] of KEGG. The reactions in the KEGG database are classified into hierarchically organized functional group categories [21]. There are eight FGC categorizations of the KEGG hierarchy, and each FGC categorization is divided into five levels. Abaka et al. pointed out that an inter-species alignment of a pair of pathways is considered biologically validated if the alignment maps reaction subsets classified under the same FGC category [21]. The biological relevance of reaction mappings is closely related to the FGC hierarchy of reactions in the mappings. More specifically, a reaction mapping is biologically more significant if it includes more reactions with higher FGC hierarchy under the same FGC category. In the experimental results of the main text, for a fixed level 5 of the hierarchy, a node mapping is called correct if there exists at least one category at the 5 th level of the FGC hierarchy that includes all the reactions in the mapping [21].

The Number of Edges of Largest Common Connected Subgraph (ELCCS)
is the number of the edges of the largest connected subgraph of the first pathway that is isomorphic to a subgraph of the second pathway [24]. ELCCS is used to evaluate the topological accuracy and biological relevance of the alignments. The larger and denser connected subgraphs are biologically more valuable [24].
4. C-1tomany is the number of correct one-to-many reaction mappings between the first pathway and the second one. To describe this measurement, we introduce some notations first. Let X, X 0 denote two species and G X , G X 0 represent their metabolic pathways corresponding to some metabolism 1.m, listed earlier in the text. Let (RS, RS 0 ) be a mapping from an alignment of G X and G X 0 where RS = {r 1 }, RS 0 = {r 1 0 , r 2 0 ,. . .,r x 0 }, and P 1 ,. . .,P i ,. . ., P x be the pathways that include reaction r 1 and are associated with metabolism 1.m in the species X [21].
Following the literatures [10] and [21], we measure the correctness of the one-to-many reaction mappings based on two aspects. On the one hand, as Ay et al. [10] reported that if both alternative pathways can complete the transformations between two given compounds through different reaction sets, then these two pathways are considered to be functionally similar. A correct one-to-many reaction mapping between different pathways should be able to identify the mapping of such alternative pathways [10]. On the other hand, Abaka et al. [21] pointed out that an alignment mapping reactions that belong back to the same original KEGG pathway is considered to be of high quality. Combining these two aspects, a one-to-many reaction mapping (RS, RS 0 ) is called correct if it satisfies the following two conditions: (1) The start reactions in RS and RS 0 share at least one input compound and the end reactions in RS and RS 0 share at least one output compound. (2) Every reaction in RS 0 is included in at least one of the pathways P 1 0 ,. . .,P i 0 ,. . ., P x 0 where each P i 0 is a pathway in metabolism 1.m of species X 0 and corresponds to P i of X [21].
5. CR-1tomany is the ratio of the number of correct one-to-many reaction mappings to the total number of mappings produced by the alignment [21]. CR-1tomany can be used to investigate the percentage of the correct one-to-many reaction mappings in the alignment. Higher values for CR-1tomany indicate that the alignments for one-to-many reaction mapping are of high quality [21].
6. C-manytomany is the number of correct many-to-many reaction mappings between the first pathway and the second one. By reference to C-1tomany, let (RS, RS 0 ) be a manyto-many mapping from an alignment of G X and G X 0 where RS = {r 1 , r 2 ,. . .,r x }, RS 0 = {r 1 0 , r 2 0 ,. . .,r y 0 }, and P 1 ,. . .,P i ,. . .,P x be the pathways that include a reaction in RS and are associated with metabolism 1.m in the species X. Similar to C-1tomany, a many-to-many reaction mapping (RS, RS 0 ) is called correct if it satisfies the following two conditions: (1) The start reactions in RS and RS 0 share at least one input compound and the end reactions in RS and RS 0 share at least one output compound. (2) Every reaction in RS 0 is included in at least one of the pathways P 1 0 ,. . .,P i 0 ,. . .,P x 0 where each P i 0 is a pathway in metabolism 1.m of species X 0 and corresponds to P i of X.
MPBR, CAMPways and SubMAP provide the options of one-to-one alignment and one-tomany alignment. We can perform one-to-one alignment of two pathways to find one-to-one reaction mappings between these two pathways. Similarly, we can perform one-to-many alignment of two pathways to find one-to-many reaction mappings between these two pathways. In the experiments, the one-to-many reaction mappings include 1-to-2 and 1-to-3 reaction mappings.
In this paper, we use EC, NC and ELCCS to measure the quality of one-to-one alignment, and use C-1tomany and CR-1tomany to measure the quality of one-to-many alignment. In addition, we use C-manytomany to evaluate the capability of MPBR for searching many-tomany reaction mappings.
In the experiments, MPBR was run using k = 3 and σ = 0.6, and CAMPways and SubMAP were run using their default parameter settings. MPBR, CAMPways and SubMAP were run on the Sugon 5000A computer system of cluster architecture at Guangxi University, using a single computing node with a quad-core Intel(R) Xeon(R) CPU E5620 @ 2.40GHz and 40GB RAM. The operating system is Linux.
The following five subsections will provide our experimental evaluations on the qualities of the alignment results computed by MPBR, CAMPways and SubMAP respectively. Subsections 'Same-domain One-to-one Alignments' and 'Same-domain One-to-many Alignments' focus on one-to-one alignment and one-to-many alignment between the species within the same domain respectively. Subsections 'Across-domains One-to-one Alignments' and 'Acrossdomains One-to-many Alignments' focus on one-to-one alignment and one-to-many alignment between the species from different domains respectively. Subsection 'Running time and memory requirements' discusses the performance of each method in terms of the running time and memory requirements. Subsection 'Many-to-many Alignments' discusses the experimental results of many-to-many alignments of the pathways. Subsection 'Case study' introduces three sample alignments to show how MPBR, CAMPways and SubMAP can be used to analyze metabolic pathways. The

Same-domain One-to-one Alignments
In this subsection, we discuss the quality of the same-domain one-to-one alignments produced by MPBR and other comparative methods. Tables 1-3 summarize the one-to-one alignment results for the same domain species with respect to distinct performance indices.
As shown in Tables 1-3, over all 28 instances, MPBR has the highest values of EC, NC and ELCCS for 19, 18 and 18 instances respectively, whereas all three methods obtain equal values of EC, NC and ELCCS for 5, 6 and 7 instances respectively. Additionally, MPBR and CAMPways obtain equal values of EC, NC and ELCCS for 4, 4 and 3 instances respectively. These experimental results emphasize that, for the same-domain one-to-one alignment, in most The best performer for the relative item is marked in bold.
cases, our MPBR method outperforms other comparative methods not only in topological accuracy but also in biological relevance of the results. Thanks to the use of structural similarity among the neighbors of reactions, MPBR is able to improve the topological and biological accuracy of the alignments.

Same-domain One-to-many Alignments
In this subsection, we compare the quality of the same-domain one-to-many alignments produced by MPBR and other comparative methods. The values of C-1tomany and CR-1tomany The best performer for the relative item is marked in bold.
doi:10.1371/journal.pone.0168044.t003 Table 4. C-1tomany and CR-1tomany of one-to-many alignment results for eco-atc. The best performer for the relative item is marked in bold. The asterisk "*" denotes that the program is unable to generate a result under our current computing environment.
of the one-to-many alignment results for the same domain species are shown in Tables 4  and 5.
From Tables 4 and 5, we can see that, MPBR performs the best with the highest values of C-1tomany and CR-1tomany in 22 and 15 out of all 28 instances respectively. For 4 instances, all three methods obtained the same values of C-1tomany and CR-1tomany, while the value of C-1tomany of MPBR is lower than CAMPways for 2 instances and is lower than SubMAP for 1 The best performer for the relative item is marked in bold. The asterisk "*" denotes that the program is unable to generate a result under our current computing environment.
doi:10.1371/journal.pone.0168044.t005 instance. The value of CR-1tomany of MPBR is lower than SubMAP for 4 instances, and some values of CR-1tomany of MPBR are lower than CAMPways for 7 instances. This means that, in most cases, MPBR is able to return more correct one-to-many reaction mappings than CAMPways and SubMAP in the same-domain one-to-many alignment. On the other hand, when the size of the pathway becomes large, SubMAP is unable to generate a result for 9 instances under the current computing environment while MPBR and CAMPways are not restricted to the size of the pathway in the same-domain one-to-many alignment.

Across-domains One-to-one Alignments
This subsection discusses the quality of the across-domains one-to-one alignments produced by MPBR and other comparative methods. Tables 6-11 present the one-to-one alignment results for different domain species with respect to distinct performance indices. As can be seen from Tables 6-11, over all 56 instances, MPBR performs better than the other two methods with the highest values of EC, NC and ELCCS for 47, 42 and 51 instances respectively. Some values of EC and NC of MPBR are a bit lower than those of CAMPways for 4 and 3 instances respectively. This demonstrates that, in most cases, MPBR is also capable of achieving better topological accuracy and biological relevance of the alignment results than other two comparative methods in across-domains one-to-one alignment. In addition, from Tables 1-3 and Tables 6-11 we can find that the values of EC, NC and ELCCS of the same-domain alignments are obviously higher than those values of acrossdomains alignments. This is also consistent with the analysis that the biological relevance of the species within the same domain is much stronger [10]. Thus, we can employ the alignments of metabolic pathways to analyze the evolution of species.

Across-domains One-to-many Alignments
In this subsection, we compare the quality of the across-domains one-to-many alignments produced by MPBR and other comparative methods. The values of C-1tomany and CR-1tomany of the one-to-many alignment results for different domain species are shown in Tables 12-15. The best performer for the relative item is marked in bold. The asterisk "*" denotes that the program is unable to generate a result under our current computing environment.
From Tables 12-15, we can see that, MPBR achieves the best values of C-1tomany and CR-1tomany compared to other comparative methods in 44 and 32 out of 56 instances respectively. For only 10 out of 56 instances MPBR fails to be the best. On the other hand, combining Tables  12-15 and S1 Table, we can find that, for the across-domains one-to-many alignment, when the size of the pathway is large enough with thousands of reactions, possibly due to the exhaustive search of reaction sets, CAMPways and SubMAP are unable to generate a result for 34 and 4 instances respectively under the current computing environment. In contrast, for these instances, the values of C-1tomany and CR-1tomany of MPBR are still high enough without being affected by the size of pathway. While the comparative methods suffer from the size of large-scale pathway, MPBR overcomes this problem and returns more correct one-to-many reaction mappings.
The above analysis of C-1tomany and CR-1tomany shows that, in most instances, MPBR also performs better than the other two methods in across-domains one-to-many alignment. In conclusion, the results from subsection 'Same-domain One-to-many Alignments' and subsection 'Across-domains One-to-many Alignments' demonstrate that MPBR is an effective method in retrieving one-to-many reaction mappings in the alignment of metabolic pathways.

Running time and memory requirements
In the experiments of one-to-one and one-to-many alignments, we have tested a total of 168 instances. In some cases, SubMAP and CAMPWays consumed an unusually long time until running out of memory, Table 16 summaries the percentage of the instances that can be solved by MPBR, CAMPWays and SubMAP respectively, and the average running time for the solved instances. In Table 16, PSI represents the percentage of the solved instances of each method, in one-to-one alignment, ART1 denotes the average running time for the 64 instances solved by all three methods and ART2 represents the average running time for the 84 instances solved by MPBR and CAMPWays, and in one-to-many alignment, ART3 denotes the average running time for the 41 instances solved by all three methods, and ART4 represents the average running time for the 80 instances solved by MPBR and CAMPWays. The best performer for the relative item is marked in bold. The asterisk "*" denotes that the program is unable to generate a result under our current computing environment. doi:10.1371/journal.pone.0168044.t015 As can be seen from Table 16, for the one-to-one alignment, although the average running time for the 64 solved instances of SubMAP is shorter than CAMPWays and MPBR, SubMAP failed to solve 20 out of 84 instances since it took an unusually long time until running out of memory in these unsolved instances, whereas both CAMPWays and MPBR solved all the instances. Meanwhile we can see that, for the one-to-one alignment, MPBR consumed less time than CAMPWays for the 84 instances solved by MPBR and CAMPWays.
For the one-to-many alignment, we can observe from Table 16 that, MPBR spent less time for the 41 instances solved by all three methods in comparison to CAMPWays and SubMAP; in addition, compared with CAMPWays, MPBR took less time for the 80 solved instances of MPBR and CAMPWays, and MPBR solved all the 84 instances while both CAMPWays and SubMAP did not.

Many-to-many Alignments
In addition to one-to-many alignments, we can also reveal alternative pathways that have similar functions by finding many-to-many reaction mappings between different pathways. This subsection discusses whether MPBR can accurately find such mappings in the alignment of metabolic pathways. The many-to-many reaction mappings include 2-to-2, 2-to-3 and 3-to-3 reaction mappings. Both CAMPways and SubMAP do not implement the functionality of many-to-many alignment. Table 17 shows the C-manytomany of many-to-many alignment results of MPBR. Both Escherichia coli (eco) and Agrobacterium tumefaciens (atc) are single-cell microorganisms, Homo sapiens (hsa) and Mus musculus (mmu) are complex organisms with cell membranes. Table 17 demonstrates that there are a number of many-to-many reaction mappings between the species among the same domain and among different domains. These results suggest that many-to-many reaction mappings frequently appear in nature. MPBR has the capability in finding many-to-many reaction mappings between different pathways to obtain biologically meaningful alignments. Table. NC of one-to-one alignment results for the fourth level of the FGC hierarchy. (DOC) S3 Table. NC of one-to-one alignment results for the third level of the FGC hierarchy. (DOC) S4 Table. NC of one-to-one alignment results for the second level of the FGC hierarchy. (DOC) S5 Table. NC of one-to-one alignment results for the first level of the FGC hierarchy. (DOC)