Loading [MathJax]/jax/output/HTML-CSS/jax.js
Full
Published Online: March 2013
Accepted: January 2013
Chaos 23, 013142 (2013); doi: http://dx.doi.org/10.1063/1.4790830

We describe techniques for the robust detection of community structure in some classes of time-dependent networks. Specifically, we consider the use of statistical null models for facilitating the principled identification of structural modules in semi-decomposable systems. Null models play an important role both in the optimization of quality functions such as modularity and in the subsequent assessment of the statistical validity of identified community structure. We examine the sensitivity of such methods to model parameters and show how comparisons to null models can help identify system scales. By considering a large number of optimizations, we quantify the variance of network diagnostics over optimizations (“optimization variance”) and over randomizations of network structure (“randomization variance”). Because the modularity quality function typically has a large number of nearly degenerate local optima for networks constructed using real data, we develop a method to construct representative partitions that uses a null model to correct for statistical noise in sets of partitions. To illustrate our results, we employ ensembles of time-dependent networks extracted from both nonlinear oscillators and empirical neuroscience data.
Many social, physical, technological, and biological systems can be modeled as networks composed of numerous interacting parts.11. M. E. J. Newman, Networks: An Introduction (Oxford University Press, 2010). As an increasing amount of time-resolved data has become available, it has become increasingly important to develop methods to quantify and characterize dynamic properties of temporal networks.22. P. Holme and J. Saramäki, Phys. Rep. 519, 97 (2012). https://doi.org/10.1016/j.physrep.2012.03.001 Generalizing the study of static networks, which are typically represented using graphs, to temporal networks entails the consideration of nodes (representing entities) and/or edges (representing ties between entities) that vary in time. As one considers data with more complicated structures, the appropriate network analyses must become increasingly nuanced. In the present paper, we discuss methods for algorithmic detection of dense clusters of nodes (i.e., communities) by optimizing quality functions on multilayer network representations of temporal networks.3,43. P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, Science 328, 876 (2010). https://doi.org/10.1126/science.1184819 4. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 We emphasize the development and analysis of different types of null-model networks, whose appropriateness depends on the structure of the networks one is studying as well as the construction of representative partitions that take advantage of a multilayer network framework. To illustrate our ideas, we use ensembles of time-dependent networks from the human brain and human behavior.
Myriad systems have components whose interactions (or the components themselves) change as a function of time. Many of these systems can be investigated using the framework of temporal networks, which consist of sets of nodes and/or edges that vary in time.22. P. Holme and J. Saramäki, Phys. Rep. 519, 97 (2012). https://doi.org/10.1016/j.physrep.2012.03.001 The formalism of temporal networks is convenient for studying data drawn from areas such as person-to-person communication (e.g., via mobile phones5,65. J.-P. Onnela, J. Saramäki, J. Hyvönen, G. Szabó, D. Lazer, K. Kaski, J. Kertész, and A. L. Barabási, Proc. Natl. Acad. Sci. U.S.A. 104, 7332 (2007). https://doi.org/10.1073/pnas.0610245104 6. Y. Wu, C. Zhou, J. Xiao, J. Kurths, and H. J. Schellnhuber, Proc. Natl. Acad. Sci. U.S.A. 107, 18803 (2010). https://doi.org/10.1073/pnas.1013140107 ), one-to-many information dissemination (such as Twitter networks77. S. Gonzaález-Bailón, J. Borge-Holthoefer, A. Rivero, and Y. Moreno, Sci. Rep. 1, 197 (2011). https://doi.org/10.1038/srep00197 ), cell biology, distributed computing, infrastructure networks, neural and brain networks, and ecological networks.22. P. Holme and J. Saramäki, Phys. Rep. 519, 97 (2012). https://doi.org/10.1016/j.physrep.2012.03.001 Important phenomena that can be studied in this framework include network constraints on gang and criminal activity,8,98. T. J. Fararo and J. Skvoretz, in Status, Network, and Structure: Theory Development in Group Processes (Stanford University Press, 1997), pp. 362–386. 9. A. Stomakhin, M. B. Short, and A. L. Bertozzi, Inverse Probl. 27, 115013 (2011). https://doi.org/10.1088/0266-5611/27/11/115013 political processes,10,1110. P. J. Mucha and M. A. Porter, Chaos 20, 041108 (2010). https://doi.org/10.1063/1.3518696 11. A. S. Waugh, L. Pei, J. H. Fowler, P. J. Mucha, and M. A. Porter, “ Party polarization in Congress: A network science approach,” arXiv:0907.3509 (2011). human brain function,4,124. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 12. K. W. Doron, D. S. Bassett, and M. S. Gazzaniga, Proc. Natl. Acad. Sci. U.S.A. 109, 18661 (2012). https://doi.org/10.1073/pnas.1216402109 human behavior,1313. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 and financial structures.14,1514. D. J. Fenn, M. A. Porter, M. McDonald, S. Williams, N. F. Johnson, and N. S. Jones, Chaos 19, 033119 (2009). https://doi.org/10.1063/1.3184538 15. D. J. Fenn, M. A. Porter, S. Williams, M. McDonald, N. F. Johnson, and N. S. Jones, Phys. Rev. E 84, 026109 (2011). https://doi.org/10.1103/PhysRevE.84.026109
Time-dependent complex systems can have densely connected components in the form of cohesive groups of nodes known as “communities” (see Fig. 1), which can be related to a system's functional modules.16,1716. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 17. H. Simon, Proc. Am. Philos. Soc. 106, 467 (1962). A wide variety of clustering techniques have been developed to identify communities, and they have yielded insights in the study of the committee structure in the United States Congress,1818. M. A. Porter, P. J. Mucha, M. Newman, and C. M. Warmbrand, Proc. Natl. Acad. Sci. U.S.A. 102, 7057 (2005). https://doi.org/10.1073/pnas.0500191102 functional groups in protein interaction networks,1919. A. C. F. Lewis, N. S. Jones, M. A. Porter, and C. M. Deane, BMC Syst. Biol. 4, 100 (2010). https://doi.org/10.1186/1752-0509-4-100 functional modules in brain networks,44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 and more. A particularly successful technique for identifying communities in networks16,2016. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 20. S. Fortunato, Phys. Rep. 486, 75 (2010). https://doi.org/10.1016/j.physrep.2009.11.002 is optimization of a quality function known as “modularity,”2121. M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004). https://doi.org/10.1103/PhysRevE.69.026113 which recently has been generalized for detecting communities in time-dependent and multiplex networks.33. P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, Science 328, 876 (2010). https://doi.org/10.1126/science.1184819
Modularity optimization allows one to algorithmically partition a network's nodes into communities such that the total connection strength within groups of the partition is more than would be expected in some null model. However, modularity optimization always yields a network partition (into a set of communities) as an output whether or not a given network truly contains modular structure. Therefore, application of subsequent diagnostics to a network partition is potentially meaningless without some comparison to benchmark or null-model networks. That is, it is important to establish whether the partition(s) obtained appear to represent meaningful community structures within the network data or whether they might have reasonably arisen at random. Moreover, robust assessment of network organization depends fundamentally on the development of statistical techniques to compare structures in a network derived from real data to those in appropriate models (see, e.g., Ref. 2222. A. Lancichinetti, F. Radicchi, and J. J. Ramasco, Phys. Rev. E 81, 046110 (2010). https://doi.org/10.1103/PhysRevE.81.046110 ). Indeed, as the constraints in null models and network benchmarks become more stringent, it can become possible to make stronger claims when interpreting organizational structures such as community structure.
In the present paper, we examine null models in time-dependent networks and investigate their use in the algorithmic detection of cohesive, dynamic communities in such networks (see Fig. 2). Indeed, community detection in temporal networks necessitates the development of null models that are appropriate for such networks. Such null models can help provide bases of comparison at various stages of the community-detection process, and they can thereby facilitate the principled identification of dynamic structure in networks. Indeed, the importance of developing null models extends beyond community detection, as such models make it possible to obtain statistically significant estimates of network diagnostics.
Our dynamic network null models fall into two categories: optimization null models, which we use in the identification of community structure; and post-optimization null models, which we use to examine the identified community structure. We describe how these null models can be selected in a manner appropriate to known features of a network's construction, identify potentially interesting network scales by determining values of interest for structural and temporal resolution parameters, and inform the choice of representative partitions of a network into communities.
A. Community detection
Community-detection algorithms provide ways to decompose a network into dense groups of nodes called “modules” or “communities.” Intuitively, a community consists of a set of nodes that are connected among one another more densely than they are to nodes in other communities. A popular way to identify community structure is to optimize a quality function, which can be used to measure the relative densities of intra-community connections versus inter-community connections. See Refs. 16, 20, and 2316. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 20. S. Fortunato, Phys. Rep. 486, 75 (2010). https://doi.org/10.1016/j.physrep.2009.11.002 23. M. E. J. Newman, Nat. Phys. 8, 25 (2012). https://doi.org/10.1038/nphys2162 for recent reviews on network community structure and Refs. 24–2724. B. H. Good, Y. A. de Montjoye, and A. Clauset, Phys. Rev. E 81, 046106 (2010). https://doi.org/10.1103/PhysRevE.81.046106 25. S. Fortunato and M. Barthélemy, Proc. Natl. Acad. Sci. U.S.A. 104, 36 (2007). https://doi.org/10.1073/pnas.0605965104 26. P. J. Bickel and A. Chen, Proc. Natl. Acad. Sci. U.S.A. 106, 21068 (2009). https://doi.org/10.1073/pnas.0907096106 27. S. Hutchings, “ The behavior of modularity-optimizing community detection algorithms,” M.Sc. Thesis (University of Oxford, 2011). for discussions of various caveats that should be considered when optimizing quality functions to detect communities.
One begins with a network of N nodes and a given set of connections between those nodes. In the usual case of single-layer networks (e.g., static networks with only one type of edge), one represents a network using an N×N adjacency matrix A. The element Aij of the adjacency matrix indicates a direct connection or “edge” from node i to node j, and its value indicates the weight of that connection. The quality of a hard partition of A into communities (whereby each node is assigned to exactly one community) can be quantified using a quality function. The most popular choice is modularity16,20,21,28,2916. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 20. S. Fortunato, Phys. Rep. 486, 75 (2010). https://doi.org/10.1016/j.physrep.2009.11.002 21. M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004). https://doi.org/10.1103/PhysRevE.69.026113 28. M. E. J. Newman, Phys. Rev. E 69, 066133 (2004). https://doi.org/10.1103/PhysRevE.69.066133 29. M. E. J. Newman, Proc. Natl. Acad. Sci. U.S.A. 103, 8577 (2006). https://doi.org/10.1073/pnas.0601602103
Q0=ij[AijγPij]δ(gi,gj),(1)
where node i is assigned to community gi, node j is assigned to community gj, the Kronecker delta δ(gi,gj)=1 if gi=gj and it equals 0 otherwise, γ is a resolution parameter (which we will call a structural resolution parameter), and Pij is the expected weight of the edge connecting node i to node j under a specified null model. The choice γ=1 is very common, but it is important to consider multiple values of γ to examine groups at multiple scales.16,30,3116. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 30. J. Reichardt and S. Bornholdt, Phys. Rev. E 74, 016110 (2006). https://doi.org/10.1103/PhysRevE.74.016110 31. J.-P. Onnela, D. J. Fenn, S. Reid, M. A. Porter, P. J. Mucha, M. D. Fricker, and N. S. Jones, Phys. Rev. E. 86, 036104 (2012). https://doi.org/10.1103/PhysRevE.86.036104 Maximization of Q0 yields a hard partition of a network into communities such that the total edge weight inside of modules is as large as possible (relative to the null model and subject to the limitations of the employed computational heuristics, as optimizing Q0 is NP-hard16,20,3216. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 20. S. Fortunato, Phys. Rep. 486, 75 (2010). https://doi.org/10.1016/j.physrep.2009.11.002 32. U. Brandes, D. Delling, M. Gaertler, R. Görke, M. Hoefer, Z. Nikoloski, and D. Wagner, IEEE Trans. Knowl. Data Eng. 20, 172 (2008). https://doi.org/10.1109/TKDE.2007.190689 ).
Recently, the null model in the quality function (1) has been generalized so that one can consider sets of L adjacency matrices, which are combined to form a rank-3 adjacency tensor A that can be used to represent time-dependent or multiplex networks. One can thereby define a multilayer modularity (also called “multislice modularity”)33. P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, Science 328, 876 (2010). https://doi.org/10.1126/science.1184819
Q=12μijlr{(AijlγlPijl)δlr+δijωjlr}δ(gil,gjr),(2)
where the adjacency matrix of layer l has components Aijl, the element Pijl gives the components of the corresponding layer-l matrix for the optimization null model, γl is the structural resolution parameter of layer l, the quantity gil gives the community assignment of node i in layer l, the quantity gjr gives the community assignment of node j in layer r, the element ωjlr gives the connection strength (i.e., an “interlayer coupling parameter,” which one can call a temporal resolution parameter if one is using the adjacency tensor to represent a time-dependent network) from node j in layer r to node j in layer l, the total edge weight in the network is μ=12jrκjr, the strength (i.e., weighted degree) of node j in layer l is κjl=kjl+cjl, the intra-layer strength of node j in layer l is kjl=iAijl, and the inter-layer strength of node j in layer l is cjl=rωjlr.
Equivalent representations that use other notation can, of course, be useful. For example, multilayer modularity can be recast as a set of rank-2 matrices describing connections between the set of all nodes across layers [e.g., for spectral partitioning29,33,3429. M. E. J. Newman, Proc. Natl. Acad. Sci. U.S.A. 103, 8577 (2006). https://doi.org/10.1073/pnas.0601602103 33. M. E. J. Newman, Phys. Rev. E 74, 036104 (2006). https://doi.org/10.1103/PhysRevE.74.036104 34. T. Richardson, P. J. Mucha, and M. A. Porter, Phys. Rev. E 80, 036111 (2009). https://doi.org/10.1103/PhysRevE.80.036111 ]. One can similarly generalize Q for higher-rank tensors, which one can use when studying community structure in networks that are both time-dependent and multiplex, through appropriate specification of inter-layer coupling tensors.
B. Network diagnostics
To characterize multilayer community structure, we compute four example diagnostics for each hard partition: the modularity Q, the number of modules n, the mean community size s (which is equal to the number of nodes in the community and is proportional to 1/n), and the stationarity ζ.3535. G. Palla, A. Barabási, and T. Vicsek, Nature 446, 664 (2007). https://doi.org/10.1038/nature05670 To compute ζ, we calculate the autocorrelation function U(t, t + m) of two states of the same community G(t) at m time steps (i.e., m network layers) apart
U(t,t+m)|G(t)G(t+m)||G(t)G(t+m)|,(3)
where |G(t)G(t+m)| is the number of nodes that are members of both G(t) and G(t + m), and |G(t)G(t+m)| is the number of nodes in the union of the community at times t and t + m. Defining t0 to be the first time step in which the community exists and t to be the last time in which it exists, the stationarity of a community is3535. G. Palla, A. Barabási, and T. Vicsek, Nature 446, 664 (2007). https://doi.org/10.1038/nature05670
ζt1t=t0U(t,t+1)tt0.(4)
This gives the mean autocorrelation over consecutive time steps.3636. Equation (4) gives the definition for the notion of stationarity that we used in Ref. 44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 . The equation for this quantity in Ref. 44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 has a typo in the denominator. We wrote incorrectly in that paper that the denominator was tt01, whereas the numerical computations in that paper used tt0.
In addition to these diagnostics, which are defined using the entire multilayer community structure, we also compute two example diagnostics on the community structures of the component layers: the mean single-layer modularity Qs and the variance var(Qs) of the single-layer modularity over all layers. The single-layer modularity Qs is defined as the static modularity quality function, Qs=ij[AijγPij]δ(gi,gj), computed for the partition g that we obtained via optimization of the multilayer modularity function Q. We have chosen to use a few simple ways to help characterize the time series for Qs, though of course other diagnostics can also be informative.
C. Data sets
We illustrate dynamic network null models using two example network ensembles: (1) 75-time-layer brain networks drawn from each of 20 human subjects and (2) behavioral networks with about 150 time layers drawn from each of 22 human subjects. Importantly, the use of network ensembles makes it possible to examine robust structure (and also its variance) over multiple network instantiations. We have previously examined both data sets in the context of neuroscientific questions.4,134. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 13. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 In this paper, we use them as illustrative examples for the consideration of methodological issues in the detection of dynamic communities in temporal networks.
These two data sets, which provide examples of different types of network data, illustrate a variety of issues in network construction: (1) node and edge definitions, (2) complete versus partial connectivity, (3) ordered versus categorical nodes, and (4) confidence in edge weights. In many fields, determining the definition of nodes and edges is itself an active area of investigation.3737. In other areas of investigation, it probably should be. See, for example, several recent papers that address such questions in the context of large-scale human brain networks38–4338. E. T. Bullmore and D. S. Bassett, Ann. Rev. Clin. Psych. 7, 113 (2011). https://doi.org/10.1146/annurev-clinpsy-040510-143934 39. D. S. Bassett, J. A. Brown, V. Deshpande, J. M. Carlson, and S. T. Grafton, Neuroimage 54, 1262 (2011). https://doi.org/10.1016/j.neuroimage.2010.09.006 40. A. Zalesky, A. Fornito, I. H. Harding, L. Cocchi, M. Yucel, C. Pantelis, and E. T. Bullmore, Neuroimage 50, 970 (2010). https://doi.org/10.1016/j.neuroimage.2009.12.027 41. J. Wang, L. Wang, Y. Zang, H. Yang, H. Tang, Q. Gong, Z. Chen, C. Zhu, and Y. He, Hum. Brain Mapp 30, 1511 (2009). https://doi.org/10.1002/hbm.20623 42. S. Bialonski, M. T. Horstmann, and K. Lehnertz, Chaos 20, 013134 (2010). https://doi.org/10.1063/1.3360561 43. A. A. Ioannides, Curr. Opin. Neurobiol. 17, 161 (2007). https://doi.org/10.1016/j.conb.2007.03.008 and in networks more generally.4444. C. T. Butts, Science 325, 414 (2009). https://doi.org/10.1126/science.1171022 Another important issue is whether to examine a given adjacency matrix in an exploratory manner or to impose structure on it based on a priori knowledge. For example, when nodes are categorical, one might represent their relations using a fully connected network and then identify communities of any group of nodes. However, when nodes are ordered—and particularly when they are in a chain of weighted nearest-neighbor connections—one expects communities to group neighboring nodes in sequence, as typical community-detection methods are unlikely to yield many out-of-sequence jumps in community assignment. The issue of confidence in the estimation of edge weights is also very important, as it can prompt an investigator to delete edges from a network when their statistical validity is questionable. A closely related issue is how to deal with known or expected missing data, which can affect either the presence or absence of nodes themselves or the weights of edges.45–4845. G. Kossinets, Soc. Networks 28, 247 (2006). https://doi.org/10.1016/j.socnet.2005.07.002 46. A. Clauset, C. Moore, and M. E. J. Newman, Nature 453, 98 (2008). https://doi.org/10.1038/nature06830 47. R. Guimerá and M. Sales-Pardo, Proc. Natl. Acad. Sci. U.S.A. 106, 22073 (2009). https://doi.org/10.1073/pnas.0908366106 48. M. Kim and J. Leskovec, SIAM International Conference on Data Mining (2011).
D. Data set 1: Brain networks
Our first data set contains categorical nodes with partial connectivity and variable confidence in edge weights. The nodes remain unchanged in time, and edge weights are based on covariance of node properties. This covariance structure is non-local in the sense that weights exist between both topologically neighboring nodes and topologically distant nodes.49,5049. O. Sporns, Networks of the Brain (MIT, 2010). 50. E. Bullmore and O. Sporns, Nat. Rev. Neurosci. 10, 186 (2009). https://doi.org/10.1038/nrn2575 This property has been linked in other dynamical systems to behaviors such as chimera states, in which coherent and incoherent regions coexist.51–5351. D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004). https://doi.org/10.1103/PhysRevLett.93.174102 52. Y. Kuramoto and D. Battogtokh, Nonlinear Phenom. Complex Syst. 5, 380 (2002). 53. S. I. Shima and Y. Kuramoto, Phys. Rev. E 69, 036213 (2004). https://doi.org/10.1103/PhysRevE.69.036213 Another interesting feature of this data set is that it is drawn from an experimental measurement with high spatial resolution (on the order of centimeters) but relatively poor temporal resolution (on the order of seconds).
As described in more detail in Ref. 44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 , we construct an ensemble of networks (20 individuals over 3 experiments, which yields 60 multilayer networks) that represent the functional connectivity between large regions of the human brain. In these networks, N = 112 centimeter-scale, anatomically distinct brain regions are our (categorical) network nodes. We study the temporal interaction of these nodes—such interactions are thought to underly cognitive function—by first measuring their activity every 2 s during simple finger movements using functional magnetic resonance imaging (fMRI). We cut these regional time series into time slices (which yield layers in the multilayer network) of roughly 3-min duration. Each such layer corresponds to a time series whose length is 80 units.
To estimate the interactions (i.e., edge weights) between nodes, we calculate a measure of statistical similarity between regional activity profiles.5454. K. J. Friston, Hum. Brain Mapp 2, 56 (1994). https://doi.org/10.1002/hbm.460020107 Using a wavelet transform, we extract frequency-specific activity from each time series in the range 0.06–0.12 Hz. For each time layer l and each pair of regions i and j, we define the weight of an edge connecting region i to region j using the coherence between the wavelet-coefficient time series in each region, and these weights form the elements of a weighted, undirected temporal network W with components Wijl=Wjil. The magnitude-squared coherence Gij between time series i and j is a function of frequency. It is defined by the equation
Gij(f)=|Fij(f)|2Fii(f)Fjj(f),(5)
where Fii(f) and Fjj(f) are the power spectral density functions of i and j, respectively, and Fij(f) is the cross-power spectral density function of i and j. We let Hij denote the mean of Gij(f) over the frequency band of interest, and the weight of edge Wijl is equal to Hij computed for layer l.
We use a false-discovery rate correction5555. Y. Benjamini and Y. Yekutieli, Ann. Stat. 29, 1165 (2001). https://doi.org/10.1214/aos/1013699998 to threshold connections whose coherence values are not significantly greater than that expected at random. This yields a multilayer network A with components Aijl (i.e., a rank-3 adjacency tensor). The nonzero entries in A retain their weights. We couple the layers of A to one another with temporal resolution parameters of weight ωjlr between node j in layer r and node j in layer l. In this paper, we let ωjlrω[0.1,40] be identical between each node j in a given layer and itself in nearest-neighbor layers. (In all other cases, ωjlr=0.)
In Fig. 3(a), we show an example time layer from A for a single subject in this experimental data. In this example, the statistical threshold is evinced by the set of matrix elements set to 0. Because brain network nodes are categorical, one can apply community detection algorithms in these situations to identify communities composed of any set of nodes. (Note that the same node from different times can be assigned to the same community even if the node is assigned to other communities at intervening times.) One biological interpretation of network communities in brain networks is that they represent groups of nodes that serve distinct cognitive functions (e.g., vision, memory, etc.) that can vary in time.12,5612. K. W. Doron, D. S. Bassett, and M. S. Gazzaniga, Proc. Natl. Acad. Sci. U.S.A. 109, 18661 (2012). https://doi.org/10.1073/pnas.1216402109 56. E. Bullmore and O. Sporns, Nat. Rev. Neurosci. 13, 336 (2012). https://doi.org/10.1038/nrn3214
E. Data set 2: Behavioral networks
Our second data set contains ordered nodes that remain unchanged in time. The network topology in this case is highly constrained, as edges are only present between consecutive nodes. (We call this “nearest-neighbor” coupling.) Another interesting feature of this data set is that the number of time slices is an order of magnitude larger than the number of nodes in a slice.
As described in more detail in Ref. 1313. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 , we construct an ensemble of 66 behavioral networks from 22 individuals and 3 experimental conditions. These networks represent a set of finger movements in the same simple motor learning experiment from which we constructed the brain networks in data set 1. Subjects were instructed to press a sequence of buttons corresponding to a sequence of 12 pseudo-musical notes shown to them on a screen.
Each node represents an interval between consecutive button presses. A single network layer consists of N  = 11 nodes (i.e., there is one interval between each pair of notes), which are connected in a chain via weighted, undirected edges. In Ref. 1313. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 , we examined the phenomenon of motor “chunking,” which is a fascinating but poorly understood phenomenon in which groups of movements are made with similar inter-movement durations. (This is similar to remembering a phone number in groups of a few digits or grouping notes together as one masters how to play a song.) For each experimental trial l and each pair of inter-movement intervals i and j, we define the weight of an edge connecting inter-movement i to inter-movement j as the normalized similarity in inter-movement durations. The normalized similarity between nodes i and j is defined as
ρijl=ˉdldijlˉdl,(6)
where dijl is the absolute value of the difference of lengths of the ith and jth inter-movement time intervals in trial l and ˉdl is the maximum value of dijl in trial l. These weights yield the elements Wijl of a weighted, undirected multilayer network W. Because finger movements occur in series, inter-movement i is connected in time to inter-movement i ± 1 but not to any other inter-movements i + n for |n|1.
To encode this conceptual relationship as a network, we set all non-contiguous connections in W to 0 and thereby construct a weighted, undirected chain network A. In Fig. 3(b), we show an example trial layer from A for a single subject in this experimental data. We couple layers of A to one another with weight ωjlr, which gives the connection strength between node j in experimental trial r and node j in trial l. In a given instantiation of the network, we again let ωjlrω[0.1,40] be identical for all nodes j for all connections between nearest-neighbor layers. (Again, ωjlr=0 in all other cases.) Because inter-movement nodes are ordered, one can apply community-detection algorithms to identify communities of nodes in sequence. Each community represents a motor “chunk.”
A. Modularity-optimization null models
After constructing a multilayer network A with elements Aijl, it is necessary to select an optimization null model P in Eq. (2). The most common modularity-optimization null model used in undirected, single-layer networks is the Newman-Girvan null model16,20,21,28,2916. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009). 20. S. Fortunato, Phys. Rep. 486, 75 (2010). https://doi.org/10.1016/j.physrep.2009.11.002 21. M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004). https://doi.org/10.1103/PhysRevE.69.026113 28. M. E. J. Newman, Phys. Rev. E 69, 066133 (2004). https://doi.org/10.1103/PhysRevE.69.066133 29. M. E. J. Newman, Proc. Natl. Acad. Sci. U.S.A. 103, 8577 (2006). https://doi.org/10.1073/pnas.0601602103
Pij=kikj2m,(7)
where ki=jAij is the strength of node i and m=12ijAij. The definition (7) can be extended to multilayer networks using
Pijl=kilkjl2ml,(8)
where kil=jAijl is the strength of node i in layer l and ml=12ijAijl. Optimization of Q using the null model (8) identifies partitions of a network into groups that have more connections (in the case of binary networks) or higher connection densities (in the case of weighted networks) than would be expected for the distribution of connections (or connection densities) expected in a null model. We use the notation Al for the layer-l adjacency matrix composed of elements Aijl and the notation Pl to denote the layer-l null-model matrix with elements Pijl. See Fig. 4(a) for an example layer Al from a multilayer behavioral network and Fig. 4(b) for an example instantiation of the Newman-Girvan null model Pl.
1. Optimization null models for ordered node networks
The Newman-Girvan null model is particularly useful for networks with categorical nodes, in which a connection between any pair of nodes can occur in theory. However, when using a chain network of ordered nodes, it is useful to consider alternative null models. For example, in a network represented by an adjacency matrix A′, one can define
Pij=ρAij,(9)
where ρ is the mean edge weight of the chain network and A′ is the binarized version of A, in which nonzero elements of A are set to 1 and zero-valued elements remain unaltered. Such a null model can also be defined for a multilayer network that is represented by a rank-3 adjacency tensor A. One can construct a null model P with components
Pijl=ρlAijl,(10)
where ρl is the mean edge weight in layer l and A′ is the binarized version of A. The optimization of Q using this null model identifies partitions of a network whose communities have a larger strength than the mean. See Fig. 4(c) for an example of this chain null model Pl for the behavioral network layer shown in Fig. 4(a).
In Fig. 4(d), we illustrate the effect that the choice of optimization null model has on the modularity values Q of the behavioral networks as a function of the structural resolution parameter. (Throughout the manuscript, we use a Louvain-like locally greedy algorithm to maximize the multilayer modularity quality function.57,5857. V. D. Blondel, J. L. Guillaume, R. Lambiotte, and E. Lefebvre, J. Stat. Mech. 10, P10008 (2008). 58. I. S. Jutla, L. G. S. Jeub, and P. J. Mucha, “ A generalized Louvain method for community detection implemented in matlab,” (2011–2012); http://netwiki.amath.unc.edu/GenLouvain. ) The Newman-Girvan null model gives decreasing values of Q for γ[0.1,2.1], whereas the chain null model produces lower values of Q, which behaves in a qualitatively different manner for γ<1 versus γ>1. To help understand this feature, we plot the number and mean size of communities as a function of γ in Figs. 4(e) and 4(f). As γ is increased, the Newman-Girvan null model yields network partitions that contain progressively more communities (with progressively smaller mean size). The number of communities that we obtain in partitions using the chain null model also increases with γ, but it does so less gradually. For γ1, one obtains a network partition consisting of a single community of size Nl=11; for γ1, each node is instead placed in its own community. For γ=1, nodes are assigned to several communities whose constituents vary with time (see, for example, Fig. 3(d)).
The above results highlight the sensitivity of network diagnostics such as Q, n, and s to the choice of an optimization null model. It is important to consider this type of sensitivity in the light of other known issues, such as the extreme near-degeneracy of quality functions like modularity.2424. B. H. Good, Y. A. de Montjoye, and A. Clauset, Phys. Rev. E 81, 046106 (2010). https://doi.org/10.1103/PhysRevE.81.046106 Importantly, the use of the chain null model provides a clear delineation of network behavior in this example into three regimes as a function of γ: a single community with variable Q (low γ), a variable number of communities as Q reaches a minimum value (γ1), and a set of singleton communities with minimum Q (high γ). This illustrates that it is crucial to consider a null model appropriate for a given network, as it can provide more interpretable results than just using the usual choices (such as the Newman-Girvan null model).
The structural resolution parameter γ can be transformed so that it measures the effective fraction of edges ξ(γ) that have larger weights than their null-model counterparts.3131. J.-P. Onnela, D. J. Fenn, S. Reid, M. A. Porter, P. J. Mucha, M. D. Fricker, and N. S. Jones, Phys. Rev. E. 86, 036104 (2012). https://doi.org/10.1103/PhysRevE.86.036104 One can define a generalization of ξ to multilayer networks, which allows one to examine the behavior of the chain null model near γ=1 in more detail. For each layer l, we define a matrix Xl(γ) with elements Xijl(γ)=AijlγPijl, and we then define cX(γ) to be the number of elements of Xl(γ) that are less than 0. We sum cX(γ) over layers in the multilayer network to construct cXml(γ). The transformed structural resolution parameter is then given by
ξml(γ)=cXml(γ)cXml(Λmin)cXml(Λmax)cXml(Λmin),(11)
where Λmin is the value of γ for which the network still forms a single community in the multilayer optimization, and Λmax is the value of γ for which the network still forms N singleton communities in the multilayer optimization. (We use Roman typeface in the subscripts in cXml and ξml to emphasize that we are describing multilayer objects and, in particular, that the subscripts do not represent indices.) In Figs. 4(g)–4(i), we report the optimized (i.e., maximized) modularity value, the number of communities, and the mean community size as functions of the transformed structural resolution parameter ξml(γ). (Compare these plots to Figs. 4(d)–4(f).) For all three diagnostics, the apparent transition points seem to be more gradual as a function of ξml(γ) than they are as a function of γ. For systems like the present one that do not exhibit a pronounced, nontrivial plateau in these diagnostics as a function of a structural resolution parameter, it might be helpful to have a priori knowledge about the expected number or sizes of communities (see, e.g., Ref. 1313. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 ) to help guide further investigation.
2. Optimization null models for networks derived from time series
Although the Newman-Girvan null model can be used in networks with categorical nodes, such as the brain networks in data set 1 (see Fig. 5(a)), it does not take advantage of the fact that these networks are derived from similarities in time series. Accordingly, we generate surrogate data to construct two dynamic network null models for community detection that might be particularly appropriate for networks derived from time-series data.
First, we note that a simple null model (which we call “Random”) for time series is to randomize the elements of the time-series vector for each node before computing the similarity matrix (see Fig. 5(b)).5959. A discrete time series can be represented as a vector. A continuous time series would first need to be discretized. However, the resulting time series do not have the mean or variance of the original time series, and this yields a correlation- or coherence-based network with very low edge weights. To preserve the mean, variance, and autocorrelation function of the original time series, we employ a surrogate-data generation method that scrambles the phase of time series in Fourier space.6060. D. Prichard and J. Theiler, Phys. Rev. Lett. 73, 951 (1994). https://doi.org/10.1103/PhysRevLett.73.951 Specifically, we assume that the linear properties of the time series are specified by the squared amplitudes of the discrete Fourier transform (FT)
|S(u)|2=|1VV1v=0svei2πuv/V|2,(12)
where sv denotes an element in a time series of length V. (That is, V is the number of elements in the time-series vector.) We construct surrogate data by multiplying the Fourier transform by phases chosen uniformly at random and transforming back to the time domain
ˉsv=1VV1v=0eiau|Su|ei2πkv/V,(13)
where au[0,2π) are chosen independently and uniformly at random.6161. The code used for this computation actually operates on [0,2π]. However, this should be an equivalent mathematical estimate to the same computation on [0,2π), which is the same except for a set of measure zero. This method, which we call the FT surrogate (see Fig. 5(c)), has been used previously to construct covariance matrices6262. H. Nakatani, I. Khalilov, P. Gong, and C. van Leeuwen, Phys. Lett. A 319, 167 (2003). https://doi.org/10.1016/j.physleta.2003.09.082 and to characterize networks.6363. A. Zalesky, A. Fornito, and E. Bullmore, Neuroimage 60, 2096 (2012). https://doi.org/10.1016/j.neuroimage.2012.02.001 A modification of this method, which we call the amplitude-adjusted Fourier transform (AAFT) surrogate, allows one to also retain the amplitude distribution of the original signal6464. J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, and J. D. Farmer, Physica D 58, 77 (1992). https://doi.org/10.1016/0167-2789(92)90102-S (see Fig. 5(d)). One can alter nonlinear relationships between time series while preserving linear relationships between time series by applying an identical shuffling to both time series; one can alter both linear and nonlinear relationships between time series by applying independent shufflings to each time series.6060. D. Prichard and J. Theiler, Phys. Rev. Lett. 73, 951 (1994). https://doi.org/10.1103/PhysRevLett.73.951
We demonstrate in Fig. 5(e) that, among the four null models that we consider, the mean coherence of pairs of FT surrogate series match that of the original data most closely. Pairs of Random time series have the smallest mean coherence, and pairs of AAFT surrogate series have the next smallest. The fact that the AAFT surrogate is less like the real data (in terms of mean coherence) than the simpler FT surrogate might stem from a rescaling step6262. H. Nakatani, I. Khalilov, P. Gong, and C. van Leeuwen, Phys. Lett. A 319, 167 (2003). https://doi.org/10.1016/j.physleta.2003.09.082 that causes the power spectrum to whiten (i.e., the step flattens the power spectral density).65,83,8465. It is also important to note that the AAFT method can allow nonlinear correlations to remain in the surrogate data. Therefore, the development of alternative surrogate data generation methods might be necessary (Refs. 83 and 84). 83. C. Räth, M. Gliozzi, I. E. Papadakis, and W. Brinkmann, Phys. Rev. Lett. 109, 144101 (2012). https://doi.org/10.1103/PhysRevLett.109.144101 84. G. Rossmanith, H. Modest, C. Räth, A. J. Banday, K. M. Gorski, and G. Morfill, Phys. Rev. D 86, 083005 (2012). https://doi.org/10.1103/PhysRevD.86.083005 In Figs. 5(f)–5(h), we show three diagnostics (optimized modularity, mean community size, and number of communities) as a function of the structural resolution parameter γ for the various optimization null models. We note that the Newman-Girvan null model produces the smallest Q value and a middling community size, whereas the surrogate time series models produce higher Q values and more communities of smaller mean size. The Random null model produces the largest value of Q and the fewest communities, which is consistent with the fact that it contains the smallest amount of shared information (i.e., mean coherence) with the real network.
B. Post-optimization null models
After identifying the partition(s) that maximize modularity, one might wish to determine whether the identified community structure is significantly different from that expected under other null hypotheses. For example, one might wish to know whether any temporal evolution is evident in the dynamic community structure (see Fig. 2(c)). To do this, one can employ post-optimization null models, in which a multilayer network is scrambled in some way to produce a new multilayer network. One can then maximize the modularity of the new network and compare the resulting community structure to that obtained using the original network. Unsurprisingly, one's choice of post-optimization null model should be influenced by the question of interest, and it should also be constrained by properties of the network under examination. We explore such influences and constraints using our example networks.
1. Intra-layer and inter-layer null models
There are various ways to construct connectional null models (i.e., intra-layer null models), which randomize the connectivity structure within a network layer (Al).4,664. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 66. In the descriptions below, we use terms like “random” rewiring to refer to a process that we are applying uniformly at random aside from specified constraints. For binary networks, one can obtain ensembles of random graphs with the same mean degree as that of a real network using Erdős-Rényi random graphs,11. M. E. J. Newman, Networks: An Introduction (Oxford University Press, 2010). and ensembles of weighted random networks can similarly be constructed from weighted random graph models.6767. D. Garlaschelli, New J. Phys. 11, 073005 (2009). https://doi.org/10.1088/1367-2630/11/7/073005 To retain both the mean and distribution of edge weights, one can employ a permutation-based connectional null model that randomly rewires network edges with no additional constraints by reassigning uniformly at random the entire set of matrix elements Aijl in the lth layer (i.e., the matrix Al). Other viable connectional null models include ones that preserve degree21,6821. M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004). https://doi.org/10.1103/PhysRevE.69.026113 68. S. Maslov and K. Sneppen, Science 296, 910 (2002). https://doi.org/10.1126/science.1065103 or strength6969. G. Ansmann and K. Lehnertz, Phys. Rev. E 84, 026103 (2011). https://doi.org/10.1103/PhysRevE.84.026103 distributions, or—for networks based on time-series data—preserve length, frequency content, and amplitude distribution of the original time series.7070. S. Bialonski, M. Wendler, and K. Lehnertz, PLoS ONE 6, e22826 (2011). https://doi.org/10.1371/journal.pone.0022826 In this section, we present results for a few null models that are applicable to a variety of temporal networks. We note, however, that this is a fruitful area of further investigation.
We employ two connectional null models specific for the broad classes of networks represented by the brain and behavioral networks that we use as examples in this paper. The brain networks provide an example of time-dependent similarity networks, which are weighted and either fully connected or almost fully connected.3131. J.-P. Onnela, D. J. Fenn, S. Reid, M. A. Porter, P. J. Mucha, M. D. Fricker, and N. S. Jones, Phys. Rev. E. 86, 036104 (2012). https://doi.org/10.1103/PhysRevE.86.036104 (The brain networks have some 0 entries in their corresponding adjacency tensors because we have removed edges with weights that are not statistically significant.44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 ) We, therefore, employ a constrained null model that is constructed by randomly rewiring edges while maintaining the empirical degree distribution.6868. S. Maslov and K. Sneppen, Science 296, 910 (2002). https://doi.org/10.1126/science.1065103 In Fig. ???, we demonstrate the use of this null model to assess dynamic community structure. Importantly, this constrained null model can be used in principle for any binary or weighted network, though it does not take advantage of specific structure (aside from strength distribution) that one might want to exploit. For example, the behavioral networks have chain-like topologies, and it is desirable to develop models that are specifically appropriate for such situations. (One can obviously make the same argument for other specific topologies.) We, therefore, introduce a highly constrained connectional null model that is constructed by reassigning edge weights uniformly at random to existing edges. This does not change the underlying binary topology. (That is, we preserve network topology but scramble network geometry.) We demonstrate the use of this null model in Fig. ???.
In addition to intra-layer null models, one can also employ inter-layer null models—such as ones that scramble time or node identities.44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 For example, we construct a temporal null model by randomly permuting the order of the network layers. This temporal null model can be used to probe the existence of significant temporal evolution of community structure. One can also construct a nodal null model by randomly permuting the inter-layer edges that connect nodes in one layer to nodes in another. After the permutation is applied, an inter-layer edge can, for example, connect node i in layer t with node ji in layer t + 1 rather than being constrained to connect each node i in layer t with itself in layer t + 1. One can use this null model to probe the importance of node identity in network organization. We demonstrate the use of our temporal null model in row 2 of Fig. 6, and we demonstrate the use of our nodal null model in row 3 of Fig. 6.
2. Calculation of diagnostics on real versus null-model networks
We characterize the effects of post-optimization null models using four diagnostics: maximized modularity Q, the number of communities n, the mean community size s, and the stationarity ζ (see Sec. II B for definitions). Due to the possibly large number of partitions with nearly optimal Q,2424. B. H. Good, Y. A. de Montjoye, and A. Clauset, Phys. Rev. E 81, 046106 (2010). https://doi.org/10.1103/PhysRevE.81.046106 the values of such diagnostics vary over realizations of a computational heuristic for both the real and null-model networks. (We call this optimization variance.) The null-model networks also have a second source of variance (which we call randomization variance) from the myriad possible network configurations that can be constructed from a randomization procedure. We note that a third type of variance—ensemble variance—can also be present in systems containing multiple networks. In the example data sets that we discuss, this represents variability among experimental subjects.
We test for statistical differences between the real and null-model networks as follows. We first compute C = 100 optimizations of the modularity quality function for a network constructed from real data and then compute the mean of each of the four diagnostics over these C samples. This yields representative values of the diagnostics. We then maximize modularity for C different randomizations of a given null model (i.e., 1 optimization per randomization) and then compute the mean of each of the four diagnostics over these C samples. For both of our example data sets, we perform this two-step procedure for each network in the ensemble (60 brain networks and 66 behavioral networks; see “Sec. II”). We then investigate whether the set of representative diagnostics for the networks constructed from real data are different from those of appropriate ensembles of null-model networks. To address this issue, we subtract the diagnostic value for the null model from that of the real network for each subject and experimental session. We then use one-sample t-tests to determine whether the resulting distribution differs significantly from 0. We show our results in Fig. 6.
Results depend on all three factors (the data set, the null model, and the diagnostic), but there do seem to be some general patterns. For example, the real networks exhibit the most consistent differences from the nodal null model for all diagnostics and both data sets (see row 2 of Fig. 6). For both data sets, the variance of single-layer modularity in the real networks is consistently greater than those for all three null models, irrespective of the mean (see the final two columns of Figs. 6(a) and 6(b)); this is a potential indication of the statistical significance of the temporal evolution. However, although optimized modularity is higher in the real network for both data sets, the number of communities is higher in the set of brain networks and lower in the set of behavioral networks. Similarly, in comparison to the connectional null model, higher modularity is associated with a smaller mean community size in the brain networks but a larger mean size in the behavioral networks (see row 1 of Fig. 6). These results demonstrate that the three post-optimization null models provide different information about the network structure of the two systems and thereby underscores the critical need for further investigations of null-model construction.
C. Structural and temporal resolution parameters
When optimizing multilayer modularity, we must choose (or otherwise derive) values for the structural resolution parameter γ and the temporal resolution parameter ω. By varying γ, one can tune the size of communities within a given layer: large values of γ yield more communities, and small values yield fewer communities. A systematic method for how to determine values of ωjlr has not yet been discussed in the literature. In principle, one could choose different ωjlr values for different nodes, but we focus on the simplest scenario in which the value of ωjlrω is identical for all nodes j and all contiguous pairs of layers l and r (and is otherwise 0). In this framework, the temporal resolution parameter ω provides a means of tuning the number of communities discovered across layers: high values of ω yield fewer communities, and low values yield more communities. It is beneficial to study a range of parameter values to examine the breadth of structural (i.e., intra-layer24,25,7124. B. H. Good, Y. A. de Montjoye, and A. Clauset, Phys. Rev. E 81, 046106 (2010). https://doi.org/10.1103/PhysRevE.81.046106 25. S. Fortunato and M. Barthélemy, Proc. Natl. Acad. Sci. U.S.A. 104, 36 (2007). https://doi.org/10.1073/pnas.0605965104 71. V. A. Traag, P. Van Dooren, and Y. Nesterov, Phys. Rev. E 84, 016114 (2011). https://doi.org/10.1103/PhysRevE.84.016114 ) and temporal (i.e., inter-layer) resolutions of community structure, and some papers have begun to make progress in this direction.3,4,13,31,723. P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, Science 328, 876 (2010). https://doi.org/10.1126/science.1184819 4. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 13. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 31. J.-P. Onnela, D. J. Fenn, S. Reid, M. A. Porter, P. J. Mucha, M. D. Fricker, and N. S. Jones, Phys. Rev. E. 86, 036104 (2012). https://doi.org/10.1103/PhysRevE.86.036104 72. K. T. Macon, P. J. Mucha, and M. A. Porter, Physica A 391, 343 (2012). https://doi.org/10.1016/j.physa.2011.06.030
To characterize community structure as a function of resolution-parameter values (and hence of system scales), we quantify the quality of partitions using the mean value of optimized Q. To do this, we examine the constitution of the partitions using the mean similarity over C optimizations, and we compute partition similarities using the z-score of the Rand coefficient.7373. A. L. Traud, E. D. Kelsic, P. J. Mucha, and M. A. Porter, SIAM Rev. 53, 526 (2011). https://doi.org/10.1137/080734315 For comparing two partitions α and β, we calculate the Rand z-score in terms of the network's total number of pairs of nodes M, the number of pairs Mα that are in the same community in partition α, the number of pairs Mβ that are in the same community in partition β, and the number of pairs wαβ that are assigned to the same community both in partition α and in partition β. The z-score of the Rand coefficient comparing these two partitions is
zαβ=1σwαβ(wαβMαMβM),(14)
where σwαβ is the standard deviation of wαβ (as in Ref. 7373. A. L. Traud, E. D. Kelsic, P. J. Mucha, and M. A. Porter, SIAM Rev. 53, 526 (2011). https://doi.org/10.1137/080734315 ). Let the mean partition similarity z denote the mean value of zαβ over all possible partition pairs for αβ.
In Fig. 7, we show both z and optimized Q as a function of γ and ω in both brain and behavioral networks. The highest modularity values occur for low γ and high ω. The mean partition similarity is high for large γ in the brain networks, and it is high for both small and large γ in the behavioral networks. Interestingly, in both systems, the partition similarity when γ=ω=1 is lower than it is elsewhere in the (γ,ω) parameter plane, so the variability in partitions tends to be large at this point. Indeed, as shown in the second row of Fig. 7, modularity exhibits significant variability for γ=ω=1 compared to other resolution-parameter values.
It is useful to be able to determine the ranges of γ and ω that produce community structure that is significantly different from a particular null model. One can thereby use null models to probe resolution-parameter values at which a network displays interesting structures. This could be especially useful in systems for which one wishes to identify length scales (such as a characteristic mean community size) or time scales4,35,74,754. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 35. G. Palla, A. Barabási, and T. Vicsek, Nature 446, 664 (2007). https://doi.org/10.1038/nature05670 74. D. S. Bassett, E. T. Owens, K. E. Daniels, and M. A. Porter, Phys. Rev. E 86, 041306 (2012). https://doi.org/10.1103/PhysRevE.86.041306 75. X. Lou and J. A. Suykens, Chaos 21, 043116 (2011). https://doi.org/10.1063/1.3655371 directly from data.
In Fig. 8, we show examples of how the difference between diagnostic values for real and null-model networks varies as a function of γ and ω. As illustrated in panels (A) and (B), the brain and behavioral networks both exhibit a distinctly higher mean optimized modularity than the associated nodal null-model network for γω1. Interestingly, this roundly peaked difference in Q is not evident in comparisons of the real networks to temporal null-model networks (see Figs. 8(c) and 8(d)), so resolution-parameter values (and hence system scales) of potential interest might be more identifiable by comparison to nodal than to temporal null models in these examples. It is possible, however, that defining temporal layers over a longer or shorter duration would yield identifiable peaks in the difference in Q.
The differences in the Rand z-score landscapes are more difficult to interpret, as the values of mean partition similarity z are much larger in the real networks for some resolution-parameter values (positive differences; red) but are much larger in the null-model networks for other resolution-parameter values (negative differences; blue). The clearest situation occurs when comparing the brain's real and temporal null-model networks (see Fig. 8(c)), as the network built from real data exhibits a much larger value of z (and hence much more consistent optimization solutions) than the temporal null-model networks for high values of γ (i.e., when there are many communities) and low ω (i.e., when there is weak temporal coupling). These results are consistent with the fact that weak temporal coupling in a multilayer network facilitates greater temporal variability in network partitions across time. Such variability appears to be significantly different than the noise induced by scrambling time layers. These results suggest potential resolution values of interest for the brain system, as partitions are very consistent across many optimizations. For example, it would be interesting to investigate community structure in these networks for high γ (e.g., γ40) and low ω (e.g., ω0.1). At these resolution values, one can identify smaller communities with greater temporal variability than the communities identified for the case of γ=ω=1.44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108
The optimization and randomization variances appear to be similar in the brain and behavioral networks (see rows 2–3 in every panel of Fig. 8) not only in terms of their mean values but also in terms of their distribution in the part of the (γ,ω) parameter plane that we examined. In particular, the variance in Q is larger in the real networks precisely where the mean is also larger, so mean and variance are likely either dependent on one another or on some common source. Importantly, such dependence influences the ability to draw statistical conclusions because it is possible that the points in the (γ,ω) plane with the largest differences in mean are not necessarily the points with the most significant differences in mean.
We also find that the dependencies of the diagnostics on γ and ω are consistent across subjects and scans, suggesting that our results are ensemble-specific rather than individual-specific.
D. Examination of data generated from a dynamical system
Real-world data are often clouded by unknown or mathematically undefinable sources of variance, so it is also important to examine data sets generated from dynamical systems (or other models). Because we are concerned with time-dependent networks, we consider an example consisting of time-dependent data generated by a well-known dynamical system.
We construct a network of Kuramoto oscillators, in which the phase θi(t) of the ith oscillator evolves in time according to
dθidt=ωi+jκAijsin(θjθi),i{1,,N},(15)
where ωi is the natural frequency of oscillator i, the matrix A gives the binary coupling between each pair of oscillators, and κ is a positive real constant that indicates the strength of the coupling. We draw the frequencies ωi from a Gaussian distribution with mean 0 and standard deviation 1. In our simulations, we use a time step of τ=0.1, a constant of κ=0.2, and a network size of N = 128.
Kuramoto oscillators have been studied in the context of various network topologies and geometries51–53,76–7851. D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004). https://doi.org/10.1103/PhysRevLett.93.174102 52. Y. Kuramoto and D. Battogtokh, Nonlinear Phenom. Complex Syst. 5, 380 (2002). 53. S. I. Shima and Y. Kuramoto, Phys. Rev. E 69, 036213 (2004). https://doi.org/10.1103/PhysRevE.69.036213 76. A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93 (2008). https://doi.org/10.1016/j.physrep.2008.09.002 77. A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Phys. Rev. Lett. 96, 114102 (2006). https://doi.org/10.1103/PhysRevLett.96.114102 78. J. Stout, M. Whiteway, E. Ott, M. Girvan, and T. M. Antonsen, Chaos 21, 025109 (2011). https://doi.org/10.1063/1.3581168 and from both the component and ensemble perspectives.7979. G. Barlev, T. M. Antonsen, and E. Ott, Chaos 21, 025103 (2011). https://doi.org/10.1063/1.3596711 We are interested in networks with dynamic community structure. Following Refs. 77 and 8077. A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Phys. Rev. Lett. 96, 114102 (2006). https://doi.org/10.1103/PhysRevLett.96.114102 80. A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Physica D 224, 27 (2006). https://doi.org/10.1016/j.physd.2006.09.029 , we impose a well-defined community structure in which each community is composed of 16 nodes. In each time step, each node has 13 connections with nodes in its own community and 1 connection with nodes outside of its community (see Fig. 9(a)).
To quantify the temporal evolution of synchronization patterns, we define a set of temporal networks given by the time-dependent correlation between pairs of oscillators
ϕij(t)=|cos[θi(t)θj(t)]|,(16)
where the angular brackets indicate an average over 20 simulations. As time evolves from time step t = 0 to t = 100, oscillators tend to synchronize with other oscillators in their same community more quickly than with oscillators in other communities (see Fig. 9(b)).
To examine the performance of our multilayer community-detection techniques in this example, we compute Aijl=Aijt=ϕi,j(t) and using the multilayer extension of the Newman-Girvan null model Pijl given in Eq. (8). We separately optimize Q for two temporal regimes: (1) regime I (with t{1,,50}), for which the synchronization within communities increases rapidly; and (2) regime II (with t{51,,100}), for which the within-community synchronization level is roughly constant but the global synchronization still increases gradually. We set ω=1 and probe the effects of the structural resolution parameter γ in regime II. In Figs. 9(c) and 9(d), we illustrate that one can identify the value of γ that best uncovers the underlying hard-wired connectivity using troughs in the optimization variance of several diagnostics (e.g., maximized modularity, number of communities, and mean partition similarity).
We probe the community structure in regime I using the value of γ that best uncovered the underlying hard-wired connectivity in regime II. We observe temporal changes of community structure at early time points, as evidenced by the large number of communities for t{1,,5} (see Figs. 9(e) and 9(f)). Importantly, the temporal dependence of community number on t is not expected from a post-optimization temporal null model (see the right panel of Fig. 9(f)). We obtain qualitatively similar results when we optimize the multilayer modularity quality function over the entire temporal network without separating the data into two regimes.
Our results illustrate that one can use dynamic community detection to uncover the resolution of inherent hard-wired structure in a data set extracted from the temporal evolution of a dynamical system and that post-optimization null models can be used to identify regimes of unexpected temporal dependence in network structure.
E. Dealing with degeneracy: Constructing representative partitions
The multilayer modularity quality function has numerous near-degeneracies, so it is important to collect many instantiations when using a non-deterministic computational heuristic to optimize modularity.2424. B. H. Good, Y. A. de Montjoye, and A. Clauset, Phys. Rev. E 81, 046106 (2010). https://doi.org/10.1103/PhysRevE.81.046106 In doing this, an important issue is how (and whether) to distill a single representative partition from a (possibly very large) set of C partitions.8181. A. Lancichinetti and S. Fortunato, Sci. Rep. 2, 336 (2012). https://doi.org/10.1038/srep00336 In Fig. 10, we illustrate a new method for constructing a representative partition based on statistical testing in comparison to null models.
Consider C partitions from a single layer of an example multilayer brain network (see Fig. 10(a)). We construct a nodal association matrix T, where the element Tij indicates how many times nodes i and j have been assigned to the same community (see Fig. 10(b)). We then construct a null-model association matrix Tr based on random permutations of the original partitions (see Fig. 10(c)). That is, for each of the C partitions, we reassign nodes uniformly at random to the n communities of mean size s that are present in the selected partition. For every pair of nodes i and j, we let Trij be the number of times these two nodes have been assigned to the same community in this permuted situation (see Fig. 10(c)). The values Trij then form a distribution for the expected number of times two nodes are assigned to the same partition. Using an example with C = 100, we observe that two nodes can be assigned to the same community up to about 30 times out of the C partitions purely by chance. To be conservative, we remove such “noise” from the original nodal association matrix T by setting any element Tij whose value is less than the maximum entry of the random association matrix to 0 (see Fig. 10(d)). This yields the thresholded matrix T, which retains statistically significant relationships between nodes.
We use a Louvain-like algorithm to perform C optimizations of the single-layer modularity Q0 for the thresholded matrix T. Interestingly, this procedure typically extracts identical partitions for each of these optimizations in our examples (see Fig. 10(e)). This method, therefore, allows one to deal with the inherent near-degeneracy of the modularity quality function and provides a robust, representative partition of the original example brain network layer (see Fig. 10(f)).
We apply the same method to multilayer networks (see Fig. 11) to find a representative partition of (1) a real network over C optimizations, (2) a temporal null-model network over C randomizations, and (3) a nodal null-model network over C randomizations. Using these examples, we have successfully uncovered representative partitions when they appear to exist (e.g., in the real networks and the temporal null-model networks) and have not been able to uncover a representative partition when one does not appear to exist (e.g., in the nodal null-model network, for which each of the 112 brain regions is placed in its own community in the representative partition). We also note that the representative partitions in the temporal null-model and real networks largely match the original data in terms of both sizes and number of communities. These results indicate the potential of this method to uncover meaningful representative partitions over optimizations or randomizations in multilayer networks.
In this paper, we discussed methodological issues in the determination and interpretation of dynamic community structure in multilayer networks. We also analyzed the behavior of several null models used for optimizing quality functions (such as modularity) in multilayer networks.
We described the construction of networks and the effects that certain choices can have on the use of both optimization and post-optimization null models. We introduce novel modularity-optimization null models for two cases: (1) networks composed of ordered nodes (a “chain null model”) and (2) networks constructed from time-series similarities (FT and AAFT surrogates). We studied “connectional,” “temporal,” and “nodal” post-optimization null models using several multilayer diagnostics (optimized modularity, number of communities, mean community size, and stationarity) as well as novel single-layer diagnostics (in the form of measures based on time series for optimized modularity). To investigate the utility of such considerations for model-generated data, we also applied our methodology to time-series data generated from coupled Kuramoto oscillators.
We examined the dependence of optimized modularity and partition similarity on structural and temporal resolution parameters as well as the influence of their variances on putative statistical conclusions. Finally, we described a simple method to address the issue of near-degeneracy in the landscapes of modularity and similar quality functions using a method to construct a robust, representative partition from a network layer.
The present paper illustrates what we feel are important considerations in the detection of dynamic communities. As one considers data with increasingly complicated structures, network analyses of such data must become increasingly nuanced, and the purpose of this paper has been to discuss and provide some potential starting points for how to try to address some of these nuances.
We thank L. C. Bassett, M. Bazzi, K. Doron, L. Jeub, S. H. Lee, D. Malinow, and the anonymous referees for useful discussions. This work was supported by the Errett Fisher Foundation, the Templeton Foundation, David and Lucile Packard Foundation, PHS Grant NS44393, Sage Center for the Study of the Mind, and the Institute for Collaborative Biotechnologies through Contract No. W911NF-09-D-0001 from the U.S. Army Research Office. M.A.P. was supported by the James S. McDonnell Foundation (Research Award #220020177), the EPSRC (EP/J001759/1), and the FET-Proactive Project PLEXMATH (FP7-ICT-2011-8; Grant #317614) funded by the European Commission. He also acknowledges SAMSI and KITP for supporting his visits to them. P.J.M. acknowledges support from Award No. R21GM099493 from the National Institute of General Medical Sciences; the content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of General Medical Sciences or the National Institutes of Health.
  1. 1. M. E. J. Newman, Networks: An Introduction (Oxford University Press, 2010). Crossref
  2. 2. P. Holme and J. Saramäki, Phys. Rep. 519, 97 (2012). https://doi.org/10.1016/j.physrep.2012.03.001 , Crossref
  3. 3. P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, Science 328, 876 (2010). https://doi.org/10.1126/science.1184819 , Crossref, CAS
  4. 4. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 , Crossref, CAS
  5. 5. J.-P. Onnela, J. Saramäki, J. Hyvönen, G. Szabó, D. Lazer, K. Kaski, J. Kertész, and A. L. Barabási, Proc. Natl. Acad. Sci. U.S.A. 104, 7332 (2007). https://doi.org/10.1073/pnas.0610245104 , Crossref, CAS
  6. 6. Y. Wu, C. Zhou, J. Xiao, J. Kurths, and H. J. Schellnhuber, Proc. Natl. Acad. Sci. U.S.A. 107, 18803 (2010). https://doi.org/10.1073/pnas.1013140107 , Crossref, CAS
  7. 7. S. Gonzaález-Bailón, J. Borge-Holthoefer, A. Rivero, and Y. Moreno, Sci. Rep. 1, 197 (2011). https://doi.org/10.1038/srep00197 , Crossref
  8. 8. T. J. Fararo and J. Skvoretz, in Status, Network, and Structure: Theory Development in Group Processes (Stanford University Press, 1997), pp. 362–386.
  9. 9. A. Stomakhin, M. B. Short, and A. L. Bertozzi, Inverse Probl. 27, 115013 (2011). https://doi.org/10.1088/0266-5611/27/11/115013 , Crossref
  10. 10. P. J. Mucha and M. A. Porter, Chaos 20, 041108 (2010). https://doi.org/10.1063/1.3518696 , Scitation
  11. 11. A. S. Waugh, L. Pei, J. H. Fowler, P. J. Mucha, and M. A. Porter, “ Party polarization in Congress: A network science approach,” arXiv:0907.3509 (2011).
  12. 12. K. W. Doron, D. S. Bassett, and M. S. Gazzaniga, Proc. Natl. Acad. Sci. U.S.A. 109, 18661 (2012). https://doi.org/10.1073/pnas.1216402109 , Crossref, CAS
  13. 13. N. F. Wymbs, D. S. Bassett, P. J. Mucha, M. A. Porter, and S. T. Grafton, Neuron 74, 936 (2012). https://doi.org/10.1016/j.neuron.2012.03.038 , Crossref, CAS
  14. 14. D. J. Fenn, M. A. Porter, M. McDonald, S. Williams, N. F. Johnson, and N. S. Jones, Chaos 19, 033119 (2009). https://doi.org/10.1063/1.3184538 , Scitation
  15. 15. D. J. Fenn, M. A. Porter, S. Williams, M. McDonald, N. F. Johnson, and N. S. Jones, Phys. Rev. E 84, 026109 (2011). https://doi.org/10.1103/PhysRevE.84.026109 , Crossref
  16. 16. M. A. Porter, J.-P. Onnela, and P. J. Mucha, Not. Am. Math. Soc. 56, 1082 (2009).
  17. 17. H. Simon, Proc. Am. Philos. Soc. 106, 467 (1962).
  18. 18. M. A. Porter, P. J. Mucha, M. Newman, and C. M. Warmbrand, Proc. Natl. Acad. Sci. U.S.A. 102, 7057 (2005). https://doi.org/10.1073/pnas.0500191102 , Crossref, CAS
  19. 19. A. C. F. Lewis, N. S. Jones, M. A. Porter, and C. M. Deane, BMC Syst. Biol. 4, 100 (2010). https://doi.org/10.1186/1752-0509-4-100 , Crossref
  20. 20. S. Fortunato, Phys. Rep. 486, 75 (2010). https://doi.org/10.1016/j.physrep.2009.11.002 , Crossref
  21. 21. M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004). https://doi.org/10.1103/PhysRevE.69.026113 , Crossref, CAS
  22. 22. A. Lancichinetti, F. Radicchi, and J. J. Ramasco, Phys. Rev. E 81, 046110 (2010). https://doi.org/10.1103/PhysRevE.81.046110 , Crossref
  23. 23. M. E. J. Newman, Nat. Phys. 8, 25 (2012). https://doi.org/10.1038/nphys2162 , Crossref, CAS
  24. 24. B. H. Good, Y. A. de Montjoye, and A. Clauset, Phys. Rev. E 81, 046106 (2010). https://doi.org/10.1103/PhysRevE.81.046106 , Crossref
  25. 25. S. Fortunato and M. Barthélemy, Proc. Natl. Acad. Sci. U.S.A. 104, 36 (2007). https://doi.org/10.1073/pnas.0605965104 , Crossref, CAS
  26. 26. P. J. Bickel and A. Chen, Proc. Natl. Acad. Sci. U.S.A. 106, 21068 (2009). https://doi.org/10.1073/pnas.0907096106 , Crossref, CAS
  27. 27. S. Hutchings, “ The behavior of modularity-optimizing community detection algorithms,” M.Sc. Thesis (University of Oxford, 2011).
  28. 28. M. E. J. Newman, Phys. Rev. E 69, 066133 (2004). https://doi.org/10.1103/PhysRevE.69.066133 , Crossref, CAS
  29. 29. M. E. J. Newman, Proc. Natl. Acad. Sci. U.S.A. 103, 8577 (2006). https://doi.org/10.1073/pnas.0601602103 , Crossref, CAS
  30. 30. J. Reichardt and S. Bornholdt, Phys. Rev. E 74, 016110 (2006). https://doi.org/10.1103/PhysRevE.74.016110 , Crossref
  31. 31. J.-P. Onnela, D. J. Fenn, S. Reid, M. A. Porter, P. J. Mucha, M. D. Fricker, and N. S. Jones, Phys. Rev. E. 86, 036104 (2012). https://doi.org/10.1103/PhysRevE.86.036104 , Crossref
  32. 32. U. Brandes, D. Delling, M. Gaertler, R. Görke, M. Hoefer, Z. Nikoloski, and D. Wagner, IEEE Trans. Knowl. Data Eng. 20, 172 (2008). https://doi.org/10.1109/TKDE.2007.190689 , Crossref
  33. 33. M. E. J. Newman, Phys. Rev. E 74, 036104 (2006). https://doi.org/10.1103/PhysRevE.74.036104 , Crossref, CAS
  34. 34. T. Richardson, P. J. Mucha, and M. A. Porter, Phys. Rev. E 80, 036111 (2009). https://doi.org/10.1103/PhysRevE.80.036111 , Crossref
  35. 35. G. Palla, A. Barabási, and T. Vicsek, Nature 446, 664 (2007). https://doi.org/10.1038/nature05670 , Crossref, CAS
  36. 36. Equation (4) gives the definition for the notion of stationarity that we used in Ref. 44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 . The equation for this quantity in Ref. 44. D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, Proc. Natl. Acad. Sci. U.S.A. 108, 7641 (2011). https://doi.org/10.1073/pnas.1018985108 has a typo in the denominator. We wrote incorrectly in that paper that the denominator was tt01, whereas the numerical computations in that paper used tt0.
  37. 37. In other areas of investigation, it probably should be.
  38. 38. E. T. Bullmore and D. S. Bassett, Ann. Rev. Clin. Psych. 7, 113 (2011). https://doi.org/10.1146/annurev-clinpsy-040510-143934 , Crossref
  39. 39. D. S. Bassett, J. A. Brown, V. Deshpande, J. M. Carlson, and S. T. Grafton, Neuroimage 54, 1262 (2011). https://doi.org/10.1016/j.neuroimage.2010.09.006 , Crossref
  40. 40. A. Zalesky, A. Fornito, I. H. Harding, L. Cocchi, M. Yucel, C. Pantelis, and E. T. Bullmore, Neuroimage 50, 970 (2010). https://doi.org/10.1016/j.neuroimage.2009.12.027 , Crossref
  41. 41. J. Wang, L. Wang, Y. Zang, H. Yang, H. Tang, Q. Gong, Z. Chen, C. Zhu, and Y. He, Hum. Brain Mapp 30, 1511 (2009). https://doi.org/10.1002/hbm.20623 , Crossref
  42. 42. S. Bialonski, M. T. Horstmann, and K. Lehnertz, Chaos 20, 013134 (2010). https://doi.org/10.1063/1.3360561 , Scitation
  43. 43. A. A. Ioannides, Curr. Opin. Neurobiol. 17, 161 (2007). https://doi.org/10.1016/j.conb.2007.03.008 , Crossref, CAS
  44. 44. C. T. Butts, Science 325, 414 (2009). https://doi.org/10.1126/science.1171022 , Crossref, CAS
  45. 45. G. Kossinets, Soc. Networks 28, 247 (2006). https://doi.org/10.1016/j.socnet.2005.07.002 , Crossref
  46. 46. A. Clauset, C. Moore, and M. E. J. Newman, Nature 453, 98 (2008). https://doi.org/10.1038/nature06830 , Crossref, CAS
  47. 47. R. Guimerá and M. Sales-Pardo, Proc. Natl. Acad. Sci. U.S.A. 106, 22073 (2009). https://doi.org/10.1073/pnas.0908366106 , Crossref, CAS
  48. 48. M. Kim and J. Leskovec, SIAM International Conference on Data Mining (2011).
  49. 49. O. Sporns, Networks of the Brain (MIT, 2010).
  50. 50. E. Bullmore and O. Sporns, Nat. Rev. Neurosci. 10, 186 (2009). https://doi.org/10.1038/nrn2575 , Crossref, CAS
  51. 51. D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004). https://doi.org/10.1103/PhysRevLett.93.174102 , Crossref
  52. 52. Y. Kuramoto and D. Battogtokh, Nonlinear Phenom. Complex Syst. 5, 380 (2002).
  53. 53. S. I. Shima and Y. Kuramoto, Phys. Rev. E 69, 036213 (2004). https://doi.org/10.1103/PhysRevE.69.036213 , Crossref
  54. 54. K. J. Friston, Hum. Brain Mapp 2, 56 (1994). https://doi.org/10.1002/hbm.460020107 , Crossref
  55. 55. Y. Benjamini and Y. Yekutieli, Ann. Stat. 29, 1165 (2001). https://doi.org/10.1214/aos/1013699998 , Crossref
  56. 56. E. Bullmore and O. Sporns, Nat. Rev. Neurosci. 13, 336 (2012). https://doi.org/10.1038/nrn3214 , Crossref, CAS
  57. 57. V. D. Blondel, J. L. Guillaume, R. Lambiotte, and E. Lefebvre, J. Stat. Mech. 10, P10008 (2008). Crossref
  58. 58. I. S. Jutla, L. G. S. Jeub, and P. J. Mucha, “ A generalized Louvain method for community detection implemented in matlab,” (2011–2012); http://netwiki.amath.unc.edu/GenLouvain.
  59. 59. A discrete time series can be represented as a vector. A continuous time series would first need to be discretized.
  60. 60. D. Prichard and J. Theiler, Phys. Rev. Lett. 73, 951 (1994). https://doi.org/10.1103/PhysRevLett.73.951 , Crossref, CAS
  61. 61. The code used for this computation actually operates on [0,2π]. However, this should be an equivalent mathematical estimate to the same computation on [0,2π), which is the same except for a set of measure zero.
  62. 62. H. Nakatani, I. Khalilov, P. Gong, and C. van Leeuwen, Phys. Lett. A 319, 167 (2003). https://doi.org/10.1016/j.physleta.2003.09.082 , Crossref, CAS
  63. 63. A. Zalesky, A. Fornito, and E. Bullmore, Neuroimage 60, 2096 (2012). https://doi.org/10.1016/j.neuroimage.2012.02.001 , Crossref
  64. 64. J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, and J. D. Farmer, Physica D 58, 77 (1992). https://doi.org/10.1016/0167-2789(92)90102-S , Crossref
  65. 65. It is also important to note that the AAFT method can allow nonlinear correlations to remain in the surrogate data. Therefore, the development of alternative surrogate data generation methods might be necessary (Refs. 83 and 84).
  66. 66. In the descriptions below, we use terms like “random” rewiring to refer to a process that we are applying uniformly at random aside from specified constraints.
  67. 67. D. Garlaschelli, New J. Phys. 11, 073005 (2009). https://doi.org/10.1088/1367-2630/11/7/073005 , Crossref
  68. 68. S. Maslov and K. Sneppen, Science 296, 910 (2002). https://doi.org/10.1126/science.1065103 , Crossref, CAS
  69. 69. G. Ansmann and K. Lehnertz, Phys. Rev. E 84, 026103 (2011). https://doi.org/10.1103/PhysRevE.84.026103 , Crossref
  70. 70. S. Bialonski, M. Wendler, and K. Lehnertz, PLoS ONE 6, e22826 (2011). https://doi.org/10.1371/journal.pone.0022826 , Crossref, CAS
  71. 71. V. A. Traag, P. Van Dooren, and Y. Nesterov, Phys. Rev. E 84, 016114 (2011). https://doi.org/10.1103/PhysRevE.84.016114 , Crossref, CAS
  72. 72. K. T. Macon, P. J. Mucha, and M. A. Porter, Physica A 391, 343 (2012). https://doi.org/10.1016/j.physa.2011.06.030 , Crossref
  73. 73. A. L. Traud, E. D. Kelsic, P. J. Mucha, and M. A. Porter, SIAM Rev. 53, 526 (2011). https://doi.org/10.1137/080734315 , Crossref
  74. 74. D. S. Bassett, E. T. Owens, K. E. Daniels, and M. A. Porter, Phys. Rev. E 86, 041306 (2012). https://doi.org/10.1103/PhysRevE.86.041306 , Crossref
  75. 75. X. Lou and J. A. Suykens, Chaos 21, 043116 (2011). https://doi.org/10.1063/1.3655371 , Scitation
  76. 76. A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93 (2008). https://doi.org/10.1016/j.physrep.2008.09.002 , Crossref
  77. 77. A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Phys. Rev. Lett. 96, 114102 (2006). https://doi.org/10.1103/PhysRevLett.96.114102 , Crossref
  78. 78. J. Stout, M. Whiteway, E. Ott, M. Girvan, and T. M. Antonsen, Chaos 21, 025109 (2011). https://doi.org/10.1063/1.3581168 , Scitation
  79. 79. G. Barlev, T. M. Antonsen, and E. Ott, Chaos 21, 025103 (2011). https://doi.org/10.1063/1.3596711 , Scitation
  80. 80. A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente, Physica D 224, 27 (2006). https://doi.org/10.1016/j.physd.2006.09.029 , Crossref
  81. 81. A. Lancichinetti and S. Fortunato, Sci. Rep. 2, 336 (2012). https://doi.org/10.1038/srep00336 , Crossref
  82. 82. A. L. Traud, C. Frost, P. J. Mucha, and M. A. Porter, Chaos 19, 041104 (2009). https://doi.org/10.1063/1.3194108 , Scitation
  83. 83. C. Räth, M. Gliozzi, I. E. Papadakis, and W. Brinkmann, Phys. Rev. Lett. 109, 144101 (2012). https://doi.org/10.1103/PhysRevLett.109.144101 , Crossref, CAS
  84. 84. G. Rossmanith, H. Modest, C. Räth, A. J. Banday, K. M. Gorski, and G. Morfill, Phys. Rev. D 86, 083005 (2012). https://doi.org/10.1103/PhysRevD.86.083005 , Crossref
  85. © 2013 American Institute of Physics.