Identifying modular flows on multilayer networks reveals highly overlapping organization in social systems

Unveiling the community structure of networks is a powerful methodology to comprehend interconnected systems across the social and natural sciences. To identify different types of functional modules in interaction data aggregated in a single network layer, researchers have developed many powerful methods. For example, flow-based methods have proven useful for identifying modular dynamics in weighted and directed networks that capture constraints on flow in the systems they represent. However, many networked systems consist of agents or components that exhibit multiple layers of interactions. Inevitably, representing this intricate network of networks as a single aggregated network leads to information loss and may obscure the actual organization. Here we propose a method based on compression of network flows that can identify modular flows in non-aggregated multilayer networks. Our numerical experiments on synthetic networks show that the method can accurately identify modules that cannot be identified in aggregated networks or by analyzing the layers separately. We capitalize on our findings and reveal the community structure of two multilayer collaboration networks: scientists affiliated to the Pierre Auger Observatory and scientists publishing works on networks on the arXiv. Compared to conventional aggregated methods, the multilayer method reveals smaller modules with more overlap that better capture the actual organization.


Introduction
Multifaceted relationships between numerous components in social and biological systems make them inherently complex to analyze 1,2 . Data about these interactions have become increasingly available and network analysis have emerged as an essential tool for studying their function [3][4][5] . For large networks, detailed modeling of individual components and their interactions is unfeasible and researchers instead seek to simplify and highlight important large-scale functional structures in the networks. Depending on the system under study and the research question at hand, researchers use methods that either operate on the plain topology of the network itself 6,7 or, to capture flow processes through the real system, on dynamics modeled on the network 8,9 . In any case, an important objective is to detect so-called communities 10 , topological groups of nodes with higher internal than external density of links compared to null models [11][12][13] or, alternatively, modules that capture flows for a relatively long time [14][15][16] .
However, community-detection methods generally assume that a single type of static link, weighted and directed at best, can account for all types of interactions between nodes in the network. This assumption oversimplifies the multifaceted nature of relationships in real systems with important consequences. Aggregating multiple types of relationships into a single weighted and directed network can distort both the topology of the network and the dynamics on the network 17 . Take social relationships as an example, where the way the same individual interacts with her relatives, friends, and colleagues may depend on location, time, or means of interaction. Is she at home or at the office? Is it weekday or weekend? Is she communicating by phone or by Facebook?
If all contact events are aggregated into a single network layer, important temporal 18 and structural 19 information is inevitably lost. Recently it has been shown that multilayer networks provide an effective framework to capture different types of interactions between nodes 20-26 , including a generalization of modularity to identify groups in multilayer networks 20 . While the generalized null models of modularity are based on Laplacian dynamics 20 , they nevertheless favor topological groups with high link density 27 , both within and between network layers 28 .

Modules in multilayer networks: a flowapproach
To identify modular flows on multilayer networks, here we introduce a method based on compression of network flows. The information-theoretic method generalizes the so-called map equation 15 for networks with memory 17 to take advantage of modular flows in multilayer networks. The framework generalizes straightforwardly, because the information-theoretic machinery remains the same and only the flow model changes, with memory of present layer rather than of previous step. This approach therefore suggests a natural concept of communities in multilayer networks as groups of nodes that capture flows within and across layers for a relatively long time.
We begin by describing how we model the dynamics and then introduce the multiplex map equation. We measure the performance on benchmark networks and contrast with results obtained with the generalization of modularity. Finally we analyze the modular flow dynamics on two multilayer collaboration networks. We have integrated the method in the Infomap software package available online 29 .

Flow dynamics on multilayer networks
A multilayer network is an efficient representation of a connected system of agents that may interact in different roles, at different times, or by different means, for example. For simplicity, we refer to them as different modes. Each physical node represents an agent, and each network layer represents the constraints on flow among the agents in a given mode. Figure 1 illustrates a multilayer network with three layers and four physical nodes. We use Latin letters to enumerate the physical nodes, Greek letters to enumerate the network layers, and pairs of Latin and Greek letters to identify node-layer tuples 30 , corresponding to physical nodes in specific network layers, which we in the following refer to as state nodes (Fig. 1B). Sometimes empirical data allow us to assign weights to both intra-and inter-layer links between state nodes. In such interconnected networks, we have complete information to model dynamics with a random walker that follows links proportional to their weights within and between network layers. In general, with intra-layer adjacency matrix W β ij of layer β and inter-layer adjacency matrix D αβ i of physical node i, the transition probabilities are where S α i = β D αβ i are the inter-layer out-strengths and s β i = j W β ij are the intra-layer out-strengths of node i in layer α and β, respectively 26 (see Methods). In practice, however, often data about inter-layer link weights are scarce. That is, information about the probability of switching layer is incomplete.
In absence of empirical inter-layer weights, we use random walker dynamics with relax rate r to model movements between layers. In a given step, with probability 1 − r, the random walker moves according to the intralayer links of the state node, and with probability r, the constraint to move in the current layer is relaxed and the random walker moves along any link of the physical node. In this way, the random walker switches from layer α to layer β with probability s β i /S α i . These dynamics are described by the transition probabilities with S i = β s β i independent of layer. 1 A relaxed step on a multilayer network resembles a teleportation step in the PageRank algorithm 8 , which allows a random surfer to move freely to a random website and explore 1 It is worth noting that Eq. (2) is equivalent to Eq. (1) when the full network. However, a relaxed step only frees the constraints set by the current network layer and allows the random walker to follow a link from node i to node j in any network layer (Fig. 1B). Accordingly, changing the relax rate from 0 to 1 modifies the constraints on the random walker from those that force it to be stuck in disconnected network layers to those that allow it to move more freely on the fully aggregated network. In this way, we can model the important interplay between interconnected network layers.

Communities in multilayer networks
There are, in principle, many ways to define communities in multilayer networks 20,31 . The challenge is to construct an effective framework. The challenge may seem daunting, since it is still debated how to define communities in single-layer networks 10 , and multilayer networks are inherently more complex with simultaneous and nonlinear coupling between the layers. However, by using the fact that many networks represent constraints on flow in social and biological systems, and that multilayer networks are just a more complete description of those constraints, a generalization of flow-based community detection methods follows straightforwardly.
We begin by exemplifying how we identify communities in a multilayer network. As an illustrative example, we use a social system in which nodes represent individuals and network layers represent family, friendship, and work relations, respectively. The constraints on flow in a network layer may give rise to modules with long flow persistence times. Moreover, and importantly, the modules in each network layer may or may not depend on other network layers. For example, if some friends run a business together, their module in the friendship-relations layer will correlate with their module in the work-relations layer, such that they form a single reinforced module across the two layers. Contrarily, all members of a family may not work together or even hang out as friends, such that the family module does not extend across layers. However, if some of the family members run a business together or hang out as friends, modules may overlap. That is, identifying modular flows on multilayer networks captures that individuals can belong to multiple highly interactive communities with limited information transfer between, such that information has long persistence times within communities. The schematic multilayer network in Fig. 1 illustrates. Each layer has a triangle of connected nodes that trap flow for a long time and form a module. The red network layer has very little overlap with the two identical blue network layers (Fig. 1A). By first representing the multilayer network as a multiplex network with state nodes (Fig. 1B), and then releasing a random walker on the multiplex network, the community structure with two overlapping modules appear (Fig. 1C). In next the section, we make this concept of modular flow in multilayer and interconnected networks precise by gen-eralizing the map equation. The three layers represented as a multiplex network with physical nodes in black and state nodes i, α in red, blue, and dashed blue connected with intra-layer link weights W α ij . (C) A random walker on the multiplex network moving between the state nodes, twice relaxing the layer constraint and following any link from the physical node of the currently visited state node (first time in j, γ and second time in i, β). While the random walker moves according to the weights between the state nodes, only the physical nodes are considered to be observables, as illustrated by the sequence of physical nodes that the random walker has visited. When the random walker moves along links of the red layer, it is trapped in the lower right triangle. When the random walker moves along links of the blue or dashed blue layer, it is trapped in the upper left triangle. As a consequence, the multilayer network has two overlapping modules with respect to flow.

The multiplex map equation
The map equation measures the length required to communicate dynamics on a network with a modular description 15 . Unlike the maximum compression achieved by the entropy rate of a random walk process on the network 32,33 , the coding scheme of the map equation grants unique names to important structures of the network. In this modular coding scheme, each entry to a module, and each node visit and module exit of a module, is assigned a unique code word. With these constraints on the coding scheme, maximum compression is achieved when nodes that capture the random walk process for a relatively long time are assigned to the same module. As a result of the duality between compressing data and finding regularities in the data, identifying the assignment of nodes into modules with maximum compression simultaneously answers: How many modules are present? And which nodes are members of which modules?
The original map equation can be generalized by modifying the constraints on the coding scheme. Explicitly, the three constraints on the original two-level version of the map equation for hard partitions are: (i ) A modular code structure with unique names on important structures, (ii ) movements in no more than two levels, modules and nodes, and (iii ) that each node can only belong to one module. Constraint (i ) is essential and cannot be relaxed, but relaxing constraint (ii ) allows for multilevel solutions 34 and relaxing constraint (iii ) allows for overlapping-module solutions 35 . But how should the map equation be generalized for multilayer structures?
The natural generalization is to capture the modified dynamics of the multilayer network while maintaining the essential constraint (i ) of a modular code structure with unique names on important structures. This principle is particularly well suited to capture the fundamental notion of multilayer networks, namely that it is the very same physical object that is represented by its state nodes in each layer. Therefore, the generalization follows straightforwardly. Whenever two state nodes that represent the same physical object are assigned to the same module, they should use a common code word (see Fig. 1C and Methods). In colloquial terms for the example above, the colleagues who are also friends will refer to each individual by a single name that may be different from what family members use. In this way, the multilayer network modules will naturally overlap if the dynamics have such properties. This generalization can be taken one step further by relaxing constraint (ii ) and allow for multilevel solutions with nested modules. In fact, we have integrated both the two-level and the multilevel multiplex map equation in the Infomap search algorithm available online 29 , but here we focus on two-level modular structures, communities.

Results and Discussion
In this section, we first validate our framework on novel multilayer benchmark networks and then analyze two inherently multilayer collaboration networks.

Performance tests on multilayer benchmark networks
To test the performance of the information-theoretic and flow-based method, we developed multilayer benchmark networks with modular structure across layers. We followed the standard approach and obtained benchmark networks from a generative model in which nodes are assigned to communities and the probability of drawing a link between two nodes depends on their community assignments 36,37 . While the multiplex map equation can identify modules that independently span across any number of layers, here we consider benchmark networks with community structure in entire layers that either correlate or not. This more easily tractable scenario nevertheless highlights salient features of modular flows. As schematically illustrated in Fig. 1, the scenario corresponds to systems that can be in different modes with dependent network layers. Using the example from above, colleagues would also be friends such that the two layers would have almost the same community structure, yet different from the community structure associated with family relations. Such redundant information is common in many social and biological networks that represent systems that can be in different modes as a whole or slowly change over time 38 .
Specifically, we first generated T independent LFR benchmark networks 37 for the different modes of the system, and then sampled L network layers from each of the mode networks. We sampled the network layers by including each link with probability 1/L to allow a link of the mode network to be sampled once on average. Each multilayer benchmark network thus comprises T × L layers, with T sets of L dependent layers. Figure 2 schematically illustrates a multilayer benchmark network with T = L = 2. The challenge is to reveal the community structure of each mode network, which corresponds to simultaneously reveal the community structure in each layer and identify the mode network from which the layer was sampled. To make the test realistic, we only provided the algorithm with the T × L network layers, and neither input any information about the number of mode networks T , nor about how or in which order the layers were sampled. In the small example illustrated in Fig. 2, generalized modularity 20 correctly identifies the communities in each layer but fails to identify the communities in the two original mode networks. Contrarily, the multiplex map equation, here with relax rate r = 0.15, both identifies the communities in each layer and the communities in the mode networks. We generate the multilayer networks in two steps. (A) First we generate T LFR benchmark networks with well-defined communities, here illustrated with two network modes in blue and red. (B) Then we sample L network layers from each mode network, here illustrated with four network layers in total. (C) Each state node in the multilayer network is classified in a community, such that communities of physical nodes map overlap. In partition 1, each state node is correctly classified. In partition 2, however, the light and dark color shades are assigned to the same module, respectively. While these communities provide the correct partition of each slice, they fail to capture the communities of the two original mode networks. Generalized modularity favors partition 2 whereas the multiplex map equation favors partition 1.
The multiplex map equation can accurately identify multilayer communities. To test the performance more systematically, we generated multilayer benchmark networks with different number of mode networks T and network layers per mode network L. For the mode networks, we used LFR benchmark networks with 128 nodes and 4 communities, each with 32 nodes with average degree 16, and the fraction of inter-community links set to 0.05. We varied T between 1 and 3, and L between 1 and 7. To quantify the performance, we applied the normalized mutual information (NMI) 39,40 to state nodes. In this way, we can quantify how well the method captures the multilayer communities. Figures 3A and B show the results for relax rate r = 0.15. Optimization of the multiplex map equation with Infomap, Multiplex Infomap for short, accurately identifies the communities of the original mode networks for up to 5-6 network layers per mode network. Contrarily, standard Infomap applied on each layer separately or on the supra-adjacency representation of the multilayer network with all state nodes interpreted as physical nodes 24 only succeed for one layer per mode network. That is, only by acknowledging the multiplex nature of the benchmark networks is it possible to accurately identify their multilayer communities. The results are only weakly dependent on the relax rate (see Fig. 3C), although the exact range depends on the relative constraints on flow manifested in network layers of the same and different mode networks (see Figs. 1 and 2 in the SI). When nothing else is stated, we use r = 0.15 throughout our analysis. With this relax rate, the random walker stays in the same network layer for about six steps.
Generalized modularity 20 does not identify this type of planted communities across network layers (see Fig. 3D) because it uses a null model only for intra-layer links and merely a coupling parameter between layers 20,28 . As a result, merging different communities across layers always improves the modularity score, as illustrated in Fig. 2C.
Overall, we were not able to recover multilayer commu-nities by treating the multilayer network as one large network, as multiple disconnected networks, or as multiple networks connected with a coupling parameter without a proper null model. We conclude that the key discriminating factor is the map equation's ability to capture the important notion of multiplex networks that sets of state nodes across layers represent the very same physical objects.
Multilayer community structure of collaboration networks We analyzed two inherently multilayer collaboration networks, the Pierre Auger Collaboration of physicists and a sample from the arXiv of researchers working on networks. The Pierre Auger Collaboration is a group of hundreds of theoretical and experimental scientists worldwide working at the largest observatory of ultrahigh energy cosmic rays 41 . The collaborators work together in different research topics on specific tasks, e.g., source detection, mass composition, experimental enhancements, shower reconstruction, etc. Scientists within the Collaboration may work on one or more tasks, and every year hundreds of internal technical reports are submitted to the repository. With access to author lists and keywords, we reconstructed the inherently multilayer collaboration network in which nodes represent scientists, links indicate collaboration between scientists, and layers represent tasks (see Table 1 in the SI). We considered all submissions between 2010 and 2012, and assigned each report to L layers according to its keywords and its content, with manual disambiguation to avoid spurious results from an automated process. For each report with more than one author, for each layer in which the report was classified, and between any pair of the N co-authors, we added a weight 1/(L(N − 1)) to the weighted, undirected, and multilayer network (see Table 1 for details). In this way, the sum of all link weights of an author across all layers simply is the number of reports written by the author. We built the arXiv multilayer network in the very same way, but instead of tasks we used arXiv categories for layers (see Table 2 in the SI). To restrict the analysis to a well-defined topic of research, we only included papers with "networks" in the title or abstract (see Table 1 for details). Because some categories or tasks are more related than others, communities naturally emerge across layers when groups of scientists work on interdisciplinary projects or several tasks simultaneously.
The collaboration networks show a highly overlapping modular organization. In Fig. 4A, we show the largest connected component of the Auger network, including more than 90% of the scientists, and their assignments into highly overlapping modules. Truly multilayer nodes, i.e., those ones corresponding to scientists active in more than one task, dominate the core of the network in this visualization, whereas single-task scientists are more peripherals nodes. For example, the multilayer analysis  reveals strong groups of collaboration across the tasks of "point source," "anisotropy," and "magnetic," (see Fig. 5A). Fig. 4B shows that essential information about the overlapping modular organization is washed out when dynamics are modeled with r = 1.0 or in the aggregated network (not shown in the figure because, qualitatively, it provides the same results), and scientists are assigned to a few overlapping communities (r = 1.0) or one community only (aggregated network). Without mentioning their names here, we find scientists who indisputable are active in several different tasks with variegated collaboration patterns captured only when dynamics are modeled with r = 0.15, whereas for r = 1.0 the scientists are grouped in single non-overlapping communities. In another case, we find two colleagues who work at nearby institutions within the same city and with highly overlapping interests and collaborations. For r = 0.15, they are assigned to highly overlapping modules across tasks, whereas for r = 1.0, they are assigned to different non-overlapping partitions. Only by maintaining the multilayer structure were we able to reveal the actual collaboration structure. Similarly, Fig. 5B shows that communities extend across layers also in the arXiv collaboration network. Whereas communities typically only extend a few layers in the Auger network, communities in the arXiv network can extend over multiple layers. This means that scientists are rather task specific in the Pierre Auger Collaboration, whereas researchers working on networks often are involved in interdisciplinary projects, although computer vision and mathematics seem to be less interdisciplinary topics. In any case, the multilayer networks analyzed with the map equation captures that scientists can simultaneously work in different groups on different topics. Table 1 summarizes the multilayer effects of community detection with the map equation framework. For easy comparison, we contrast multilayer results obtained with dynamics modeled with relax rate r = 0.15, with results obtained with relax rate r = 1.0 (see Fig. 3 in the SI for full comparison). The latter maximum relax rate corresponds to completely washed out multilayer information, but, unlike the aggregated networks, it allows Multiplex Infomap to assign nodes to multiple modules. For both the Auger and the arXiv networks, we find that flow are confined in smaller and more overlapping modules. We also measure this effect in terms of the persistence gain in modules. For modules obtained with r = 0.15, the persistence gain quantifies how much longer a random walker stays within the modules with dynamics modeled with r = 0.15 compared with r = 1.0. When a random walker only moves freely between layers in one step out of about six compared with free movements between layers in any step, we find that its chance to stay within the same module increases by 25 and 13 percent in the Auger and arXiv network, respectively. As a result of this persistence gain, the modular description of a random walker's trajectory can be significantly com-pressed in both networks. Since compressing data is dual to finding regularities in the data 32,42 , the multiplex map equation applied to the multilayer representation allows us to discover patterns that are absent in the aggregated network. Evidently, these patterns contain essential information about the constraints on flow through the systems.
In summary, compared with conventional network analysis, Multiplex Infomap applied to the studied multilayer networks reveals smaller modules with more overlap that better capture the actual organization. Shoehorning multiplex networks into conventional communitydetection algorithms can obscure important structural information and earlier attempts of generalizing conventional community-detection methods to identify communities across layers have proven difficult. In contrast, thanks to the map equation's intrinsic ability to capture that sets of nodes across layers represent the very same physical objects in multiplex networks, the framework generalizes straightforwardly. In absence of empirical inter-layer links, here we have modeled the dynamics between layers. However, inter-layer interaction data would provide further important information about the organization of social and biological systems, and calls for more empirical work.

Dynamics on multilayer networks
The rationale behind the multiplex map equation is simple: encode the trajectory between physical nodes of a random walker that itself navigates between state nodes in different layers (see Fig. 1C). For a modular description, the trajectory is encoded with unique codewords on all modules and all physical nodes within each module, respectively. We are only interested in the codelength and can derive them from the stationary distribution of the random walker. The stationary distribution on the state nodes can be derived from the transition probabilities P αβ ij described in Eq. (1) for interconnected networks with empirical inter-layer link weights and in Eq. (2) for multilayer networks with inter-layer link weights modeled with relaxation parameter r. Assuming that the stationary distribution of state node i, α is p α i , it can in principle be derived from the recursive system of equations However, to guarantee a unique ergodic solution in directed networks, we use teleportation at a low rate τ to state nodes proportional to their intra-layer instrength. 43 . To reduce the smoothening effect of teleportation and make the results more robust to the teleportation parameter τ , we use unrecorded teleportation steps and recorded steps along links 43 . We obtain the recorded visit rates by first calculating the stationary distribution with teleportation to state nodes proportional to their out-strength, with the power-iteration method 44 . Then we derive the recorded steps along links q βα ji and nodes p α i in a subsequent step We use teleportation rate τ = 0.15 throughout our analysis of directed networks, but the results are robust to variation of τ in a wide range around this value. For undirected networks, results are independent of τ .

The multiplex map equation
We seek to minimize the description length L(M) given by the the map equation over possible network partitions M, with each state node i, α assigned to a mod-ule ı = 1, 2, . . . , m. The network partition that gives the shortest description length best captures the community structure with respect to the dynamics on the multilayer network. The map equation uses m module codebooks and one index codebook to describe the random walker's movements within and between modules, respectively, and weights them by how often they are used. From Shannon's source coding theorem 32 , the average description length of each codebook is given by the entropy H(·) of the associated probability distribution of events. Therefore, both the description length and the rate of use of each codebook can be expressed in terms of the visit rates of the state nodes p α i and the transition rates at which the random walker enters q ı and exits q ı each module.
Module codebook ı has one codeword for all state nodes of each physical node assigned to the module and one exit codeword. The codeword lengths are derived from the rates at which the random walker visits each of the physical nodes in the module, and exits the module, q ı . We use p ı to denote the sum of these rates, and P ı = {p i∈ı /p ı } to denote the normalized probability distribution. Similarly, the index codebook has codewords for module entries. The codeword lengths are derived from rates at which the random walker enters each module, q ı . We use q to denote the sum of these rates, and Q = {q ı /q } to denote the normalized probability distribution. We want to express average length of codewords from the index codebook and the module codebooks weighted by their rates of use. Therefore, the map equation is This is the standard formulation of the map equation 15 with one important difference: state nodes of a physical node can be assigned to multiple modules, but if they are assigned to the same module they are assigned a common codeword derived from their total visit rate given by Eq. (9).

OUTLINE
The supporting information is organized into two sections. First, we provide additional information about the synthetic benchmark graphs. In particular, we study the impact of the relax rate r on a special set of multilayer networks with different fraction of overlapping communities across layers. Moreover, we provide additional details about the definition of Normalized Mutual Information for multilayer networks. Finally, in Second, we provide additional information about the real datasets that we study. We have made code available at www.mapequation.org

Community structure of synthetic multilayer networks
In this section we provide additional information about the synthetic benchmark graphs and associated performance tests.

NMI for multilayer networks
As explained in the main text, we use Normalized Mutual Information to assess the similarity of the multilayer partitions. Let us call X the community assignments of one partition and Y the ones of the other partition. The community assignments are the clusters that the nodelayer tuples belong to.
If we draw a tuple at random (with uniform probability), the probability of observing a certain community x is proportional to the number of tuples assigned to it: p(X = x) = n x /N , where n x is the number of tuples in community x and N is the total number of tuples. We can also define the joint probability p(x, y), which is proportional to the number of tuples assigned to community x in one partition and community y in the other.
For the NMI, we used the following definition: where H is the Shannon entropy. As a final remark, in our synthetic benchmark graphs, when the number of layers is high, it can happen that some tuples are never sampled. As a consequence, some tuples are in the reference partition but cannot be in the partition returned by the algorithm. To resolve this issue, we only consider the tuples that appear in both partitions.

Parameters for the LFR benchmark graphs
Here we give further details about how the LFR benchmark networks 37 have been generated. As mentioned in the main text, our synthetic multiplexes comprises T independent LFR benchmark networks, which we call modes. Each of these networks is generated specifying a certain number of parameters, which set the degree distribution, the community sizes and how many links are inter and intra the communities. Inside each community, the links are randomly inserted via the configuration model 45 , and the same model is used to insert the links between the communities.
In the main text, we showed the results obtained for 128 nodes and 4 communities of 32 nodes. The degree of all nodes was set equal to 16, and the mixing parameter, i.e. the average fraction of inter-community links, was set to 5%. These parameters are the same of the popular GN benchmark 36 . However, we also analysed networks of 1, 000 nodes with heterogenous degree and community size distribution, and found the same qualitative behaviour.

Benchmarks with overlapping communities across layers
We consider a set of synthetic networks where each layer consists of 256 nodes grouped in 8 clique communities poorly interconnected with each other. The communities in all layers are assigned in order to obtain a specific fraction of overlapping nodes across layers, i.e., sub-sets of nodes connected with each other in different layers. In our numerical experiments, we consider different realizations with overlapping fraction ranging from 0 to 1. Moreover, we consider different values of the relax rate, ranging from 0 to 1, to understand how the interplay between structure (i.e., overlapping fraction) and dynamics (i.e., relax rate) affects the detection of the planted partitions.
In Fig. S1, we show the resulting phase diagrams for a network with two layers, reporting the number of modules and the resulting NMI of state nodes (Multiplex) and averaged across network layers (Average) against both overlapping fraction and relax rate. First, it is interesting to note that for a wide range of small overlapping fractions and relax rates, the number of detected modules is 16. For overlapping fraction smaller than 50 percent, the networks in the two layers are significantly different. For relax rate lower than 50-60 percent, the layers do not couple and the flow stays preferentially in the cliques within each layer separately. In this region, the two layers behave like if they are not part of the same multiplex and the NMI (Multiplex) estimator is able to detect this behavior, as shown in the right panel of the figure. Conversely, the average NMI suggests that the found partitions per layer are correct: this is equivalent to perform the standard community detection in each layer separately. Outside this region, the two layers are more coupled, because the relax rate is sufficient to allow more information to flow within them and, regardless of overlapping fraction, the number of detected modules is 8.
For overlapping fraction larger than 50%, the NMI (Multiplex) is 1 almost regardless of the relax rate. The average NMI is also 1 in the same range, suggesting that the two measures are equivalent in absence of multilayer communities as defined in the main text (see Fig. 2) and in presence of high overlapping fraction across layers. This result is easily understood in terms of persistent flows across layers and shows that for these multilayer networks, with this specific topology, the average NMI is a more suitable indicator of similarity between the planted partitions and the detected ones.
Moreover, regardless of the overlapping fraction of the underlying multilayer network, 8 modules are always revealed from the analysis of the aggregated network. This result shows the limitations of community detection in aggregated networks: only when the two layers have highly overlapping partitions is the aggregated network a good proxy for the whole multilayer structure. However, this is not the case for the majority of empirical multilayer networks and aggregation could cause a significant loss of information, often leading to a misleading partitioning of the network.

Effect of the relax rate
For multilayer benchmark networks considered in the main text, the optimal range of relax rates depends on the number of mode networks. For example, with only one mode network and relax rate r = 1, it is possible to accurately identify the communities in the mode network for any number of network layers (see Fig. S3). For more than one mode network, on the other hand, too high relax rate washes out the constraints set by each mode network, whereas too low relax rate overstates the constraints set by each network layer. Without access to actual inter-layer link weights, the relax rate should be chosen appropriately for the system under study. Here we show the NMI for the synthetic benchmarks as a function of the relax rate for a single-mode system: T = 1. We already showed the same diagram for a threemode system, T = 3, in Fig. 3 in the main paper (also shown here for comparison). If there is a single mode, the optimal solution is achieved for high values of the relax rate, because the data are aggregated. However, as already shown, aggregating the data cannot find the correct partition if multiple modes are present.
A T=1 B T=3 Figure S2: NMI for the synthetic benchmark graphs as a function of the relax rate r. (A) Performance on a single state system: higher relax rates have a better performance because are similar to the aggregated. (B) For a three state system, only the multilayer solution can detect the correct partition, as the aggregated returns a single module.

Modularity Optimization
We performed the same performance tests on generalized modularity 20 and Fig. 3D in the main text shows that optimizing modularity can only accurately identify communities in each network layer separately in this specific benchmark test. For no value of inter-layer coupling DX did the method accurately identify the multilayer communities of the mode networks. Specifically, for DX < 1, the average NMI shows that the method is capable of detecting the correct partition of each layer, but the multiplex NMI shows that it is not capable of iden-tifying which network layers correspond to which mode networks. For DX > 1, the performance drops further, because inter-layer link weights dominate over intra-layer link weights. We conclude that generalized modularity does not identify this type of multilayer communities, because it uses a null model only for intra-layer links and merely a coupling parameter between layers 20,28 . As a result, merging different communities across layers always improves the modularity score, as illustrated in Fig. 3C.
Comparing the aggregated and multilayer networks of Auger and ArXiv Figure S3 shows how differences between the partitions found with the multilayer and the aggregated approach. At increasing relax rate, the random walker becomes less and less localized in a specific layer. Accordingly, the NMI between the multilayer and the aggregated solutions increases. For r = 1, the walker moves freely between layers, but the NMI is not one because the multilayer solution still allows for overlap (the optimization algorithm we used on the aggregated does not identify overlaps by construction).
In both datasets, with increasing relax rate, we find bigger modules (module size increases) and fewer community assignments per physical node (overlap decreases). The module size is defined as the number of nodes divided by the effective number of modules. The effective number of modules is defined as 2 H , where H is the Shannon entropy of the partition (see previous Section).
A Auger C ArXiv B Auger D ArXiv Figure S3: Aggregation is responsible for significant changes in the community structure. Module sizes, number of assignments per node (overlap) and NMI for communities revealed from the multilayer and the aggregated networks in the Pierre Auger Collaboration (top panels) and the ArXiv (bottom panels) networks. For a given relax rate, the NMI measures the similarity between the obtained partition and the partitions obtained from the aggregated topology (blue curve) and the aggregated dynamics at relax rate r = 1.0 (red curve), respectively.

The Pierre Auger Collaboration Dataset
We considered all internal technical reports submitted to the Pierre Auger 2 GAP repository between 2010 and 2012, assigning each report to one or more tasks according to its keywords and its content and performing disambiguation manually, to avoid unavoidable spurious results due to automated processes. The final number of distinct authors in this dataset is 514, with 9,209 collaborations, classified in 16 layers. Finally, we built the international co-authorship network for each layer, excluding only those papers with only one author to avoid self-loops in the corresponding adjacency matrices. Privacy policies have been considered and we anonymized the data by assigning a random numerical integer to each author. The list of layers is shown in Table S1 together with the corresponding tasks.
Detailed community structure.
We show in Fig. S4 and Fig. S5 more detailed maps of the community structures shown in the main text for the multilayer and the aggregated networks, respectively. Figure S4: Community detection in the Pierre Auger Collaboration network. Detailed map of the partitions obtained applying the map equation with relax rate r = 0.15 to the multilayer network. The size of a node is proportional to the multilayer activity of the corresponding author: larger the number of tasks where he or she collaborates larger the size of the node. Colors within the pie chart code the different communities where the author is classified into, where the area of each slice is proportional to the number of times the author is classified in the corresponding community.