Combinatorial algorithm for counting small induced graphs and orbits

Graphlet analysis is an approach to network analysis that is particularly popular in bioinformatics. We show how to set up a system of linear equations that relate the orbit counts and can be used in an algorithm that is significantly faster than the existing approaches based on direct enumeration of graphlets. The approach presented in this paper presents a generalization of the currently fastest method for counting 5-node graphlets in bioinformatics. The algorithm requires existence of a vertex with certain properties; we show that such vertex exists for graphlets of arbitrary size, except for complete graphs and a cycle with four nodes, which are treated separately. Empirical analysis of running time agrees with the theoretical results.


Introduction
Analysis of networks plays a prominent role in various fields, from learning patterns [1] and predicting new links in social networks [2,3], inferring gene functions from protein-protein interaction networks [4] in bioinformatics, to predicting various properties of chemical compounds (mutagenicity, boiling point, anti-cancer activity) [5] from their molecular structure in chemoinformatics. Many methods rely on the concept of node similarity, which is typically defined in a local sense, e.g. two nodes are similar if they share a large number of neighbours. Such definitions are insufficient for detecting the role of the node. A typical social structure includes hubs, followers, adversaries and intermediaries between groups. While local similarity definitions treat the hub and its adjacent nodes as similar, a role-based similarity would consider the hubs as similar disregarding their distance in the graph.
A popular approach in bioinformatics extracts the node's local topology by counting the small connected induced subgraphs (called graphlets) [6], which the node touches, and, when a more detailed picture is required, the node's position (orbit) [7] in those graphs. See the following paragraphs for a more formal definition. Fig 1 illustrates all four-node graphlets and orbits of their nodes. Most applications of graphlet and orbit counts are based on the assumption that the node's local network topology is somehow related to the functionality or some other property of the observed node in the network. Therefore, we can assume that nodes with similar signatures will have similar observed properties. This is the foundation for methods such as clustering of nodes, inference of certain node's properties etc. The nodes often a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 the algorithm's time complexity on sparse graphs is lower by an order of magnitude in comparison with enumeration-based approaches.

Preliminaries
Referring to graphlets, orbits, neighbours, etc., requires some notation which we summarize in Table 1 and use throughout this paper.
Let K = (V K , E K ) be a subgraph of J = (V J , E J ), and let v 2 V K . We will denote J's vertex that corresponds to v by v J . If there are multiple isomorphic embeddings of K in J, v J refers to one of them. Similarly, if S V K , then the corresponding vertices in J are denoted by S J .

Related work
The most basic case of counting induced patterns in graphs is that of counting triangles. Itai and Rodeh [8] showed that this can be done faster than by exhaustive enumeration in O(n 3 ) time. Raising the graph's adjacency matrix A to the third power gives the number of paths of length 3 between pairs of nodes. Element A 3 x;x represents the number of paths of length 3 that start and finish in the node x, which corresponds to the number of triangles that include x. The total number of triangles is then 1 6 P x2G A 3 x;x . Note that the same triangle is counted twice for each of its three nodes. The time complexity of this procedure equals that of multiplying two matrices, which is faster than exhaustive enumeration of triangles in dense graphs. A natural extension of this result is to larger cliques. Nesetril and Poljak [9] studied the problem of detecting a clique of size k in a graph with n nodes. They showed that this problem can be solved faster than with the straight-forward O(n k ) solution. Their approach reduces the original problem to detection of triangles in a graph with O(n k/3 ) nodes. Since we can detect triangles faster than in O(n 3 ) with fast matrix multiplication algorithms, we can also detect cliques of size k faster than O(n k ).
Counting all non-induced subgraphs is as hard as counting all induced subgraphs because they are connected through a system of linear equations. Despite this it is sometimes beneficial to compute induced counts from non-induced ones. Rapid Graphlet Enumerator (RAGE) [10] NðSÞ common neighbours of nodes in the set S & V; NðSÞ ¼ \ v2S NðvÞ number of common neighbours of vertex v, of vertices v 1 , v 2 , . . ., v j , and of vertices from set S, respectively; that is, takes this approach for counting four-node graphlets. Instead of counting induced subgraphs directly, it reconstructs them from counts of non-induced subgraphs. For computing the latter, it uses specifically crafted methods for each of the 6 possible subgraphs (P 4 , claw, C 4 , paw, diamond and K 4 ). The time complexity of counting non-induced cycles and complete graphs is O (e Á d + e 2 ), while counting other subgraphs runs in O(e Á d). However, the run-time of counting cycles and cliques in real-world networks is usually much lower. Some approaches exploit the relations between the numbers of occurrences of induced subgraphs in a graph. Kloks et al. [11] showed how to construct a system of equations that allows computing the number of occurrences of all six possible induced four-node subgraphs if we know the count of any of them. The time complexity of setting up the system equals the time complexity of multiplying two square matrices of size n. Kowaluk et al. [12] generalized the result by Kloks to counting subgraph patterns of arbitrary size. Their solution depends on the size of the independent set in the pattern graph and relies on fast matrix multiplication techniques. They also provide an analysis of their approach on sparse graphs, where they avoid matrix multiplications and derive the time bounds in terms of the number of edges in the graph.
Floderus et al. [13] researched whether some induced subgraphs are easier to count than others as is the case with non-induced subgraphs. For example, we can count non-induced stars with k nodes, P x2V cðxÞ kÀ 1 À Á , in linear time. They conjectured that all induced subgraphs are equally hard to count. They showed that the time complexity in terms of the size of G for counting any pattern graph H on k nodes in graph G is at least as high as counting independent sets on k nodes in terms of the size of G.
Vassilevska and Williams [14] studied the problem of finding and counting individual noninduced subgraphs. Their results depend on the size s of the independent set in the pattern graph and rely on efficient computations of matrix permanents and not on fast matrix multiplication techniques like some other approaches. If we restrict the problem to counting small patterns and therefore treat k and s as small constants, their approach counts a non-induced pattern in O(n k−s+2 ) time. This is an improvement over a simple enumeration when s ! 3. Kowaluk et al. [12] also improved on the result of Vassilevska and Williams when s = 2. Alon et al. [15] developed algorithms for counting non-induced cycles with 3 to 7 nodes in O(n ω ), where ω represents the exponent of matrix multiplication algorithms.
Alon et al. [16] introduced the color-coding technique for finding simple paths and cycles in graphs. Their technique is applicable not just to paths and cycles but also to other patterns with a small treewidth. The authors of [17] used such color-coding approach to approximate a 'treelet' distribution (frequency of non-induced trees) for trees with up to 10 nodes.
Recently, Melckenbeeck et al. [18] published a paper that describes how to generate systems of equations similar to those used in the ORCA algorithm [19] for arbitrarily large graphlets. However, the resulting equations do not satisfy the requirements needed for an efficient counting algorithm. Consider for example the equation o 50 + o 55 = ∑ P 7 (x,a,b,c) (c(a, b, c) − 1). There can be as many as O(ed 2 ) sets of nodes {a, b, c} with a nonzero number of common neighbours, which makes the computation of common neighbours the limiting factor in terms of space and time requirements. The method we present in this paper and the related proofs show how to avoid this issue and construct an efficient algorithm for arbitraty graphlet sizes.

Outline of the proposed algorithm
We will derive a system of linear equations that relate the orbit counts of a fixed node for graphlets with k vertices. The coefficients on the left-hand sides reflect the symmetries in the graphlets and do not depend on the host graph, so they are derived in advance. The right-hand sides are computed as sums over graphlets with k − 1 vertices induced in the host graph H, and the sums include terms that represent the number of common neighbours of certain vertices in the embeddings of graphlets in H.
The resulting system of equations will be triangular and have a rank of jO k j À 1. We can efficiently enumerate the complete graphlet, after which the system of equations for the remaining orbit counts can be solved using integer arithmetic, thus avoiding any numerical errors.

Original contributions
We already presented the original idea of the algorithm in a recent article in Bioinformatics [19], in which we focused on its use in genetics and avoided formal descriptions and analysis. In this paper we 1. present the algorithm more formally; 2. describe a general method for derivation of the system of equations relating the orbit counts; 3. generalize it to induced subgraphs of arbitrary size; in particular, we prove that the system of equations with the properties required for the efficient implementation of the algorithm exists for any k ! 4; 4. provide worst time-complexity analysis and the analysis of the expected time complexity on random graphs; 5. empirically explore the efficiency of the orbit counting algorithm and compare it with the theoretical results.
The remainder of the paper is composed of two parts. In the next section we show a technique for building the system of equations with desired properties, and in the following section we present an algorithm based on them and analyze its time-and space-complexity.

Relations between orbit counts
We will show how to construct linear relations between a chosen orbit count o i and some orbits belonging to graphlets with a larger number of edges. We will illustrate the procedure on figures showing the derivation of the following Eq (1) that relates the count for orbit 59 and counts for orbits 65, 68 and 70.

Derivation of general relations between orbit counts
Let orbit O i appear in a connected simple k-node graphlet G a = (V G , E G ) (a = m(i)). We denote the G a 's node that is in orbit O i by x; if there are multiple such nodes, we pick one. Next, we choose a node y 6 ¼ x, such that G 0 = G a \{y} is still a connected graph; we will impose additional constraints on y later to ensure an efficient implementation of the algorithm. According to our notation, x G 0 is the node in G 0 that corresponds to x in G a ; let O m be its orbit. We label the remaining k − 2 nodes with x 1 , x 2 , . . ., x k−2 (Fig 3).
We now go in the opposite direction: starting with G 0 , we consider its possible extensions to G a . Let E & V G 0 be a set of nodes such that adding a new vertex y connected to all vertices in E yields G a with x in orbit O i (Fig 4). Let E be a set of all such subsets E.
Let G 0H be some particular occurrence of G 0 in H. To count o i for the node x H (the node in H to which x maps), we need to explore the extensions of G 0H to G H a . A necessary (but insufficient) condition to put x H into O i is that the additional node y is a common neighbour of all vertices E H for one of E H 2 E H (with respect to the particular occurrence of G 0 in H). There are at most candidate nodes y; c(E G 0 ) represents the number of neighbours of E that are already in G 0 (i.e. x i ) and cannot be mapped to y. Eq (2) represents the term in the sum in the right side of the  Extensions of the k − 1-node graphlet G 0 to G a . G 7 can be extended to G 24 by attaching y to either x 1 and x 3 or to x 2 and x 3 , hence E ¼ ffx 1 ; x 3 g; fx 2 ; x 3 gg. With respect to Eq (2), x 1 and x 3 (as well as x 2 and x 3 ) have one common neighbour in G 0 , so c(E G 0 ) = 1 and relation. To compute the total orbit count o i for x, we sum Eq (2) over all occurences of G 0 in H (Fig 4). Condition Node y can also be connected to any of the other k − 1 − |E| (E 2 E) nodes in G 0H , resulting in 2 k−1−|E| possible graphlets and orbits for x. The counts for these orbits are summed on the left side of the relation (Fig 5).
While adding the orbit counts on the left-hand side, we need to account for the over-counts, that is, the number of times that Eq (2) counts the same occurrence of G Ã = G m(p) (the graphlet containing the orbit p) within H. G Ã is obtained by extending G 0 with y 2 N(E). The coefficients on the left-hand side thus equals the number of ways in which G 0 can be extended to G Ã with a fixed node x (Fig 6).
In general, we have to consider all induced occurrences of G 0 in G Ã (with a fixed point x), which is the same as considering nodes z 2 V G Ã whose removal results in G 0 with x G 0 in orbit O p . For every such case we increase the coefficient by the number of extensions E 2 E such that node z is connected to the extension nodes, i.e. N(z) E G Ã . The general procedure for relating the orbit count o i with counts of orbits with higher indices is outlined in Algorithm 1. (2)) . The right side of equation sums over.

Algorithm 1 Derive an equation for orbit
. all occurrences of G 0 in the host graph.
. Is z in the same orbit as y . given a fixed point x? f p f p þ jfE 2 E : NðzÞ Egj . By how many extensions?
. Left side is a weighted sum of orbit counts. return equation l = r end function Additional constraints on selection of y In the preceding derivation, the only limitation on selection of vertex y was that the remaining graphlet is still connected. Different choices of y yield different equations. With the coefficients independent of the host graph and known in advance, the time consuming part of using these equations to calculate orbit counts is the computation of the right-hand side terms. To speed it up, we impose some additional constraints on the choice of the node y: the restraints will be such that the right-hand sides will contain only the counts c(S) in which either |S| < k − 2, or equal |S| = k − 2 with the nodes in S forming a connected subgraph of G k . This will allow precalculation and caching of all c(S) needed for computation of right-hand sides.
For efficient precomputation, vertex y 6 ¼ x must meet the following criteria: G\{y} is a connected graph, 3. if d(y) = k − 2, the neighbours of y induce a connected graph, where d(y) represents the degree of y. We will prove that such a vertex exists in any graphlet k ! 4 and all possible x, except for complete graphlets (all vertices violate the first condition) and for the cycle on four points, C 4 (all vertices violate the last condition).
Let L i represent the set of vertices at a distance i from x (see Fig 7). Let l i be the vertex in L i with the smallest degree. Let L u be the last non-empty set, and, accordingly, l u the vertex with the smallest degree among the vertices farthest from x. We will show that l u fulfils the conditions in most cases, except in some for which we can use l u−1 .
Algorithm 2 Selection of node y.
function SELECTY(x) u distance to the furthest node from x L u set of nodes at distance u L u−1 set of nodes at distance u − 1 if l u satisfies the criteria then y lowest degree node from L u else y lowest degree node from L u−1 endif return y end function Each node v 2 L i (i > 0) has at least one neighbour in L i−1 , since the first node in the shortest path from v to x belongs to L i−1 . Consequently, all L i for i u are non-empty. Note also A vertex in L 2 is not adjacent to x by definition of L 2 , and there are no loops, so to have a degree of k − 2 it must be adjacent to all other vertices. Lemma 3 If G is not a complete graph, then d(l u ) k − 2 For u > 1, the lemma follows directly from Lemma 1, so we only need to prove it for u = 1. For contrapositive, assume that d(l 1 ) = k − 1. Since l 1 has the smallest degree in L 1 , all vertices in L 1 have a degree of k − 1. Furthermore, x has a degree of k − 1 since all vertices in L 1 are adjacent to it by definition of L 1 . Hence, G is a complete graph.
The last lemma ensures that the farthest vertex with the lowest degree, l u , fulfills the first condition. It also fulfills the second one: all vertices are connected to x with the shortest paths of lengths at most u, which cannot include l u , thus the removal of l u keeps them connected (at least) via x.
We will prove that l u also fulfills the third condition, except for one special case (d(l u ) = k − 2 and u = 2 and |L 2 | = 1 and k ! 5 and d(l 1 ) k − 2), in which we choose another suitable vertex. We will consider six different cases, which are (except for the trivial first case) illustrated in Fig 7. 1. d(l u ) < k − 2: Condition (iii) does not apply.
2. d(l u ) = k − 2 and u = 1: Since all vertices except x are in L 1 , they are adjacent to x (Fig 7a). x itself is among the neighbours of l 1 , hence neighbours of l 1 are connected through x. doi:10.1371/journal.pone.0171428.g007 3. d(l u ) = k − 2 and u = 2 and |L 2 | > 1: Since l 2 is the vertex with the smallest degree in L 2 , all vertices in L 2 must have a degree of k − 2 and are adjacent to all vertices except x by Lemma 2 (Fig 7b). The neighbour set of l 2 is L 1 [ L 2 \l 2 . Since |L 2 | > 1, there exists a vertex v 2 L 2 s.t. v 6 ¼ l 2 . v is adjacent to all nodes from L 2 [ L 1 , therefore L 1 [ L 2 \l 2 is connected. 4. d(l u ) = k − 2 and u = 2 and |L 2 | = 1 and k = 4: L 1 contains two vertices; both are adjacent to x by definition of L 1 and to l 2 since d(l 2 ) = k − 2 = 2. l 2 is not adjacent to x by definition of L 2 . This leaves only two possible graphs, the cycle C 4 and a diamond (Fig 7c). For the former, the vertex with the required properties does not exist. For the diamond, l u fulfills all three conditions. 5. d(l u ) = k − 2 and u = 2 and |L 2 | = 1 and k ! 5 and d(l 1 ) = k − 1: The neighbour set of l 2 is the entire L 1 (Fig 7d). Since the smallest degree in L 1 is k − 1, L 1 is a complete graph and therefore connected.
6. d(l u ) = k − 2 and u = 2 and |L 2 | = 1 and k ! 5 and d(l 1 ) k − 2: The graph consists of L 0 = {x}, L 2 = {l 2 }, and of L 1 with at least 3 vertices since k ! 5 (Fig 7e). All nodes in L 1 are adjacent to x by definition of L 1 and to l 2 by Lemma 2 since we assume d(l 2 ) = k − 2.
In this case, l u does not always fulfil the conditions, so we choose the lowest degree vertex from L 1 , l 1 . It fulfils the condition (i) by assumptions of this special case. As for condition (ii), the graph G\l 1 is still connected since all points in L 1 are adjacent to x. Since |L 1 | ! 3 and d(l u ) = d(l 2 ) = k − 2, vertices x and l 2 are connected through the remaining vertices in L 1 \l 1 . Condition (iii) needs to be verified just for the case when d(l 1 ) = k − 2 (Fig 7f). The neighbours of l 1 include x, l 2 and all vertices from L 1 except one. Since |L 1 | ! 3, L 1 must include at least one other neighbour of l 1 , which thus connects x and l u .
We have covered all possible cases: the degree of l u cannot exceed k − 2 due to Lemma 3 (assuming the graph is not complete), and when d(l u ) = k − 2, u cannot exceed 2 due to Lemma 1.
We have proven that the vertex with the smallest degree in L u , l u , fulfills the given conditions in all cases except when d(l u ) = k − 2 and u = 2 and |L 2 | = 1 and k ! 5 and d(l 1 ) k − 2. In the latter case, the conditions are fulfilled by l 1 . Complete graphlets and C 4 are handled differently.

Equation for a cycle on 4 nodes
A cycle on 4 nodes, C 4 , is treated separately since there is no suitable node y with the required properties. For C 4 (O 8 ) we choose one of the nodes adjacent to x for the role of y, resulting in x; x 2 2 Nðx 1 Þ; Note that this choice violates the third condition that the neighbours of y should induce a connected graph. The Eq (3) contains a term on the right side that corresponds to the number of common neighbours of node x and some other node x 2 at distance 2 from x. The algorithm stores precomputed values for all such pairs x and x 2 , which would require O(nd 2 ) space and increase the algorithm's space complexity. However, we can still handle this case without consequences for the time and space complexity. We achieve this by reusing O(n) space and recomputing the number of common neighbours every time the algorithm starts a computation of orbit counts for a different node of interest x. This optimization is necessary to keep the space requirement at O(nd) for counting four-node graphlets and does not impact the time complexity.

System of equations
In the constructed system of equations, each orbit is related to orbits from graphlets with higher number of edges. This yields a triangular system of equations: we have one equation for every orbit O and these equations include as terms only the orbit O and other orbits belonging to graphlets with a larger number of edges (e.g., the orbit 59 in Eq (1) is related to orbits 65, 68 and 70).
The system has O k À 1 linear equations for O k orbit counts. To solve it, one orbit count must be enumerated directly. The networks that we encounter in practical applications are usually sparse, which makes the complete graphlet (clique) a good candidate. Because of its rarity and symmetricity, we can efficiently restrict the enumeration.
Enumerating the orbit in the graph with the largest number of edges also simplifies solving the given triangular system of equations.

Extension to edge orbits
Edge orbits (Fig 8) can be defined in a similar way as node orbits. We can use the same approach for setting up the corresponding system of equations. Since the system does not refer to a single x but to an edge (x a , x b ), the selected node y must not coincide with either of these endpoints. We can set x = x a and show that we can always choose a node y 6 ¼ x b with the required properties.
Since x b is always in L 1 , we have to analyze only the cases where we choose y from L 1 . This happens in cases 1, 2, and 6 from the proof in section about additional constraints on y. We need to check that there at least two suitable vertices for y, so if one of them is x b , the other is chosen as y.
1. d(l u ) < k − 2: We need to consider only the case when u = 1. Since d(l u ) < k − 2, all vertices in L 1 satisfy condition (i). The remaining graph is connected through x (condition (ii)), and condition (iii) does not apply.
2. d(l u ) = k − 2 and u = 1: Recall that the graph is not complete. Since all vertices in L 1 are connected to x, there must be at least one pair in L 1 that is not connected and thus has a degree of k − 2. These two vertices satisfy condition (i). Conditions (ii) and (iii) again hold since all vertices are connected through x.
6. d(l u ) = k − 2 and u = 2 and |L 2 | = 1 and k ! 5 and d(l 1 ) k − 2: Since the vertex in L 2 has a degree of k − 2, it is connected to all vertices in L 1 ; nodes in L 1 are connected to x. The nodes in L 1 do not induce a complete graph (d(l 1 ) k − 2), so there must again exist a disconnected pair in L 1 , which satisfies all conditions like in above case 2.
For C 4 , one of the nodes in L 1 is x b and the other is y.

Algorithm
Coefficients on the left-hand side of the relations are related to symmetry properties of the graphlets and not to the graph H. The terms on the right sides of equations depend on the host graph H. Their computation requires enumeration of all graphlets of size k − 1 and adding up their possible extensions. The first step is pre-computation and storing of cðSÞ for all subsets S with up to k − 3 vertices and for all connected subsets of k − 2 vertices. These conditions match the criteria for selection of y, so the precomputed values cðSÞ represent the terms in the sum on the right-hand sides of equations. This is followed by direct enumeration of cliques with k vertices. This enumeration does not have to be extremely fast, but just fast enough not to dominate the time complexity of the entire graphlet counting algorithm. For this purpose we can employ some of the approaches to listing cliques [20,21].
Following this precomputation, the next two steps are repeated for each vertex x 2 V.

Computation of sums on the right-hand sides of equations.
Computation is implemented as enumeration of (k − 1)-node graphlets touching x, as specified by the conditions under the sums. For each graphlet, the terms in the sum consist of the counts cðSÞ precomputed in the previous step. For k 5, the number of graphlets with k − 1 nodes is small, so it is feasible to design efficient individual procedures for enumerating them. These procedures involve early pruning of non-viable candidates and completely avoiding any isomorphism testing. Medium-sized graphlets (k = 5 or 6) require graphlet recognition of enumerated connected subgraphs, however these patterns can be efficiently distinguished with the use of some trivial invariants such as a degree sequence. Enumeration of larger graphlets would benefit from efficient methods for isomorphism testing.
Solving the system of equations. The system is triangular, with each equation relating one orbit to those with larger number of edges, Since the count for the highest orbit, which belongs to the clique, is computed by direct enumeration, the system can be solved by decreasing orbit indices. All orbit counts, coefficients and free terms are integers, thus the computation is numerically stable.

Time-and space-complexity
We will analyze the worst-case complexity and the expected complexity on random Erdős-Rényi graphs, followed by empirical verification. Worst-case complexity. We will evaluate the worst-case time complexity of the algorithm in terms of the number of nodes (n) and the maximum degree of a node (d) in the host graph. We treat the size of the graphlets, k, as a constant. We assume that the graph is stored as a list of adjacent nodes together with a hash table for checking whether two nodes are connected in constant time. The algorithm consists of four steps.
Precomputation of common neighbours. We need to precompute the number of common neighbours of sets of k − 3 or fewer nodes and of connected sets of k − 2 nodes to efficiently construct right sides of our equations. To achieve this we enumerate all subsets of k − 2 or fewer neighbours for every node. This results in time complexity O(nd k−2 ). Storing the number of common neighbours of sets of at most k − 3 nodes with the above method requires O(nd k−3 ) space. Because we request that in the case of k − 2 nodes, these nodes induce a connected subgraph, we can limit their number to the number of (k − 2)-node induced connected subgraphs, which is also O(nd k−3 ).

Enumeration of cliques.
We will refer to the time complexity of counting k-node cliques in this step as O(T k ). A worst-case time complexity is O(nd k−1 ) and requires constant space. However, this enumeration can be implemented very efficiently in practical applications on sparse networks that contain few cliques.
Enumerating all (k − 1)-node graphlets and counting their extensions. This step computes the right sides of the system of equations. It requires constant space, since the space is reused for each vertex, and runs in O(nd k−2 ) time needed for enumeration of k − 1-node graphlets.
Solving the system of equations. The system of equations is independent of the host graph and requires constant time and space.
Overall, the algorithm has a O(nd k−2 + T k ) time complexity while requiring O(nd k−3 ) space. In the worst case, the time complexity is the same as that of a simple exhaustive enumeration method, O(T k ) = O(nd k−1 ). However, the term T k is much smaller in practice.
Expected time complexity in random graphs. Although the worst-case time complexity of the algorithm is equal to that of brute-force enumeration, the actual performance on realworld networks and on random graphs is much better. We analyzed the expected time complexity on random Erdős-Rényi graphs with n nodes and edge probability p. Throughout this analysis we will assume that np > 1, otherwise the graph is likely to have more than one component which can be processed independently.
The precomputation consists of iterating over central nodes, enumerating all sets of l k − 2 neighbours and incrementing the number of common neighbours of the leaf nodes. The l nodes have to be connected to the central node, which happens with probability p l . The expected time complexity of this step is Oðn P kÀ 2 l¼1 n l p l Þ. Assuming np > 1, we can simplify it to O(n k−1 p k−2 ).
In the second step, the algorithm enumerates all subgraphs with k − 1 nodes. It does so incrementally by first enumerating smaller connected subgraphs of size l and extending them to larger connected subgraphs. The expected time complexity is therefore proportional to the expected number of connected subgraphs with l k − 1 nodes. We need to estimate the probability that a set of l nodes induces a connected subgraph. We can view the process of building every such subgraph by consecutively attaching a new node to at least one of the existing nodes. This of course will overestimate the number of connected subgraphs by some constant because every such subgraph can be built in several different orders of attaching nodes. The probability that an edge exists from some newly added node to at least one of the i existing nodes is 1 − (1 − p) i . The expected number of enumerated subgraphs is therefore Oð P kÀ 1 l¼1 n l Q lÀ 1 i¼1 ð1 À ð1 À pÞ i ÞÞ ¼ Oð P kÀ 1 l¼1 n l p lÀ 1 Þ. Assuming np > 1, the expected time complexity is O(n k−1 p k−2 ).
The total expected time complexity for setting-up the system of equations in Erdős-Rényi graphs with n nodes and edge probability p is thus O(n k−1 p k−2 ). In practice, we observe graphlets with 4 and 5 nodes. The expected time complexities for these cases are O(n 3 p 2 ) and O(n 4 p 3 ), respectively.

Empirical evaluation of time complexity.
We evaluated the performance of our algorithm for counting 4-and 5-node graphlets on random Erdős-Rényi graphs.
We measured the time needed for counting node-and edge-orbits with the Orca algorithm and compared it to a bruteforce enumeration. For the latter we used an implementation from GraphCrunch. Orca outperforms exhaustive enumeration by an order of magnitude (Tables 2  and 3). The running times for counting node-orbits and edge-orbits are practically the same in the case of counting 4-or 5-node graphlets in random graphs with 1 000 nodes and of increasing density. The size of the graphs (n = 1000) was chosen arbitrarily to put the run times in the range of a couple of seconds. In the remainder of this section we focus on counting nodeorbits.
Second, we compare the running times of counting orbits of 4-node and 5-node graphlets on random graphs with 10 000 nodes and up to 800 000 edges. These graphs are sparse as the number of edges represents only about 1.6% of all possible edges; as such they represent a realistic case of large network analysis. A logarithmic plot of execution times in Fig 9 shows a polynomial dependence on the size of graphlets. Dotted lines correspond to a bruteforce enumeration method and solid lines to our Orca algorithm. The line corresponding to bruteforce counting of 4-node graphlets aligns nicely with counting 5-node graphlets using Orca. This reflects the enumeration of 4-node graphlets used to compute 5-node graphlet counts. The only difference is by a constant factor.
The slope of the Orca's lines should be 2 and 3, respectively, according to the expected time complexities (O(n 3 p 2 ), O(n 4 p 3 )). However, this is clearly not the case in Fig 9. Further experiments show that this is the result of CPU cache misses when accessing the precomputed lookup tables. We performed a similar experiment with disabled CPU cache. Because of the slowdown, we decreased the number of nodes to 1 000 and maintained average degree of nodes. The measurements with disabled CPU cache in Fig 10 line up with the expected slopes of 2 and 3.
Finally, we probed for the region in which the enumeration of cliques begins dominating the time complexity. We performed the experiment for counting 4-node graphlets in graphs with 1 000 nodes and increasing edge probabilities p. In Fig 11 the plot follows a straight line up to around p = 0.07 and another steeper line from p = 0.3 onwards. This is consistent with the contribution of the step of enumerating cliques. Random sparse graphs contain fewer cliques whose enumeration is efficient and does not significantly affect the running time. However, as the graphs become denser, this becomes the bottleneck of the algorithm.

Final remarks
The source code of the algorithm in C++, which computes the node and edge orbits for k = 4 and k = 5, along with the randomly generated data used in experiments is available at https:// github.com/thocevar/orca. The corresponding R package orca is also available on CRAN.
Parts of this algorithm that have been presented previously are also already included in the GraphCrunch package [22].

Author Contributions
Conceptualization: TH JD.
Formal analysis: TH JD.