Extracting the hierarchical organization of complex systems

Background

Complex networks are convenient representations of the interactions within complex systems . Here, we focus on the identification of inclusion hierarchies in complex neworks, that is, to the unraveling of the nested organization of the nodes in a network into modules, which are comprised of sub-modules and so onWe do not consider other hierarchical schemes that classify nodes according to, for instance, their importance . Another issue that we do not address here is that of “overlapping” modules. In the literature, some authors refer to the existence of “soft” boundaries between communities . However, there has been so far no rigorous connection between the soft boundaries and the overlap between communities. Moreover, at present, there is no theoretical model that includes overlapping modules, that is, modules that share nodes, as opposed to communities that share edges..

A method for the identification of the hierarchical organization of nodes in a network must fulfill two requirements: (i) it must be accurate for many types of networks, and (ii) it must identify the different levels in the hierarchy as well as the number of modules and their composition at each level. The first condition may appear as trivial, but we make it explicit to exclude algorithms that only work for a particular network or family of networks, but that will otherwise fail. The second condition is more restrictive, as it excludes methods whose output is subject to interpretation. Specifically, a method does not fulfill the second condition if it organizes nodes into a tree structure, but it is up to the researcher to find a “sensible” criterion to establish which are the different levels in that tree. An implication of the previous two requirements is that any method for the identification of node organization must have a null output for networks, such as Erdős-Rényi random graphs, which do not have an internal structure.

To our knowledge, there is no procedure that enables one to simultaneously assess whether a network is organized in a hierarchical fashion and to identify the different levels in the hierarchy in an unsupervised way. Ravasz et al. studied the hierarchical structure of metabolic networks, but in their analysis the authors put emphasis on detecting “global signatures” of a hierarchical network architecture. Specifically, they reported that, for the metabolic networks studied and for certain hierarchical network models, the clustering coefficient of nodes appears to scale with the connectivity as $C(k)\sim k^{-1}$. This scaling, however, is neither a necessary nor a sufficient condition for a network to be hierarchical .

More direct methods to investigate the hierarchical organization of the nodes in a network have also been recently proposed . Although useful in some contexts, these methods do not clearly identify hierarchical levels and thus fail to satisfy condition (ii) above. Furthermore, all these methods yield a tree even for networks with no internal structure.

In the following, we define inclusion hierarchies in complex networks and describe an ensemble of hierarchically nested random graphs. We then introduce a method that is able to accurately extract the hierarchical organization of hierarchical random graphs. Finally, we apply our method to several real-world networks.

Inclusion hierarchies

Let $\mathcal{G}_{i}\subset\mathcal{G}$ be the set of groups in which node $n_{i}$ holds membership. Here, we consider that node $n_{i}$ holds membership in only one group per level, and that membership to groups follows a nested hierarchy. Therefore, for node $n_{i}$ to hold membership in group $g_{11}$, node $n_{i}$ must also hold membership in group $g_{1}$.

We assume that the probability $p_{ij}$ of the edge $(n_{i},n_{j})$ being present in a network is a function solely of the set of co-memberships $\mathcal{M}_{ij}=\mathcal{G}_{i}\cap\mathcal{G}_{j}$ of the two nodes. Note that our assumptions imply that: (i) $\mathcal{M}_{ij}$ obeys transitivity, so that if $\mathcal{M}_{ij}=\mathcal{M}_{ik}$, then $\mathcal{M}_{ij}=\mathcal{M}_{jk}$; and (ii) node memberships in groups {$\,g_{k_{1}k_{2}}\,$} at the second level are uniquely and completely defined by the sub-network of connections of all nodes holding membership in group $g_{k_{1}}$, that is, information at deeper levels in the hierarchy is totally decoupled from the information at higher levels in the hierarchy.

In the simplest scenario, $p_{ij}$ is a non-decreasing function of the cardinality $x$ of $\mathcal{M}_{ij}$, which implies that groups of nodes holding membership in the same groups will be more densely connected than a randomly selected group of nodes. This is precisely the underlying assumption in many algorithms aiming to detect the top level community structure of complex networks assuming a flat organization of the nodes .

Let us now introduce an ensemble of random networks which are constructed following hierarchical node membership assignment: hierarchically nested random graphs. We restrict our ensemble to networks with a homogeneous hierarchical organization of the nodes (see Supplementary Information for other kinds of hierarchical organization) that have the same degree distribution as Erdős-Rényi graphs .

Extracting the hierarchical organization of networks

Our method consists of two major steps (Fig. 1): (i) measuring the “proximity” in the hierarchy between all pairs of nodes, which we call node affinity; and (ii) uncovering the overall hierarchical organization of node affinities, or, in other words, detecting the underlying organization of group memberships.

A standard approach for quantifying the affinity between a pair of nodes in a network is to measure their “topological overlap” , which is defined as the ratio between the number of common neighbors of the two nodes and the minimum degree of the two nodes. This measure identifies affinity between nodes with a dense pattern of local connections. Because topological overlap is a local measure, it will fail to detect any structure when a network is not locally dense (Fig. 2).

We propose a new affinity measure based on surveying of the modularity landscape , a collective property of the network. Our definition of affinity between nodes draws upon the idea that modules correspond to sets of nodes which are more strongly interconnected than one would expect from chance alone . We show below that our affinity measure detects the modular structure even in the absence of a dense pattern of local connections.

Consider the ensemble $\cal P$ of all partitions of a network into modules , and assign to each partition $P$ the modularity

where $L$ is the total number of links in the network, $l_{i}$ is the number of links within module $i$, $d_{i}$ is the sum of degrees of all the nodes inside module $i$, and the sum is over all the $m$ modules in partition $P$ (Fig. 1A). The modularity of a partition is high when the number of intra-module links is much larger than expected for a random partition.

Let $\cal P_{\rm max}$ be the set of partitions for which the modularity $M$ is a local maxima, that is, partitions for which neither the change of a single node from one module to another nor the merging or splitting of modules will yield a higher modularity . Let $B_{\rm max}=\{\,b(P):P\in\cal P_{\rm max}\,\}$ be the sizes of the “basin of attraction” of those maxima. The affinity $A_{ij}$ of a pair of nodes $(i,j)$ is then the probability that when local maxima $P\in\cal P_{\rm max}$ are sampled with probabilities $b(P)$, nodes $(i,j)$ are classified in the same module.

Note that, in contrast to other affinity measures proposed in Refs. , the measure we propose does not necessarily coincide with the “optimal” division of nodes into modules, that is, the partition that maximizes $M$ . In fact, the modules at the top level of the hierarchy do not necessarily correspond to the best partition found for the global network, even for relatively simple networks (Fig. 2C).

Statistical significance of hierarchical organization—

Given a set of elements and a matrix of affinities between them, a commonly used tool to cluster the elements and, presumably, uncover their hierarchical organization is hierarchical clustering . Hierarchical clustering methods have three major drawbacks: (i) They are only accurate at a local level—at every step a pair of units merge and some details of the affinity matrix are averaged with an inevitable loss of information; (ii) the output is always a hierarchical tree (or dendogram), regardless of whether the system is indeed hierarchically organized or not; (iii) there is no statistically sound general criterion to determine the relevant levels on the hierarchy.

In order to overcome the first caveat of agglomerative methods such as hierarchical clustering, one necessarily has to follow a top to bottom approach that keeps the details of the matrix. That is the spirit of divisive methods such as k-means or principal component analysis , which group nodes into “clusters” given an affinity matrix. However, these methods have a significant limitation: the number of clusters is an external parameter, and, again, there is no sound and general criterion to objectively determine the correct number of clusters.

Because of the caveats of current agglomerative and divisive methods, we propose a “box-clustering” method that iteratively identifies in an unsupervised manner the modules at each level in the hierarchy. Starting from the top level, each iteration corresponds to a different hierarchical level (Fig. 2).

In order to assess whether the network under analysis has an internal organization we need to compare with the appropriate null model, which in this case is an ensemble of “equivalent” networks with no internal organization. These equivalent networks must have the same number of nodes and an identical degree sequence. A standard method for generating such networks is to use the Markov-chain switching algorithm . Despite their having no internal structure, these networks have numerous partitions with non-zero modularity . Thus, to quantify the level of organization of a network, one needs to compare the modularities of the sampled maxima for the original network and its corresponding random ensemble; if the network has a non-random internal structure, then local maxima in the original landscape should have larger modularities than local maxima in the landscapes of the randomized networks.

Specifically, for a given network, we compute the average modularity $M_{\rm av}$ from { $M(P)\>:\>P\in\cal P_{\rm max}\>$}. Then, we compute the same quantity $M^{i}_{\rm av}$ for each network in the equivalent random ensemble. In virtue of the central limit theorem, the set of average modularities for the whole ensemble {$\>M^{i}_{\rm av}\>$} is normally distributed with mean $M_{\rm rand}$ and variance $\sigma^{2}_{M_{\rm rand}}$. To quantify the level of organization of a network, we thus compute the z-score of the average modularity

If $z$ is larger than a threshold value $z_{t}$, then the network has internal structure and we proceed to identify the different modules, otherwise we conclude that the network has no structure. In what follows, we show results for $z_{t}=2.3267$, which corresponds to a 1% significance level (Supplementary Material)Results for real networks at a 5% significance level are identical, however, the more stringent threshold is more efficient at detecting the last level in the hierarchy for model networks. Only for a 1-3% of the cases—depending on the cohesiveness of the levels—do we find that algorithm finds one more level than expected..

Building the hierarchical tree—

In networks organized in a hierarchical fashion, nodes that belong to the same module at the bottom level of the hierarchy have greater affinity than nodes that are together at a higher level in the hierarchy. Thus, if a network has a hierarchical organization, one will be able to order the nodes in such a way that groups of nodes with large affinity are close to each oder. With such an ordering, the affinity matrix will then have a “nested” block-diagonal structure (Fig. 1). This is indeed what we find for networks belonging to the ensemble of hierarchically nested random graphs (Fig. 2).

For real-world networks, we do not know a priori which nodes are going to be co-classified together, that is, we do not know which is the ordering of the nodes for which the affinity matrix has a nested block-diagonal structure. To find such an ordering, we use simulated annealing to minimize a cost function that weighs each matrix element with its distance to the diagonal

where $N$ is the order of the affinity matrix (see Fig. 1A and Supplementary Information for alternative ordering schemes).

This problem belongs to the general class of quadratic assignment problems . Other particular cases of quadratic assignment problems have been suggested to uncover different features of similarity matrices Our algorithm is able to find the proper ordering for the affinity matrix and to accurately reveal the structure of hierarchically nested random graphs (Fig. 2).

Importantly, we want to balance the ability of the model to accurately describe the data with its parsimony, that is, we do not want to over-fit the data. Thus, we use the Bayesian information criterion in order to determine the best set of boxes We have also applied Akaike’s information criterion , obtaining the same results for most of the cases..

To find the modular organization of the nodes at the top level (level 1), we fit the block diagonal model to the global affinity matrix. As we said previously, we assume that the information at different levels in the hierarchy is decoupled, thus in order to detect sub-modules beyond the first level, one needs to break the network into the sub-networks defined by each module and apply the same procedure (Fig. 1). The algorithm iterates these steps for each identified box until no sub-networks are found to have internal structure.

Method validation

We validate our method on hierarchically nested random graphs with one, two, and three hierarchical levels. We define the accuracy of the method as the mutual information between the empirical partition and the theoretical one . Figure 2C shows that the algorithm uncovers the correct number of levels in the hierarchy.

Moreover, our method always detects the top level, even for the networks with three hierarchical levels. In contrast, because the partition that globally maximizes $M$ corresponds to the sub-modules in the second level, even the more accurate module identification algorithms based on modularity maximization would fail to capture the top level organization (Joshi et al. 2007, ).

The hierarchically nested random graphs considered above have a homogeneous hierarchical structure; however, real-world networks are not likely to be so regular. In particular, for real-world networks one expects that some modules will have deeper hierarchical structures than others. We thus have verified that our method is also able to correctly uncover the organization of model networks with heterogeneous hierarchical structures (Supplementary Information).

Testing on real world networks

Having validated our method, we next analyze different types of real-world networks for which we have some insight into the network structure: the world-wide air-transportation network , an e-mail exchange network of a Catalan university , and an electronic circuit .

In the air transportation network, nodes correspond to airports and two nodes are connected if there is a non-stop flight connecting them. In the email network, nodes are people and two people are connected if they send emails to each other. In the electronic network, nodes are transistors and two transistors are connected if the output of one transistor is the input of the other (Table 1).

We find that the air-transportation network is strongly modular and has a deep hierarchical organization (Fig. 3). This finding does not come as a surprise since historical, economic, political, and geographical constraints shape the topology of the network . We find eight main modules that closely match major continents and sub-continenets, and major political divisions and thus truly represent the highest level of the hierarchyThe ability of the present method to detect the top level is significant. A previous study co-authored by two of us identified 19 modules in the world-wide air-transportation network using the most accurate module detection algorithm in the literature ..

The electronic circuit network is comprised of eight D-flipflops and 58 logic gates . Our method identifies two levels in the network (Fig. 4A). At the top level, modules are groups of logic gates, all the logic gates comprising a D-flipflop being in the same module. At the second level, the majority of modules comprise single gates.

For the email network, five of the seven major modules at the top level (Fig. 4B) correspond to schools in the university, with more than 70% of the nodes in each of those modules affiliated to the corresponding school. The remaining two major modules at the top level are a mixture of schools and administration offices (often collocated on campus), which are distinctly separated at the second level. The second level also identifies major departments and groups within a school, as well as research centers closely related to a school.

Application to metabolic networks

Finally, we analyze the metabolic networks of E. coli obtained from two different sourcesIn the Supplementary Material we also show the organization obtained for the metabolic network for E. coli from the Ma-Zeng database , and for the metabolic network of H. pylori developed at UCSD . (Fig. 5): the KEGG database , and the reconstruction compiled by Palsson’s Systems Biology Lab at UCSD . In these networks, nodes are metabolites and two metabolites are connected if there is a reaction that transforms one into the other .

To quantify the plausability of our classification scheme, we analyze the within-module consistency of metabolite pathway classification for the top and the second levels of the metabolic network for E. coli reconstructed at UCSD . For each module, we first identify the pathways represented; then, we compute the fraction of metabolites that are classified in the most abundant pathway. We find that there is a clear correlation between modules and known pathways: At the top level, for all the modules except one, we find that the most abundant pathway comprises more than 50% of the metabolites in the module.

For the second level, we find that for most of the modules all the metabolites are classified in the same pathway. We also detect smaller pathways that are not visible at the top level (such as those for polyketides and nonribosomal peptides, and for secondary metabolites).

Our results thus provide an objective description of cellular metabolism that, while not affected by human subjectivity, captures our current understanding of these networks. Interestingly, “known” pathways do not correspond to a single module at the top level, implying that large pathways are in fact comprised of smaller units. Intriguingly, these units are not necessarily uniform in “pathway composition” but are a mixture of sub-modules associated to different pathways. Thus, an important question is how the modules we identify relate to metabolism evolution .

References