Definitions and notation

We represent a given biological network using a graph denoted with G=(V, E). Here, the set of nodes V denotes the set of interacting molecules, and the set of edges E denotes the interactions among them. In the rest of this paper, we use the term graph to denote a biological network. Here, we focus on undirected graphs. Figure 1 a represents a graph that contains seven nodes and eight edges.

We say that a graph is connected if there is a path between all pairs of its nodes. We say that a graph S=(V _S ,E _S ) is a subgraph of G if V _S ⊆V and E _S ⊆E. In the rest of this paper, we only consider connected subgraphs. Thus, to simplify our terminology, we use the term subgraph instead of connected subgraph. Notice that a subgraph of a given graph can be uniquely determined by the set of edges E _S of that subgraph as all of its nodes are connected.

We say that two subgraphs \(\phantom {\dot {i}\!}S_{1} = (V_{S_{1}}, E_{S_{1}})\) and \(\phantom {\dot {i}\!}S_{2} = (V_{S_{2}}, E_{S_{2}})\) of G are identical if they have the same set of edges. A less constrained association between two subgraphs is isomorphism. Two subgraphs S ₁ and S ₂ are isomorphic if the following condition holds: There exists a bijection \(f: V_{S_{1}} \rightarrow V_{S_{2}}\) such that \(\forall ({u,v}) \in E_{S_{1}}, \iff (f(u), f(v)) \in E_{S_{2}}\).

We say that two subgraphs S ₁ and S ₂ overlap if they share at least one edge (i.e., \(E_{S_{1}} \cap E_{S_{2}} \neq \emptyset \)). In Fig. 1 a, consider the four subgraphs S ₁, S ₂, S ₃, and S ₄ defined by the set of edges {(a,b), (a,c), (b,c), (b,e)}, {(e,f), (f,g), (e,g), (e,d)}, {(e,f), (f,g), (e,g), (b,e)}, and {(b,e), (d,e), (e,f), (e,g)} respectively. S ₁ and S ₂ are disjoint as they do not share any edges. S ₁ and S ₃ overlap as they share the edge (b,e). Similarly S ₂ and S ₃ overlap. All three subgraphs S ₁, S ₂, and S ₃ are isomorphic as they have the same topology. S ₁ and S ₄ are non-isomorphic as they do not satisfy the bijection function defined above.

Notice that isomorphism is a transitive relation. Thus, for a given subgraph S of G, the set of all subgraphs of G which are isomorphic to S defines an equivalence class. We represent the subgraphs in each equivalence class with a graph isomorphic to those in that equivalence class and call it a pattern. Figure 1 c shows the pattern that represents the equivalence class { S ₁, S ₂, S ₃}.

There are alternative definitions of the frequency of a pattern in a given graph. The classical frequency definition is the number of all subgraphs of the target graph which are isomorphic to the given pattern. This definition, also known as the F1 measure [16], counts all the subgraphs regardless of whether they overlap with each other or not. There are two other frequency definitions which avoid overlaps between different subgraphs. F2 measure counts the largest subset of subgraphs in a given equivalence class which do not share any edges with the rest of the subgraphs in that subset. It however allows them to share nodes. F3 measure is more stringent as it requires that no two subgraphs can share a node. Consider the pattern in Fig. 1 c and the target graph in Fig. 1 a. The frequency of this pattern in the target graph according to the F1 measure is three as it has three embeddings ({ S ₁, S ₂, S ₃}). On the other hand F2 is two { S ₁, S ₂}, and F3 is one (S ₁ or S ₂ or S ₃). From here on, we denote the F1, F2, and F3 counts of a motif M in graph G using the notations F1 _G (M), F2 _G (M), and F3 _G (M) respectively.

The downward closure property states that the frequency of a pattern should monotonically decrease as this pattern grows (by inserting new nodes or edges to it). More specifically, consider a function f() that operates on a pattern and returns a real number. Let us denote two patterns with P ₁ and P ₂. We say that the function f() has downward closure property if and only if f(P ₂)≤f(P ₁) for all (P ₁, P ₂) pairs where P ₁ is a subgraph of P ₂.

Under the light of these definitions, next we show that F1 measure is not downward closed. Consider the pattern P ₁ in Fig. 1 b. The frequency of P ₁ is two in the target graph in Fig. 1 a. Now consider the pattern P ₂ in Fig. 1 c which contains P ₁. Although P ₁ is a subgraph of P ₂, the frequency of P ₂ is three in the same graph (i.e., more than that of P ₁). Next, consider the pattern P ₃ in Fig. 1 d. P ₃ contains P ₂, and its frequency is only one (i.e., less than that of P ₂). This example demonstrates that the F1 measure not only fails to monotonically decrease, but it also fluctuates (i.e., its value may go up or down) as we grow the pattern (see [17, 18] for further discussions on this issue).

Unlike the F1 measure, F2 is downward closed. In the following, we formally prove this.

Similarly, we say that F3 measure which also counts non-overlapping embeddings, is also downward closed.

Failure to satisfy the downward closure property has major implications on the correctness of motif identification. Traditional motif identification algorithms often grow a motif starting from an initial motif of a small number of edges (see Section “Summary of existing methods”). Should they employ the F1 measure, these algorithms cannot have an early stopping criteria as they grow motifs. This is because the frequency can go up as we grow motif even when the current motif frequency is low. Next, we formally define the problem considered in this paper.

Problem definition Given an input graph G=(V, E), the number of nodes in the target motif μ, and frequency threshold α, we aim to find all patterns of μ nodes which have frequency at least α in G under the frequency measure F2. The method we develop in this paper can however be easily extended to F3 as well (see Section “Finding MIS: Going from F1 to F2”).

Summary of existing methods

We classify the literature on motif identification and counting, based on the underlying frequency measure. This is because the frequency measure dramatically changes the cost of counting motifs as well as how we can interpret the frequency of the underlying pattern. Most of the existing studies use F1 frequency measure to count the embeddings of a pattern in a given graph (e.g., [19–24]). These methods carry the drawbacks inherent in the F1 measure. First, F1 ignores the fact that different copies of the same motif can overlap due to the nodes and the edges they share. This can lead to artificially massive number of motif embeddings as the same node or edge can participate in multiple embeddings. To understand this better, consider the pattern and the graph in Figs. 1 c and 1 a respectively. F1 counts three copies of the pattern (S ₁, S ₂, and S ₃). Different nodes and edges however contribute to this count at different numbers. The edge (a, b) appears only in S ₁ while (b, e) appears in both S ₁ and S ₃.

Second and more importantly, the F1 measure is not downward closed. This is because as we grow a pattern by including new edges or nodes, its count as computed by F1 is not monotonic; it may decrease, stay the same, or increase. Lack of downward closure property makes it nearly impossible to decide if the motif found is the largest one in size while growing a pattern. Thus, using F2 is essential for the tractability of identifying frequent patterns. We use the F2 measure in this paper. Thus, the studies limited to the F1 measure are out of the scope of this paper.

Several algorithms tackle the problem of finding frequent patterns in multiple graphs. FSG [25] is one of the key methods in this class. These methods, however, do not count the number of occurrences of a pattern in each graph. They rather check if the given pattern appears at least once in each graph. Vanetik et al. [17] also addressed the same problem.

Finding frequent patterns or counting them without overlaps (i.e., using F2 or F3 measures) have received little attention in the literature. One of the existing algorithms in this category is SUBDUE [15]. Flexible Pattern Finder Algorithm (FPF) [16] detects frequent patterns using both F2 and F3. Two algorithms were proposed by Kuramochi and Karypis [26], named hSiGraM, vSiGraM. However, these algorithms are computationally expensive and do not scale to large graphs or motifs. We evaluate SUBDUE and FSG experimentally in Section “Results and discussion”.

Algorithm overview

In this section, we provide an overview of our method for discovering motifs. At the heart of our method lie four unique graph patterns. We call them the basic building patterns for we use them as guide to construct larger motifs of arbitrary sizes and topologies. Figure 2 presents these basic building patterns. We explain why we use these four specific patterns in Section “Joining patterns to find larger patterns” in detail.

Algorithm 1 presents the pseudo-code of our method. We elaborate on each key step of our method in subsequent sections. The algorithm takes a graph G, the number of nodes of the target motif μ, and the minimum acceptable motif frequency as input α. For each of the four basic building patterns, it first locates all subgraphs in G that are isomorphic to that pattern (Line 1). Let us denote the set of instances of the ith pattern (i∈ {1, 2, 3, 4}) with S _i . In each set S _i , it is possible to have overlapping subgraps. It then extracts the maximum set of edge-disjoint subgraphs in each set S _i (Line 2) (see Section “Finding MIS: Going from F1 to F2” for details). Let us denote the resulting set with S i′ for the ith pattern. Notice that the cardinalities of the sets S _i and S i′ are the F ₁ and F ₂ measures of the ith pattern respectively. The union of all the sets S i′ constitutes the current motif instances as well as the basic building pattern instances at this point (Line 3). The algorithm then iteratively grows the current motif set. At each iteration, it joins the current motif set with the basic building pattern set (Line 9). More specifically, a motif instance and a basic building pattern join if they share at least one edge. Joining two such subgraphs either creates a pattern which already exists in the current set (Line 10) or a new pattern (Line 12). At each iteration, after growing the current set, it filters the overlapping subgraphs to identify MIS for each pattern (Line 18). The algorithm removes all patterns with frequency lower than the user supplied cutoff (Line 21). It reports the frequent subgraphs that have as many edges as the target motif size (Line 23). The algorithm terminates when the current set can not be grown to have any other patterns which satisfy the target motif (i.e. each pattern in the current set is either larger than the target motif size or its frequency is lower than the user specified frequency).

Joining patterns to find larger patterns

Here, we describe one join iteration of our method; the process of joining the subgraphs of current set of patterns with the subgraphs of the basic building patterns to construct larger patterns. At the end of the iteration, the resulting set of subgraphs becomes the current set of subgraphs for the next join iteration.

Recall that we join two subgraphs only if they share at least one edge. Joining two such subgraphs either yields a pattern that is isomorphic to one of the existing patterns or a new one. In the former case, we consider the set of subgraphs S isomorphic to that pattern. We check if the new subgraph is already in S. If it is in S, we discard it. Otherwise, we store it in S. In the latter case (i.e., the pattern is observed the first time), we save this as a new pattern and also keep the corresponding subgraph.

Notice that, although the subgraphs in S do not overlap prior to join, this may no longer hold after new subgraphs are inserted into S. At the end of each join iteration, we select the MIS for each pattern. We defer the discussion on how we do this to Section “Finding MIS: Going from F1 to F2”. We then remove the patterns with F2 values below the user supplied frequency threshold, α. This eliminates non-promising patterns, and thus, reduces the number of candidate patterns for the next join iteration. Using the F2 measure ensures that patterns maintain downward closure property. Thus, non-frequent patterns will never grow to yield frequent patterns.

Why do we need different equivalence classes? If the motif frequency is measured using F1, it is sufficient to join the subgraphs belonging to existing patterns with only those which belong to the same equivalence class of the simple pattern with two edges (see Fig. 2 a) to construct any larger pattern. This however is not true when F2 (or F3) is used to count the motif frequency. To understand the rationale behind this, recall that each equivalence class represents a set of disjoint isomorphic subgraphs. As a result, no two subgraphs from the same equivalence class join for they do not share any edges. Therefore we need more than one equivalence class to construct new and larger patterns.

Given that we need multiple patterns, next, we seek the answer to the following question: What is the smallest set of patterns which can be used to produce arbitrary large topologies by joining them? Here we outline the key steps of the proof that the four basic building patterns, presented in Fig. 2, suffice to construct any larger pattern. That said, we do not guarantee to find all copies of such patterns in the target network.

Before we discuss our induction steps, we explain our strategy on a specific motif size of four to improve the clarity of the discussion on induction. Figure 3 shows all the possible patterns which can be constructed with undirected four edges. A careful inspection shows that each one is an overlapping combination of two of the basic building patterns. For instance, the pattern in Fig. 3 a can result from joining the basic pattern in Fig. 2 a with the basic pattern in Fig. 2 c. It is worth noting that we can construct some of the patterns in Fig. 3 by joining two different pairs of basic building patterns. This redundancy ensures we can still locate a specific pattern even if one of those pairs does not exist. Therefore, our method can construct any pattern with four edges from patterns with three or two edges.

We conduct our proof for the arbitrary pattern size by induction.

Basis The four basic patterns in Fig. 2 constitute all possible graph topologies with two or three edges.

Induction step We assume that our method can construct any pattern with up to k edges (k≥ 3). We next show that any pattern with k+1 edges can be constructed by joining a pattern with k edges with one of the basic building patterns.

Recall that the downward closure property states that those smaller patterns have at least as much frequency as the larger one according to F2 (see Theorem 1.1). This means that if a pattern with k+1 edges is frequent, then so is any of the k edge patterns obtained by removing an edge from that pattern.

Consider a graph G and a copy of a pattern P1 of size k edges in G, S ₁. Also, consider a copy of a pattern P2 with k+1 edges such that P2 contains P1 and one additional edge. Let us denote this additional edge with (a, b). We need to show that P2 can be obtained from P1 by joining it with at least one of the basic patterns.

Since both P1 and P2 are connected graphs, at least one of the two nodes a and b has an edge in P1. Without violating the generality of the proof, let us assume that b has an edge (b,c) in P1. Figure 4 a illustrates the two edges (a, b) and (b, c).

First, we consider using the basic pattern M1 in Fig. 2 a in the join operation. In this case, a copy of M1, {(a,b),(b,c)} will join with S ₁ having a common edge (b,c) which will result in the pattern P2 with k+1 edges. This join however occurs only if the subgraph {(a,b),(b,c)} is included in the F2 counts of M1 (i.e. within the chosen non-overlapping copies of M1).

If this condition fails, we consider the degrees of the two nodes b and c in pattern P1. We start with node c. Let us denote the degree of a node with function d e g() (e.g. d e g(c) is the degree of node c in pattern P1).

If d e g(c)>1, then c has at least one more edge on top of (b,c). Let us denote this edge with (c, d) (see Fig. 4 b). In this scenario, we join a copy of the motif M4 (Fig. 2 d), {(a,b),(b,c),(c,d)} (if this copy exists in the F2 count of M4) to obtain P2.

Finally, if d e g(c)=1, it is guaranteed that d e g(b)>1. This is because if both nodes b and c have degree one, S1 cannot be a connected subgraph. Let us denote one of the additional edges of b with (b,d) (see Fig. 4 c). In this case, we join the subgraph that isomorphic to the pattern M3, {(a,b),(b,c),(b,d)}, with S ₁ to obtain P2. We can do this if this copy exists in the F2 count of M3.

In summary, we conclude that any pattern P2 with k+1 edges can be constructed by joining a pattern P1 with k edges (or k−1 edges) and one of the basic building patterns to obtain the additional edge (or edges) if at least one of the many possible scenarios hold. We however cannot guarantee that the joins will find all of the instances of the k+1 edge pattern on the target graph.

Recall that as we aim to calculate the frequency of a given motif using F2, there is no self join of any pattern. Thus, the basic building patterns set is the smallest set of patterns as we can not construct one of those four patterns using the three other patterns. More specifically, this means that we can not use only one of those four basic building patterns to construct larger patterns by joining pairs of subgraphs belong to that pattern’s equivalence class. This is because if we join the embeddings of a single motif topology (such as the first pattern in Fig. 2 a) we cannot get any larger pattern as they do not share any edge(s).

Finding MIS: Going from F1 to F2

Here, we explain how we compute the F2 frequency for a given pattern. We use two algorithms for this purpose. We explain why we have two separate algorithms later in this section after describing the two algorithms. The first one is a heuristic used in the literature [16]. This algorithm constructs a new graph, called the overlap graph for each pattern as follows. Each node in the overlap graph of a pattern denotes an embedding of that pattern in the target graph. We add an edge between two nodes of the overlap graph if the corresponding embeddings represented by those nodes overlap in the original graph. Once the overlap graph is constructed, the algorithm starts by selecting the node with the minimum degree (i.e. overlaps with the minimum number of embeddings) in the overlap graph. We include the subgraph represented by this node in the edge-disjoint set. We then delete that node along with all of its neighboring nodes in the overlap graph. We update the degree of the neighbors of the deleted nodes. We repeat this process of picking the smallest degree node and shrinking the overlap graph until the overlap graph is empty.

The algorithm described above works well for patterns with small number of embeddings. It however becomes computationally impractical as the number of embeddings of the underlying pattern gets large. This is because both constructing the overlap graph (particularly identifying its edges) and updating it are computationally expensive tasks. Therefore, we use this algorithm for all patterns except for the basic building patterns (where number of embeddings are often too large).

The second algorithm addresses the scalability issue of the the first one. This scalability issue is imposed by the expensive task of calculating the degree of each node in the overlap graph (i.e. the number of overlaps of each embedding). Recall from the previous algorithm that this number is considered as a loss value when selecting the node (i.e. embedding) with minimum degree (i.e. number of overlaps) to include in the final MIS of the pattern under consideration. Briefly, the second algorithm we introduce here avoids the expensive task of calculating number of overlaps for each embedding. The algorithm performs this by algebraically computing such numbers instead of performing actual overlapping tests. Once we compute node degrees of the overlap graph, this algorithm selects the disjoint embeddings the same way as the former algorithm described before. More specifically, the algorithm selects the node with the minimum degree and includes its corresponding embedding in the final MIS. It then removes neighboring nodes to that node from the overlap graph. It repeats this process until the overlap graph is empty. Next, we explain how we compute the degree of a node in the overlap graph for the pattern M1 in Fig. 2 a. Our computation is similar for the other three basic building patterns, yet tailored towards their specific topologies (derivation is shown in Additional file 1: Appendix). Figure 5 shows a hypothetical subgraph S ₁={ (a,c), (b,c)} in the input graph G which is isomorphic to M1. This subgraph is represented by a node in the overlap graph of M1’s embeddings. Let us denote the degree of a node in the original graph G with function d() (e.g. d(v _i ) is the degree of node v _i ). Another embedding of M1 in G overlaps with S ₁ only if it contains the edge (a,c), or (b,c). Any edge in G connected to the middle node c forms two overlapping embeddings, one with the subgraph that has edge the (a,c) and the other with the subgraph that has the edge (b,c). We exclude the edges belong to S ₁ (i.e. the embedding we want to calculate its number of overlaps) itself from the potential edges of G that considered in the overlapping embeddings with S ₁. Thus, by excluding the two edges (a,c) and (b,c) from c’s degree, node c yields 2 × (d(c) - 2) overlaps. In addition, any edge that belongs to node a forms an embedding when combined with the edge (a,c). Excluding the edge (a,c), node a yields d(a) - 1 overlaps. Similarly, node b produces d(b) - 1 overlaps. Thus, the total number of overlaps for the embedding S ₁ = { (a,c), (b,c)} combined from edges of its three nodes { (a,b,c)} is

$$ 2(d(c) - 2) + d(a) - 1 + d(b) - 1 = 2d(c) + d(a) + d(b) - 6 $$

Notice that unlike the first algorithm, the second one requires a unique derivation for each pattern. Thus, we apply it only to the basic building patterns, for their topologies do not depend on the input graph. Also, it is worth noting that typically the basic building blocks have much larger number of embeddings as compared to the patterns derived by joining them. Thus, the efficiency of the second algorithm is needed for them more than the patterns obtained in subsequent iterations (see experimental results).

To adapt our method to count non-overlapping embeddings of each pattern according to F3 instead of F2, we only need to change how we calculate the MIS of this pattern. More specifically, we change the criteria which states that “two subgraphs overlap if they share at least one edge” to “two subgraphs overlap if they share at least one node” (see Section “Definitions and notation”). This will result in changing the overlap graph constructed using the first method we explain in this section. In addition, it will also have slight change in calculating the total number of overlap of each embedding using the second method we discuss in this section. Practically, we expect the overlap graph to be denser when we use the F3 measure as compared to that for the F2 measure. To illustrate this, consider the graph G in Fig. 1 a and the pattern in Fig. 1 c. This patter have 3 embeddings in G which are S ₁, S ₂, and S ₃ defined by the set of edges {(a,b), (a,c), (b,c), (b,e)}, {(e,f), (f,g), (e,g), (e,d)}, {(e,f), (f,g), (e,g), (b,e)} respectively. Figure 6 a and Fig. 6 b represent the overlap graph of this pattern based on F2 and F3 measures respectively.

Accelerating our algorithm through efficient filters

Recall that at each iteration, our algorithm generates new subgraphs. For each of these subgraphs, it checks if this subgraph is isomorphic to one of the patterns constructed till that iteration. Isomorphism test is a computationally expensive task. Next, we describe how we avoid a large fraction of these tests.

We develop two canonical labeling strategies for patterns. Canonical labeling assigns unique labels to the nodes of a given pattern [27]. If two patterns are isomorphic, then they have the same canonical labeling. The inverse is however not true. Unlike isomorphism test, comparing the canonical labeling is a trivial task. Following from this observation, when we construct a new subgraph, we first compare its canonical labeling to those of existing patterns. We then limit the costly isomorphism test to only those patterns which have the same canonical labeling as the new subgraph.

The first canonical labeling counts the degree (i.e. number of incident edges) of each node in the given pattern. It then sorts those degrees and keeps them as a vector we call the degree vector. If two patterns have different degree vectors, then they are guaranteed to have different topologies. Despite its simplicity, this labeling filters out a large fraction of patterns. To test its efficiency, we have tested it on random graphs generated using Barabási −Albert model [28]. We generate 1000 pairs of graphs where each pair is non-isomorphic and have the same number of nodes and edges. The degree vector successfully filters 85 % of the 1000 experiments.

The second canonical labeling extends on the first one. It was first introduced by [29]. Consider a pattern P=(V, E). Let us define the distance between two nodes v _i , v _j ∈V as the number of edges on the shortest path that connects v _i and v _j and denote it with x _ij . Let us define the diameter of P as the maximum distance between any two nodes, and denote it with x. Using this notation, we assign label to node v _i as: \(\sum _{j}^{j \in V} 2^{x-x_{ij}-d(v_{j})}\). Once we compute the labels of all the nodes in the given pattern, we sort them. We call the resulting vector the nodes vector. Similar to the first labeling above, two isomorphic graphs are guaranteed to yield the same labeling. We compute and compare the nodes vector with only the patterns which cannot be eliminated using the first canonical labeling. We then consider the patterns with identical canonical labels for graph isomorphism.

Complexity analysis

Notice that x _i here depends significantly on the topology and the density of the given network G. To the best of our knowledge, there is no closed formula that calculates x _i (i.e. the number of unique topologies of certain size in a given graph G).

Evaluation of running time

In this experiment, we evaluate the running time of our motif discovery algorithm. Our goal here is to observe the effect of varying parameters; graph size, graph density, and motif size on the running time of our algorithm.

Effect of graph and motif size

We evaluate the running time of our method under varying graph and motif sizes using both synthetic and real datasets.

Results on synthetic graphs

We generate synthetic graphs of varying size (i.e. number of nodes) from 100 to 1000 at increments of 100. We fix the graph density to two edges per node on the average (i.e., mean node degree is set to four). We set the minimum desired motif frequency, α= 10. We run experiments for motif sizes μ= 5, 10, and 15 and report the running time. Figure 7 presents the results.

The results demonstrate that our method scales well with growing graph and motif sizes. The running time grows with increasing graph and motif sizes, yet it remains practical for very large graphs. For motif sizes of 5 and 10, it runs in only several minutes even for the largest input graph. As the motif size grows, the cost increases. However, our method can identify very large motifs in a little over a day for massive networks. We observe that the motif size has more influence on the performance of our method than the input graph size. This is because the number of alternative motif topologies grow exponentially with the motif size. This is an inherent characteristic of the underlying computational problem. However, even when the motif size is 15 our method remains to have a practical running time.

Results on real graphs

Next, we test our method on real dataset. We set the minimum desired motif frequency, α= 5. We run experiments for motif sizes μ= 5, 10, and 15 and report the running time. Figure 8 presents the results. Similar to the synthetic dataset results, our method scales to large graph and motif sizes on the real dataset. Note that the number of alternative motif topologies grows exponentially with the motif size. Furthermore, the cost of subgraph isomorphiosm also grows exponentially with the motif size. Despite these two major complicating factors, the running time of our method increases only by about an order of magnitude when we increase the motif size by five. Finally, the parallel between these results and those in Fig. 7 suggests that synthetic graphs generated by Barabási −Albert model have similar structural properties as the real PPI graphs.

Effect of graph size and density

Here, we evaluate the effect of varying input graph size and density on the running time of our algorithm. We use synthetic dataset in order to control the graph density in this experiment. We generate synthetic graphs varying network size from 100 to 1000 at increments of 100. We set the desired motif frequency α= 5 and the motif size μ= 10. We vary graph density from one to four which covers broad range of biological networks [34]. For each input graph and density value, we report the total running time. Figure 9 presents the results.

We observe that the running time increases with growing graph density. As the graph density increases, the number of alternative embeddings of a given motif grows as well. This also increases the number of overlapping subgraph pairs, which in turn increases the cost of finding MIS for each pattern to calculate its F2 frequency (see Section “Finding MIS: Going from F1 to F2”). Despite these major complications inherent in the nature of the motif counting problem, our method remains scalable with respect to growing density. These results suggest that our method is reliable and computationally feasible for a broad range of networks with different sizes and densities.

Comparison with existing methods

Here, we compare our method against two methods in the literature which are tailored towards a problem similar to the one considered in this paper, namely SUBDUE and FSG. We measure the running time and accuracy. We compute accuracy in terms of three parameters, the number of unique motifs found, the average frequency per motif in the target graph, and the frequency of the most abundant motif.

Of these two methods, for SUBDUE, we only report the accuracy of the result as we observe that for most datasets and motif sizes, SUBDUE fails to identify motifs (results shown later in this section). For FSG, we only report the running time. This is because FSG finds motifs in multiple graphs, limited to at most one embedding per graph. In other words, it cannot find multiple embeddings of the same motif in a single graph. Therefore, FSG would yield very low accuracy when applied to a single graph. In the rest of the paper, we will refer to our method as MD (Motif Discovery) for simplicity.

Comparison with SUBDUE

In this experiment, we analyze the effect of varying input graph and motif sizes on the accuracy of our method as compared to that of SUBDUE. We use real dataset in this experiment (see Table 1). SUBDUE does not allow the user to set a minimum allowable motif frequency parameter. It finds all subgraph topologies of a given size even for those subgraphs that appear only once. Due to this limitation of SUBDUE, to have a fair comparison, we set α= 1 for our method as well. We follow our earlier definition (see “Definitions and notation”), and use motif size μ to denote the number of nodes in the given motif topology. We run both methods on each input graph using motif sizes μ= 5, 10, and 15. We report the accuracy of our method as well as SUBDUE. Figures 10, 11, and 12 present the results of μ= 5, 10, and 15 respectively.

Our results for μ= 5 (Fig. 10) demonstrate that both methods identify similar number of unique motifs, yet our method outperforms SUBDUE significantly in terms of the average frequency per motif in all cases (see Fig. 10 b). When we focus on the most abundant topology of each method, we observe a similar pattern; our method always finds patterns with much higher frequency than SUBDUE in all the experiments (see Fig. 10 c). It is worth nothing that motif discovery problem gets exponentially harder with growing motif size. As a result, we expect most algorithms tailored for motif identification to perform well for small motif sizes such as μ= 5. Next, we observe how our method and SUBDUE perform for large values of μ.

As we grow the motif size to μ= 10 (Fig. 11), the results suggest that the gap between our method and SUBDUE grows rapidly in terms all three accuracy measures. More importantly, the results also show that in half of the cases, particularity where the input graph size is large, SUBDUE could not find any motifs while our method continue to locate patterns with high frequency. For example, our method was capable of finding motif topologies with frequency over 100 while SUBDUE could not locate any motif (see Fig. 11 c).

For few cases (see Fig. 11 b), (hhv-8, cje, and tpa), the average frequency per motif of SUBDUE is slightly higher than that of our method. This is because, we set the minimum frequency α= 1. Our method locates many topologies which exist only once while SUBDUE fails to locate them. For example, our algorithm finds thousands of unique motif topologies while subdue outputs only 8 motif topologies for the hhv-8 organism (see Fig. 11 a). As a result, these unique topologies pull the average frequency down. That said, Fig. 11 c confirms that our method can identify motifs which are more frequent than those found by SUBDUE even for those organisms.

As we further increase the motif size to μ= 15 (Fig. 12), the significance of our method becomes more prevalent. We observe that SUBDUE could not find any motifs in any of the graphs accept for tpa’s PPI network. On the other hand, our algorithm not only identifies a massive number of patterns (see Fig. 12 a), but also some of these patterns have very large frequencies (see Fig. 12 c).

In summary, the results demonstrate that our method scales to large input graph and motif sizes and continue to locate patterns with high frequency for a broad range of motif and input graph sizes while SUBDUE fails to do so.

Comparison with FSG

In this experiment, we compare the effect of different input graph and motif sizes to the running time of our algorithm and that of FSG. We use real dataset in this experiment (see Table 1). FSG method requires multiple graphs as input. It defines the frequency of the motif topology as number of different graphs that this motif appears within. Since our method operate on one input graph, we set the desired motif frequency α= 1 to be consistent with FSG. FSG defines motif size as the number of edges in the given motif. To be consistent with FSG, we use μ to denote the number of edges in the motif in this experiment. We run both methods on each input graph using motif sizes μ= 7, 8, and 9. We report the running time of our method (MD) as well as FSG. We do not run experiments for μ> 9 as FSG fails to scale to large motif sizes unlike our method. Figure 13 presents the results.

We observe that our method (MD) is orders of magnitude faster than FSG, particularly in large motif sizes. The running time of our method increases slowly with both motif size and the graph size. On the other hand, the running time of FSG increases slowly with the input graph size, but very rapidly with the motif size. Only for a few cases of small motif sizes (i.e ≤7 edges) FSG performs better than our method. This is due the overhead of calculating F2 for the basic building patterns where number of overlapped embeddings is huge. That said, the running time difference in those cases are negligible. These results suggest that our method outperforms FSG in terms of running time for a broad range of input real biological networks with different sizes. This performance advantage is further magnified by the fact that our method can find multiple embeddings of each motif while FSG finds only one. The two main reasons behind the fact that our method is significantly faster than FSG is that our method (i) does not calculate the frequency of the each new pattern by locating the copies of this pattern in the network using subgraph isomorphism as FSG does, and (ii) it ensures that every generated pattern exists at least once in the underlying graph.

Evaluation of statistical significance

Table 2 suggests that, for small motif size (i.e. μ= 5), the most abundant motif is not significantly more frequent than other motifs of the same size. However, as motifs get large in size (i.e. μ= 10 and 15), the gap between the frequency of the most abundant motif and the rest of the motifs becomes highly significant. This implies that larger motifs characterize topological properties of PPI networks better than small motifs. This is because when motif size is small different motifs have similar frequency values, and this cannot be statistically different in abundance than each other. On the other hand, for large motif size, although the number of unique motif topologies is large, they vary a lot in their abundances; the most frequent one gets significantly more abundant than the rest.

Table 3 shows that, for most of the PPI network and motif size combinations, the most abundant motif is highly over-represented in the original network compared to random networks. In three cases (Rattus norvegicus, μ= 10 and 15, and Human herpesvirus8, μ= 15), we observe that the most abundant is significantly under-represented. These results demonstrate that the motif abundance in PPI networks is not random for nearly all combinations we tested. Thus, studying these structures has great potential to help understand how these networks function. Among the six PPI networks, Rattus norvegicus stands out to be the one with consistently under-represented or random motif abundance. The PPI of Helicobacter pylori consistently has the most significant motif abundance for all motif sizes. This indicates that the interactions in this network follow a regular pattern repeating themselves at different locations of the network. Finally, notice that the two z-score values reported in Tables 2 and 3 do not follow the same pattern (that is a high z-score according to one measure does not imply a high value for the other). This implies that the frequencies of different motifs (i.e., including the ones which are not most abundant) in these PPIs differ from those in random networks. In other words, the PPI networks topologically deviate from random networks.

Case study on Human herpesvirus

Here we briefly analyze the motifs identified by our method on the hhv-8 PPI network which causes Kaposi’s sarcoma disease. We choose this organism in our case study as it has the smallest PPI network among the organisms in our database (see Table 1). Notice from Fig. 11 c that despite its small size (48 nodes and 82 edges), hhv-8 has four disjoint embeddings of a very large motif with 10 nodes, covering a significant fraction of its PPI network. This begs the question whether there is a fundamental recurring function that hhv-8 serves and is covered through evolutionary process with high redundancy. Figure 14 presents the structure of those four embeddings. Each row of Table 4 lists the Uniprot ids of the ten proteins that contribute to each of those embeddings. Analysis of these proteins in the Gene Ontology database [37] reveals that three of those four embeddings, each contains two proteins one responsible for viral DNA packaging (O40944 and P88919) and one responsible for virion assembly (P88954). Without either process, no infectious progeny virus could be formed [38]. Several studies use these two processes as targets to identify effective inhibitors. The existence of these two process in each of the three instances reflects the functional importance of the motif topology found. These results suggest that our algorithm can find significant and valuable motifs which can be use to detect key functions governed by the network processes.

O40944	P88947	P88935	P88951	P88960	P88940	P90489	P88918	P90495	P88902
O40910	O40944	P88947	P88929	P88920	P88925	P88927	P90486	P88918	P88954
P88918	P88919	P88929	P88948	P88920	P88950	O36551	P88942	Q98141	P88954
O40944	Q98141	P88920	P88951	P88954	P88947	P88948	P88958	P88939	P88944