Effective Identification of Conserved Pathways in Biological Networks Using Hidden Markov Models

Background The advent of various high-throughput experimental techniques for measuring molecular interactions has enabled the systematic study of biological interactions on a global scale. Since biological processes are carried out by elaborate collaborations of numerous molecules that give rise to a complex network of molecular interactions, comparative analysis of these biological networks can bring important insights into the functional organization and regulatory mechanisms of biological systems. Methodology/Principal Findings In this paper, we present an effective framework for identifying common interaction patterns in the biological networks of different organisms based on hidden Markov models (HMMs). Given two or more networks, our method efficiently finds the top matching paths in the respective networks, where the matching paths may contain a flexible number of consecutive insertions and deletions. Conclusions/Significance Based on several protein-protein interaction (PPI) networks obtained from the Database of Interacting Proteins (DIP) and other public databases, we demonstrate that our method is able to detect biologically significant pathways that are conserved across different organisms. Our algorithm has a polynomial complexity that grows linearly with the size of the aligned paths. This enables the search for very long paths with more than 10 nodes within a few minutes on a desktop computer. The software program that implements this algorithm is available upon request from the authors.


Introduction
Recent advances in high-throughput experimental techniques for measuring molecular interactions [1][2][3][4] have enabled the systematic study of biological interactions on a global scale for an increasing number of organisms [5]. Genome-scale interaction networks provide invaluable resources for investigating the functional organization of cells and understanding their regulatory mechanisms. Biological networks can be conveniently represented as graphs, in which the nodes represent the basic entities in a given network and the edges indicate the interactions between them. Network alignment provides an effective means for comparing the networks of different organisms by aligning these graphs and finding their common substructures. This can facilitate the discovery of conserved functional modules and ultimately help us study their functions and the detailed molecular mechanisms that contribute to these functions. For this reason, there have been growing efforts to develop efficient network alignment algorithms that can effectively detect conserved interaction patterns in various biological networks, including protein-protein interaction (PPI) networks [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], metabolic networks [7,12,21], gene regulatory networks [22], and signal transduction networks [23]. It has been demonstrated that network alignment algorithms can detect many known biological pathways and also make statistically significant predictions of novel pathways.
Network alignment can be broadly divided into two categories, namely, global alignment, which tries to find the best coherent mapping between nodes in different networks that covers all nodes; and local alignment, which simply tries to detect significant common substructures in the given networks. Typically, the global network alignment problem is formulated as a graph matching problem whose goal is to find the optimal alignment that maximizes a global objective function that simultaneously measures the similarity between the constituent nodes and also between their interaction patterns. This optimization problem can be solved by a number of techniques, such as integer programming [24], spectral clustering [16,17], and message passing [20]. To cope with the high complexity of the global alignment problem, many algorithms incorporate heuristic techniques, such as greedy extension of high scoring subnetwork alignments and progressive construction of multiple network alignments [9,15,17,19].
There are also many local network alignment algorithms, where examples include PathBLAST [6], NetworkBLAST [10], QPath [11], PathMatch and GraphMatch [12], just to name a few. These algorithms can effectively find conserved substructures with relatively small sizes, but many of them suffer from high computational complexity that makes it difficult to find larger substructures. Furthermore, many algorithms have limited flexibility of handling node insertions and deletions and/or rely on randomized heuristics that may not necessarily yield optimal results. In [18], we introduced an effective framework for local network alignment based on hidden Markov models (HMMs) that can effectively overcome many of these issues. The HMM framework can naturally integrate both the ''node similarity'' (typically estimated by sequence similarity) and the ''interaction reliability'' into the scoring scheme for comparing aligned paths, and it can deal with a large class of path isomorphism. Based on the HMM-based framework, we devised an efficient algorithm that can find the optimal homologous pathway for a given query pathway in a PPI network, whose complexity is linear with respect to the network size and the query length, making it applicable to search for long pathways. It was demonstrated that the algorithm can accurately detect homologous pathways that are biologically significant. However, the algorithm in [18] was mainly developed for querying pathways in a target network, hence it cannot be directly used for local alignment of general networks.
In this paper, we extend the HMM-based framework proposed in [18] to make it applicable for local alignment of general biological networks. Especially, we focus on the problem of identifying similar pathways that are conserved across two or more biological networks. Based on HMMs, we propose a general probabilistic framework for scoring pathway alignments and present an efficient search algorithm that can find the top k alignments of homologous pathways with the highest scores. The algorithm has polynomial complexity which increases linearly with the length of the aligned pathways as well as the number of interactions in each network. The aligned pathways in a predicted alignment may contain flexible number of consecutive insertions and/or deletions. By combining the high-scoring pathway alignments that overlap with another, we can also detect conserved subnetworks with a general structure. Note that the algorithm can be also used for network querying, by designating one network as the query and another network as the target network.

Methods
In this section, we present an algorithm for solving the local network alignment problem based on HMMs. For simplicity, we first focus on the problem of aligning two networks, which can be formally defined as follows: Given two biological networks G 1 and G 2 and a specified length L, find the most similar pair (p,q) of linear paths, where p belongs to the network G 1 and q belongs to G 2 , and each of them have L nodes. As we show later, the pairwise network alignment algorithm can be easily extended for aligning multiple networks in a straightforward manner.

Pairwise Network Alignment
Let G 1~( U,D) be a graph representing a biological network. We assume that G 1 has a set U~fu 1 ,u 2 , . . . ,u N1 g of N 1 nodes, representing the entities in the network, and a set D~fd ij g of M 1 edges, where d ij represents the interaction (binding or regulation) between u i and u j . When the network G 1 is undirected, we assume that both d ij and d ji are present in the set D for simplicity. For example, when G 1 represents a PPI network, u i corresponds to a protein, and the edge between u i and u j indicates that these proteins can bind to each other. For a pair (u i ,u j ) of interacting nodes such that d ij [D, we define their interaction reliability as w 1 (u i ,u j ). Similarly, let G 2~( V,E) be another graph with N 2 nodes and M 2 edges, representing a different biological network. We denote the interaction reliability between two nodes v i and v j in the graph G 2 as w 2 (v i ,v j ). Finally, we denote the similarity between two nodes u i [G 1 and v j [G 2 in the respective networks as h(u i ,v j ), which may be derived using the sequence similarity between two biological entities represented by two nodes as in our experiments.
Our goal is to find the best matching pair of paths p~p 1 p 2 . . . p L (p i [U) and q~q 1 q 2 . . . q L (q i [V) in the respective networks that maximizes a predefined pathway alignment score S(p,q). In order to obtain meaningful results, the alignment score S(p,q) should sensibly integrate the similarity score h(p i ,q i ) between aligned nodes p i and q i (1ƒiƒL), the interaction reliability scores w 1 (p i ,p iz1 ) between p i and p iz1 (1ƒiƒL{1) and w 2 (q j ,q jz1 ) between q j and q jz1 (1ƒjƒL{1), and the penalty for any gaps in the alignment. Figure 1C illustrates an example of an alignment between two similar paths p and q, where p belongs to G 1 and q belongs to G 2 as shown in Fig. 1A. The dashed lines in Fig. 1A that connect two nodes u i and v j indicate that there exist significant similarities between the connected nodes. In the example shown in Figure 1C, the optimal alignment that maximizes the alignment score S(p,q) has two gaps at q 3 and p 5 . Note that ''insertions'' and ''deletions'' are relative terms, and an insertion in p (e.g., p 5 ) can be viewed as a deletion in the aligned path q, and similarly, an insertion in q (e.g., q 3 ) can be viewed as a deletion in p.

Network Representation by HMM
To define the alignment score S(p,q), we adopt the hidden Markov model (HMM) formalism. We begin by constructing two HMMs based on the network graphs G 1 and G 2 . Let us first focus on the construction of HMM for G 1 . Each node u i [U in G 1 corresponds to a hidden state in the HMM. For convenience, we represent this hidden state using the same notation u i . For each edge d ij [D in the graph G 1 , we add an edge from state u i to state u j in the HMM. The resulting HMM has an identical structure as the network graph G 1 . The HMM for G 2 can be constructed in a similar way. Figure 2A illustrates the HMMs that correspond to the network graphs shown in Fig. 1A. In order to find the best matching pairs of paths in the given networks, we define the concept of a ''virtual'' path s~s 1 s 2 . . . s L that contains L nodes, as shown in Fig. 1B. A node s i in the virtual path can be viewed as a symbol that is emitted by a pair of hidden states p i and q i in the respective HMMs. From this point of view, the two HMMs can be regarded as generative models that jointly produce (or ''emit'') the virtual path s, and the underlying state sequence for s will be a pair of state sequences p and q in the respective HMMs. Therefore, the concept of a virtual path can naturally couple a path in G 1 with another in G 2 , providing a convenient framework for identifying conserved pathways in the original biological networks.
The described HMM-based network representation allows us to naturally integrate the interaction reliability scores and the node similarity scores into an effective probabilistic framework. We first define two mappings f 1 : w 1 (u m ,u n ).t 1 (u n ju m ) and f 2 : w 2 (v m ,v n ).t 2 (v n jv m ), which convert the interaction reliability scores w 1 (u m ,u n ) and w 2 (v m ,v n ) between two nodes in G 1 and G 2 to the following transition probabilities P(p i~un jp i{1~um )~t 1 (u n ju m )~f 1 (w 1 (u m ,u n )) ð1Þ between the corresponding hidden states in the constructed HMMs.
Similarly, the mapping f 2 follows the same constraints: To specify the emission probability of a virtual symbol s i at a pair of hidden states in the two HMMs, we define another mapping g : h(u m ,v n ).e(u m ,v n ) that converts the node similarity score h(u m ,v n ) to the following ''pairing'' probability where (p i ,q i )~(u m ,v n ) is the pair of underlying hidden states for s i . The mapping g is defined so that (i)

Ungapped Alignment
Based on the HMM framework, the problem of finding the best matching pair of paths is transformed into the problem of finding the optimal pair of state sequences in the two HMMs that jointly maximize the observation probability of the virtual path s. In an ungapped pathway alignment, the underlying state pair (p i ,q i ) of a virtual symbol s i directly corresponds to a pair of aligned nodes in the original networks. We can find the optimal pair of paths in polynomial time by using a dynamic programming algorithm defined in the following, which is conceptually similar to the Viterbi algorithm. We first define c(t,j,') as the log-probability of the most probable pair of paths for a subsequence b s s~s 1 . . . s t of length t (ƒL), where the underlying states for the virtual symbol s t are p t~uj and q t~v' . The log-probability c(t,j,') can be recursively computed as follows: We repeat the above iterations until t~L. At the end of the iterations, the maximum log-probability of the virtual path s is given by: where fp Ã ,q Ã g~arg max p,q ½log P(p,q) is the optimal pair of state sequences that correspond to the best matching paths in the original biological networks. Once we have computed  log P(p Ã ,q Ã ), it is straightforward to find fp Ã ,q Ã g by tracing the recursive equations that led to the maximum log-probability log P(p Ã ,q Ã ). Although the above algorithm only finds the topscoring pair of paths, we can easily extend it to find the top k pairs simply by replacing the max operator by an operator that finds the k largest scores. The computational complexity of the above algorithm is O(kLM 1 M 2 ) for finding the top k pairs of matching paths, where L is the length of the aligned paths that we want to find, M 1 is the number of edges in G 1 , and M 2 is the number of edges in G 2 . Note that the complexity is linear with respect to all the parameters k, L, M 1 , and M 2 .
The log-probability S(p,q)~log P(p,q) can serve as a good alignment score for the paths p and q that effectively combines node similarity and interaction reliability. In principle, we can also use non-stochastic emission (pairing) scores s em (u j ,v ' ) and transition scores s 1 tr (u j ju i ) and s 2 tr (v ' jv k ) in the recursive equation (4), in place of the log-probabilities log e(u j ,v ' ), log t 1 (u j ju i ), and log t 2 (v ' jv k ), respectively. This will yield a non-stochastic pathway alignment score instead of an observation probability.
As we can see, the concept of the ''virtual'' path provides an intuitive way of coupling states in two different HMMs. In fact, by taking a closer look at the recursive equation (4), the proposed alignment algorithm can also be viewed as a Markovian walk on a product graph, whose nodes consist of all possible pairs of hidden states in the respective HMMs and the edges between these nodes are determined by the connectivity (or transition probability) between the corresponding states in the HMMs. The algorithm searches for the optimal path (or the top-k paths) in the product graph that yields the highest score based on the parameters of the given HMMs.

Alignment with Gaps
To accommodate gaps in the aligned paths p and q, we modify the previous HMMs as follows. First, we add an accompanying stateũ u m for every state u m in G 1 , and similarly, we add an accompanying stateṽ v n for every state v n in G 2 . Next, we add an outgoing edge from each state to the corresponding accompanying state. In addition to this, we also add outgoing edges from the accompanying state to all the neighboring states of the original state. To be more precise,ũ u m will have an outgoing edge to every u k [U(m)~fu k jd mk [Dg, andṽ v n will have an outgoing edge to every v ' [V(n)~fv ' je n' [Eg. By varying the transition probabilities t 1 (ũ u m ju m ) and t 2 (ṽ v n jv n ), we can control the probabilities of having insertions and/or deletions, and thereby control the ''gap penalties'' in a pathway alignment. We adjust the outgoing transition probability from u m so that t 1 (ũ u m ju m )z P uk t(u k ju m )1 ; and for the outgoing transition probability from v n so that We can also control the probabilities of having consecutive insertions or deletions by adjusting the probabilities t 1 (ũ u m jũ u m ) and t 2 (ṽ v n jṽ v n ) for making self-transitions at eitherũ u m orṽ v n . The outgoing transition probabilities t 1 (u k jũ u m ) from an accompanying stateũ u m are chosen so that they are proportional to t 1 (u k ju m ) and satisfy t 1 (ũ u m jũ u m )z P u k t 1 (u k jũ u m )1 . The transition probabilities in G 2 can be chosen in a similar manner. The structures of the modified HMMs are depicted in Fig. 2B. Note that, in a gapped alignment, the matching paths (or state sequences) p and q will still contain L nodes each, and the only difference from an ungapped alignment is that the paths may now contain one or more accompanying nodes which represent gaps. The proposed framework does not impose any restriction on the number of gaps and their locations in the pathway alignment.
In order to find the optimal pair of paths (and their alignment) that maximize the pathway alignment score, we can apply the same dynamic programming algorithm described in the previous section. The retrieved paths can contain any of the hidden states u j (1ƒjƒ2N 1 ) and v ' (1ƒ'ƒ2N 2 ) in the modified HMMs, where we define u mzN1Dũ u m and v nzN2Dṽ v n for notational convenience. The optimal paths fp Ã ,q Ã g~arg max p,q ½log P(p,q) is the best matching pair of paths from two networks, and they may now contain insertions and/or deletions. As before, if we want to find the top k pairs instead of a single top-scoring pair, we can simply replace the max operator by an operator that finds the k largest scores. Note that the computational complexity of the algorithm is O (4kLM 1 M 2 ), which is still linear with respect to all the parameters.

Extension to Multiple Networks
It is straightforward to extend the described pairwise network alignment algorithm for aligning multiple networks. Without loss of generality, we only consider the extension to the alignment of three networks. Given three network graphs G 1 , G 2 , and G 3 , we construct the corresponding HMMs based on their structures. We again use the concept of virtual paths, and now we assume that a virtual path s is jointly emitted by these three HMMs. The emission of a virtual symbol s i is now governed by a pairing probability e(u j ,v ' ,x n ) of three hidden states u j , v ' , and x n that belong to the HMMs that correspond to G 1 , G 2 , and G 3 , respectively. We can find the best matching paths based on the following recursive equation: where e(u j ,v ' ,x n )!e(u j ,v ' )e(v ' ,x n )e(u j ,x n ) is assumed for simplicity. We repeat the above iterations until we reach t~L and compute the maximum log-probability as follows: log P(p Ã ,q Ã ,r Ã )~max p,q,r log P(p,q,r) ½ max j,',n c(L,j,',n), ð7Þ where fp Ã ,q Ã ,r Ã g~arg max q,q,r ½log P(p,q,r) corresponds to the set of best matching paths in the three networks.

Implementation of the Alignment Algorithm
It should be noted that although we fix the length of the virtual path to L, we can in fact find any top-scoring alignment with a shorter length L 0 ƒL, since we store all the alignment scores for shorter alignments while running the dynamic programming algorithm. The recursive equations in (4) and (6) do not restrict multiple occurrence of the same node in the final pathway alignment. However, when it is desirable to avoid such multiple occurrence, we can easily incorporate a ''look-back'' step into each iteration in order to prevent adding a node that is already included in the (intermediate) alignment. As this requires tracing the intermediate optimal (or top k) alignment, the computational complexity of the recursive equations (4) and (6) with a ''lookback'' step will be increased in proportion to the length of the intermediate alignment.
In order to obtain more general subnetwork alignments, not just alignments of linear paths, we can combine the overlapping paths among the top k retrieved pairs of paths. The edges that are already contained in the constructed subnetwork alignment (which correspond to the conserved molecular interactions in the biological networks) are then removed from the HMMs, and we run the dynamic programming algorithm again to find another subnetwork alignment that does not overlap with the retrieved subnetworks. By repeating this ''search and peel-off'' process, we can effectively find diverse subnetwork regions that are conserved in the given networks.
The memory complexity of the proposed algorithm is O (kLN 1 N 2 . . . N B ) for finding the top k pathway alignments for B networks. Although the required amount of memory increases only linearly with respect to each parameter, it can still make the algorithm infeasible when we want to align multiple number of large networks. To overcome this problem, we may assign nonzero pairing probabilities e( : ) to a set of nodes (in the respective networks) only if every pair in this set has considerable node similarity that exceeds a certain threshold. Assuming that there are T sets of nodes that satisfy this condition, we only need to consider these T possible node alignments, in which case the overall memory complexity reduces to O(kLT). Since T is often much smaller than N 1 N 2 . . . N B , this scheme can save significant amount of memory, thereby making the algorithm feasible.

Results
To demonstrate the effectiveness of the HMM-based network alignment algorithm, we carried out the following experiments. First, we used our algorithm to align two pairs of small synthetic networks that were used to validate the network alignment algorithm proposed in [24]. Second, we used the proposed algorithm for finding putative pathways in the fruit fly PPI network that look similar to known human pathways. Finally, we applied the algorithm for aligning microbial PPI networks to assess its ability to find conserved functional modules.

Aligning Synthetic Networks
To illustrate the potential capability of aligning different types of molecular networks, we first tested our algorithm using two small synthetic examples, which include a pair of undirected networks and another pair of directed networks. These examples were obtained from the tutorial files in the PathBLAST plugin of software Cytopscape (version 1.1, http://www.cytoscape.org/ plugins1.php) and they were used for the validation of a network alignment algorithm called MNAligner [24].
HMM parameterization. For aligning the synthetic networks, we parameterized the HMMs as follows. We set the transition scores s tr (u n ju m ) directly based on the ''adjacent matrices'' given in [24], which contain the interaction scores between two nodes in the respective networks. Every interaction score takes a value between 0 and 1, hence we can view it as the ''interaction probability''. We took the logarithm of this interaction probability as the transition score s tr (u n ju m ). When there is no interaction between two nodes, we have s tr (u n ju m )~{?. This keeps the HMM from making a direct transition from a state u m to a non-relevant state u n , thereby preventing the inclusion of irrelevant protein interactions that do not have any biological support in the network. Similarly, we obtained the emission scores s em (u m ,v n ) by taking the logarithm of the similarity scores between nodes given by the ''similarity matrices'' in [24]. The adjacent matrices and the similarity scores for the two examples can be found in the Supporting Information S1.
Example 1: Aligning undirected networks. We first used our algorithm for aligning a pair of undirected networks. To compare the alignment results with the results obtained by MNAligner [24], we looked for the top 500 alignments without gaps, where the length of the virtual path was set to L~3. By incorporating ''look-back'' steps into our dynamic programming algorithm, we restricted the multiple occurrence of the same node pair in the obtained pathway alignment. The top-scoring pathway alignment obtained from our algorithm was AjQQ<CjBB<F jHH, which is identical to the optimal alignment identified by both PathBLAST [6] and MNAligner [24]. Unlike PathBLAST, the proposed HMM-based algorithm and the MNAligner both keep the natural order of the nodes in the original networks. We also noticed that the paths A<C<F and QQ<BB<HH can be aligned with several other potential similar paths in the corresponding networks from the top 500 aligned results. After removing the interactions included in the top-scoring alignment, we searched for the next top-scoring alignment. This returned the alignment JjWW <IjDD<LjOO, which was also ranked as the second best alignment by MNAligner [24]. Repeating the experiment after removing this alignment returned BjMM<DjCC<EjZZ as the third best alignment. This is different from the alignment HjAA<GjNN<BjCC that was found by MNAligner, which got a lower score in our experiment. We noted that the alignment HjAA<GjNN<BjCC is not as significant as the three alignments that we found, as H<G<B can be aligned with many other paths with the same alignment score. By repeating the above experiments and combining the pathway alignment results, we obtained the global network alignment illustrated in Fig. 3A, where a bold line represents that the corresponding edges in the respective networks are matched, whereas a thin line indicates a mismatch. These results show that the HMM-based method can effectively identify the top matching paths in different undirected networks, and it yields better results with higher alignment scores integrating both node similarity and interaction probability compared to PathBLAST and MNAligner for this purpose. Example 2: Aligning directed networks. Without any modification, our algorithm can also be used for aligning directed networks. We demonstrate this by using the second example that contains a pair of small directed networks. In this experiment, we set the length of the virtual path to L~6, which is the length of the longest path in these two networks. As there are fewer legitimate paths in these networks, we only looked for the top 20 aligned pairs of paths. The obtained pathway alignments were combined to get the global network alignment shown in Fig. 3B. The alignment results were similar to those obtained by MNAligner [24], except that we found fewer aligned nodes and edges. This is natural since there exist only a few similar pairs of nodes in the given networks (see Supporting Information S1) and as our algorithm focuses on finding the best local alignments instead of a global alignment. Note that, unlike PathBLAST, which finds path alignments based on several heuristics, the proposed algorithm can find the mathematically optimal path alignment for the given networks.

Aligning Annotated Pathways with PPI Networks
HMM parameterization. The proposed algorithm can also be used for identifying putative pathways in a new biological network, which look similar to known pathways. To demonstrate this, we used our algorithm to search for human signaling pathways in the fruit fly PPI network. In order to compare the search results with those of the network querying algorithm in [18], the HMMs were parameterized according to the nonstochastic scoring scheme in [18] as we describe in the following. The transition score s tr (u n ju m ) was set to s tr (u n ju m )~log (1)~0 in the presence of interaction between the proteins that correspond to u n and u m , and it was set to s tr (u n ju m )~log (0)~{? in the absence of any interaction. To allow gaps in alignments, the transition score from a state u m to its accompanying stateũ u m was set to s tr (ũ u m ju m )~0, and we set the self-transition score atũ u m to s tr (ũ u m jũ u m )~0 to allow consecutive gaps. Furthermore, the score for making a transition fromũ u m to a regular state u n was set to s tr (u n jũ u m )~0 for u n [U(m)~fu n jd mn [Dg and s tr (u n jũ u m )~{? for u n 6 [U(m). The emission score s em (u m ,v n ) for two proteins u m and v n in different networks (where the query network is simply a linear path in this case) was computed based on their sequence similarity. For each protein pair (u m ,v n ), we computed its E-value using the PRSS routine in the FASTA package [25,26], which is known to yield more accurate E-values compared to BLASTP [27]. We regarded a protein pair (u m ,v n ) as a ''match'' if its E-value E v (u m ,v n ) was below a threshold l th . Otherwise, we regarded the pair as a ''mismatch'', which implies that the proteins do not bear significant similarity. Based on this criterion, we set the emission score s em (u m ,v n ) as follows: The value D can be viewed as the mismatch penalty, and is selected so that {D%{ log 10 l th . We set the insertion and deletion penalty also to {D. Finally, since two accompanying states cannot be paired with each other, we set s em (ũ u m ,ṽ v n )~{?.
Querying human pathways in the fruit fly PPI network. We first obtained the PPI network of Drosophila melanogaster from the Database of Interacting Proteins (DIP) [28] and constructed the ''target HMM''. Then we constructed a ''query HMM'' for the human hedgehog signaling pathway and another query HMM based on the human MAP kinase pathway. When constructing the query HMMs, we regarded each signaling pathway as a ''directed network'' with a linear structure, instead of a ''sequence of proteins'' as in [18]. The similarity threshold was set to l th~0 :5 and the gap penalty was set to D~12, as in [18]. The constructed query HMMs were then used to search for matching paths in the target HMM. Despite the generality and the different implementation of the proposed algorithm, the top pathways retrieved by the proposed algorithm agree with the predictions in [18], which is the direct consequence of the mathematical optimality of both methods. For the human hedgehog signaling pathway lhh-Ptch-Smo-Stk36-Gli, the topscoring pathway in the D. melanogaster network agreed well with the putative D. melanogaster hedgehog signaling pathway reported in the KEGG database [29]. In fact, the best aligned path in the fruit fly network contained shh-ptc-Smo-fu-ci, which is identical to the core portion of the putative fly hedgehog signaling pathway ( http: //www. genome. jp /dbget-bin /get_ pathway?org _ name = dme&mapno = 04340) in the KEGG database [29]. The query result of the human MAP kinase pathway Egfr-drk-Sos-Ras85D-ph1-Mekk1-ERKA was also biologically significant, and the seven proteins in the retrieved pathway matched exactly with the proteins in the putative fruit fly MAP kinase pathway (http://www.genome.jp/ dbget-bin/get_pathway?org_name = map&mapno = 04010) reported in KEGG. These results compare favorably to the results obtained by one of the state-of-the-art algorithms [11], where they found two identical proteins in the putative fly hedgehog signaling pathway and five proteins in the putative fly MAPK pathway.

Aligning Microbial PPI Networks
In order to validate the accuracy of our algorithm for predicting functional modules that are conserved in different organisms, we performed additional experiments using three microbial PPI networks obtained from [9]. In our experiments, we performed a pairwise alignment between the E. coli and the C. crescentus networks as well as a pairwise alignment between the E. coli and the S. typhimurium networks. We assessed the accuracy of our algorithm based on the consistency of the KEGG ortholog (KO) group annotations [29] of the aligned proteins. In order to measure the consistency of KO group annotations, we computed the specificity of the predictions based on a similar methodology that was used in [14]. To compute this measure, we first remove all the aligned protein pairs that do not have complete KO annotations, and then compute the total number of annotated protein pairs. An annotated protein pair is regarded as being correct if both proteins have the same KO group annotations, and incorrect if the annotations do not agree. The specificity is defined as the ratio of the number of ''correct'' protein pairs among all annotated protein pairs.
For this experiment, the parameters of the HMMs have been chosen as follows. First, the transition scores s tr (u n ju m ) have been obtained by taking the logarithm of the protein interaction probabilities in the microbial networks, which had been assigned by the SRINI algorithm [30]. The emission scores s em (u m ,v n ) have been computed based on the sequence similarity between the proteins u m and v n , as in the previous section, where the protein similarities have been estimated based on the BLASTP hit scores between protein pairs provided in [9].
Based on the constructed HMMs, we used our algorithm to find the top-scoring pathway alignment with gaps. At each iteration, we looked for the top aligned pair of paths, stored the alignment, and removed the interactions included in the alignment from the respective networks for the next iteration. By repeating this iteration, we found 200 high-scoring path alignments. This experiment has been repeated with varying virtual path length: L~6, 12, 18, 24, and 30. In all our experiments, we disallowed multiple occurrence of identical protein pairs and set the gap/ mismatch penalty to D~1:0. For each experiment, we computed the cumulative specificity for the top k alignments, which is given by where c c i is the total number of correctly aligned protein pairs in the top i alignments, and c a i is the total number of annotated protein pairs also in the top i alignments. The result from the pairwise alignment of the E. coli and the C. crescentus networks is shown in Fig. 4A, and the result from the alignment of the E. coli and the S. typhimurium networks is shown in Fig. 4B. As we can see in both Fig. 4A and Fig. 4B, the cumulative specificity cs k generally decreases when we increase the alignment length L. This is expected since the algorithm tends to recruit more protein pairs in the alignment if we increase L. Furthermore, cs k generally decreases if we increase k. This is natural, since alignments with lower scores correspond to less conserved pathways with larger variations. Although it is difficult to directly compare our results with those reported in [14], it is still worth to note that the cumulative specificity (for the top 200 alignments) of the proposed HMM-based algorithm is higher than the specificity of the alignment algorithm Graemlin 2.0 [14], for both pairwise network alignments. These results clearly indicate that our HMM-based algorithm can produce accurate network alignments that are biologically meaningful.

Discussion
In this paper, we proposed an HMM-based network alignment algorithm that can be used for finding conserved pathways in two or more biological networks. The HMM framework and the proposed alignment algorithm has a number of important advantages compared to other existing local network alignment algorithms. First of all, despite its generality, the proposed algorithm is very simple and efficient. In fact, the alignment algorithm based on the proposed HMM framework is a variant of the Viterbi algorithm. As a result, it has a very low polynomial computational complexity, which grows only linearly with respect to the length of the identified pathways and the number of edges in each network. This makes it possible to find conserved pathways with more than 10 nodes in networks with thousands of nodes and tens of thousands of interactions within a few minutes on a personal computer. Furthermore, the HMM-based framework can handle a large class of path isomorphism, which allows us to find pathway alignments with any number of gaps (node insertions and deletions) at arbitrary locations. In addition to this, the proposed framework is very flexible in choosing the scoring scheme for pathway alignments, where different penalties can be used for mismatches, insertions and deletions. We can also assign different penalties for gap opening and gap extension, which can be convenient when comparing networks that are remotely related to each other. Another important advantage of the proposed framework is that it allows us to use an efficient dynamic programming algorithm for finding the mathematically optimal alignment. Considering that many available algorithms rely on heuristics that cannot guarantee the optimality of the obtained solutions, this is certainly a significant merit of the HMM-based approach. Although the mathematical optimality does not guarantee the biological significance of the obtained solution, it can certainly lead to more accurate predictions if combined with a realistic scoring scheme for assessing pathway homology. As demonstrated in our experiments, the proposed algorithm yields accurate and biologically meaningful results both for querying known pathways in the network of another organism and also for finding conserved functional modules in the networks of different organisms. Finally, the HMM-based framework presented in this paper can be extended for aligning multiple networks. While many current multiple network alignment algorithms adopt a progressive approach for comparing multiple networks [9,[14][15][16][17], our HMMbased framework provides a potential way to simultaneously align multiple networks to find the optimal set of conserved pathways with maximum alignment score.
For future research, we plan to evaluate the performance of our HMM-based algorithm more extensively by investigating the consistency of the predicted alignments based on other available functional annotations, including the gene ontology (GO) annotations [31]. It would be also beneficial to develop a more elaborate scoring scheme that integrates additional information, such as the GO annotations and the KO group annotations, to obtain more reliable alignment results. Finally, we are currently working on simultaneous multiple network alignment based on the HMM framework, where the goal is to construct a scalable multiple alignment algorithm that yields network alignments with higher fidelity.