Mapping Multiplex Hubs in Human Functional Brain Networks

Typical brain networks consist of many peripheral regions and a few highly central ones, i.e., hubs, playing key functional roles in cerebral inter-regional interactions. Studies have shown that networks, obtained from the analysis of specific frequency components of brain activity, present peculiar architectures with unique profiles of region centrality. However, the identification of hubs in networks built from different frequency bands simultaneously is still a challenging problem, remaining largely unexplored. Here we identify each frequency component with one layer of a multiplex network and face this challenge by exploiting the recent advances in the analysis of multiplex topologies. First, we show that each frequency band carries unique topological information, fundamental to accurately model brain functional networks. We then demonstrate that hubs in the multiplex network, in general different from those ones obtained after discarding or aggregating the measured signals as usual, provide a more accurate map of brain's most important functional regions, allowing to distinguish between healthy and schizophrenic populations better than conventional network approaches.


INTRODUCTION
The brain functional network is generally built by interconnecting brain regions according to some measure of functional connectivity Sporns, 2009, 2012). Studies using functional magnetic resonance imaging (Van Den Heuvel and Pol, 2010;Poldrack and Farah, 2015) (fMRI) provided convincing evidence supporting the existence of special regions, i.e., hubs, that play a fundamental role in brain functional connectivity Power et al., 2013) by mediating interactions among other regions and favoring the brain's integrated operation. Generally, the strength of this connectivity is empirically estimated by inter-regional correlations calculated after post-processing and filtering fMRI signals with a conventional pass band, keeping components between 0.01 and 0.1 Hz (Cordes et al., 2001(Cordes et al., , 2002Fox and Raichle, 2007). The importance of each region with respect to the overall connectivity, i.e., nodal centrality in the functional network, is of particular interest in many applications (Sporns et al., 2007;Bassett et al., 2008;Bullmore and Sporns, 2009;Lynall et al., 2010;Rubinov and Sporns, 2010;Zuo et al., 2012). However, it has been shown that networks with unique hub regions can be built from different frequency ranges (Sasai et al., 2014) and that region centrality might largely fluctuate depending on frequency cuts (Thompson and Fransson, 2015), with components above 0.1 Hz also contributing to functional connectivity with unique topological information Mantini et al., 2007;Supekar et al., 2008;Chavez et al., 2010;Liao et al., 2013;Chen and Glover, 2015). Such an evidence impels the development of a novel framework to account for full information from all frequency bands separately and simultaneously, without discarding any particular component or aggregating some of them to build single networks.
In this study, we tackle this challenging issue by employing the theoretical and computational tools recently developed for analyzing and modeling multiplex networks (Mucha et al., 2010;De Domenico et al., 2013, 2015b. Multiplex architectures are special networks consisting of different layers, each encoding a different type of relationship or interaction between nodes (Boccaletti et al., 2014;Kivelä et al., 2014). Recent studies modeled and analyzed brain networks using temporal networks, a special type of multilayer system (Bassett et al., 2011(Bassett et al., , 2013Braun et al., 2015). In this context, we identify each frequency component with a distinct layer of a multiplex network whose nodes represent the brain's regions of interest and edges represent their functional connectivity in a specific frequency range.
This novel approach rises two fundamental questions, requiring to (i) verify if and how brain regions playing the role of hubs in the new multiplex functional network differ from the ones obtained using standard network approaches; and (ii) if and how we can exploit such differences to improve our understanding of brain disorders. In the following, we will provide extensive evidence demonstrating that hub regions in multiplex functional networks are different from hub regions in standard functional networks and that such differences in the nodal centrality profile allow us to identify patients affected by schizophrenia more accurately than conventional approaches based on discarding or aggregating information about brain functional activity.

Overview of the Data Set and fMRI Preprocessing
The publicly available MR data set contributed by The Center for Biomedical Research Excellence (COBRE) was used in this study. The data set was downloaded from the following repository: http://fcon_1000.projects.nitrc.org/indi/retro/cobre. html. It includes resting functional and anatomical MRI data acquired from 71 Schizophrenic patients and 74 healthy controls (age: 18-65 for both groups). Parameters of fMRI acquisition released by the provider are as follows: TR = 2000 ms, TE = 29 ms, voxel size = 3 × 3 × 4 mm 3 , total scan time = 6 min. One patients data was discarded from all analyses due to the shortness of the data length. The following pre-processing steps were applied to functional MR images by using the SPM8 package (Wellcome Department of Imaging Neuroscience, London, UK): motion-correction, slice-timing-correction, spatial smoothing with Gaussian kernel (5-mm full-width-at-half-maximum) and spatial normalization. Signal fluctuations of fMRI are driven by not only neural but also physiological effects-such as respiration and cardiac pulsation-and environmental conditions-such as scanner instabilities and subject motion. These nuisance effects can be canceled out by discarding, for instance, the signal from the ROI centered in the white matter, the signal from the ventricular ROI, and the signal from the ROI located within the soft-tissue. We have linearly removed these components as well as six motion-correction parameters after temporally shifting them by optimal time-lags yielding the highest correlation with the averaged signal of all gray matter voxels (Anderson et al., 2011).

Statistical Analysis of Age, Gender, and Handed-Ness Distribution
We performed a Wilcoxon rank-sum test, a non-parametric version of unpaired two-sample t-test, to test the null hypothesis that phenotypic details in the two groups are sampled from continuous distributions with equal medians. The test did not reject the null in the case of age (p-value = 0.4253) and gender (p-value = 0.1186), therefore the discrimination power of the method proposed in the following can not be explained by differences in age or gender between the two groups (see Supplementary Table 3 for statistical descriptors). The test rejected the null in the case of handed-ness (p-value = 0.004), suggesting that this factor might affect the discrimination power of our method. However, we show in Supplementary  Figure 3 that, by including such information in the classification procedure, discrimination accuracy and all other statistical indicators are not significantly improved with respect to the case when phenotypic data is not accounted for, suggesting that differences in handed-ness are not responsible for our findings.

Connectivity Matrices
A set of 264 spherical ROIs (5 mm radii) was used to extract the mean signal within each ROI. For each individual, the coherence between all pairs of in-ROI averaged signals was estimated in specific frequency bands, as described in the text. We kept the edges between pairs of ROIs whose weight was significantly different from a null model where observed signals were replaced by surrogates. More specifically, we used the well-known iterative amplitude-adjusted Fourier transform (IAAFT) algorithm to build surrogate time series preserving the power spectrum and the probability density of the original ones, while removing higher-order self-correlations. For each pair of ROIs i and j, we have verified that the distribution of the weights obtained from the null model corresponds to a Gaussian described by sample mean µ ij and variance σ 2 ij . Let w ij indicate the weight obtained from empirical data: we have calculated the absolute Zscore as z ij = |w ij − µ ij |/σ ij and discarded all those edges for which z ij < 3, corresponding to cross-coherence not statistically significant. It is worth remarking that the chosen threshold provides a statistical test with significance 99.8%. On one hand, lower thresholds would keep links that are more likely to be observed by chance; on the other hand, higher thresholds would dramatically reduce the density of the network making any further analysis less reliable. Our choice provides a good tradeoff between these two extremal cases. Finally, we used the values z ij as entries of the resulting connectivity matrix.
As a final remark, it is worth mentioning that using the absolute value of z-scores does not allow to distinguish between significantly correlated and anti-correlated signals, a characteristic that is typically exploited in the neuroscience literature. In a future study, we plan to take into account, separately, the information obtained from correlated and anticorrelated signals by using distinct layers. By comparing against other standard metrics, we will be able to quantify how our results might be further improved.

Multiplex Network Model
A multilayer network allows to encode different types of interactions or relationships among a set of nodes. More specifically, in the case of our study we make use of multiplex networks to model functional connectivity. In a multiplex network, the links are of different type: one can assign a different "color" to each type, thus obtained an edge-colored representation of the network. In this type of architectures, nodes exist in one or more layers, i.e., it is not required that all nodes exist in all layers. Correlation networks, as the ones used in this study, define edge-colored graphs where each layer encodes the correlations observed in a specific frequency band. However, it has been shown that by interconnecting nodes with their replicas across layers, the resulting interconnected multiplex network can be described by an adjacency tensor (De Domenico et al., 2013) with components M iα jβ , an object generalizing the well-known concept of adjacency matrix to higher orders, encoding connections between node i in layer α and node j in layer β. While a strict biological interpretation of inter-layer links in this context might be difficult, we suggest that their existence and their weight might encode how different frequency components interact each other and with which intensity. For interconnected multiplex networks, M iα jβ = 0 for i =j and simultaneously α = β. The presence of interconnections allows to exploit tensorial algebra to generalize many singlelayer network descriptors, from centrality (De Domenico et al., 2015c) to mesoscale structure (Mucha et al., 2010;De Domenico et al., 2015a). However, it is not always possible to assign a weight to inter-layer links by using the data, and it is common to parameterize the intensity of interconnections (Gomez et al., 2013;De Domenico et al., 2014), i.e., M iα jβ = D for i = j and simultaneously α = β, to study the resulting interconnected multiplex network as a function of this parameter D. This is exactly the case of the present study, where the choice of D depends on the analysis of interest.
It is worth mentioning that, in general, care must be taken when network measures are applied to functional networks, because links between nodes do not directly map physical connections between different areas of the brain. Nevertheless, links in functional networks indirectly quantify the statistical correlation or similarity between two regions of interest and are widely used in literature to gain insight about brain's activity (Bullmore and Sporns, 2009). In the following, we will not make the difference between positive and negative correlations, that will be considered, for each band frequency, as additional layers of the multiplex functional brain network in a future study.

Structural Reducibility of Brain Multiplex Functional Network
The analysis of structural reducibility of a multilayer network allows to find layers that provide redundant topological information, suggesting how to merge some layers with other ones, to obtain an optimal multilayer network (De Domenico et al., 2015b). The whole procedure can be summarized as follows: (i) compute the distance (based on quantum Jensen-Shannon divergence) between all pairs of layers; (ii) perform hierarchical clustering of layers using such distance matrix and use changes in the relative entropy q(•) as the quality function for the resulting partition; (iii) finally, choose the partition which maximizes the quality function, i.e., the distinguishability from the fully aggregated graph obtained by summing up the adjacency matrices of all layers. It is worth remarking that this analysis is independent on the choice of interconnections weight, i.e., it does not depend on D. Here, we do not enter into the details of the whole method; instead we focus on the Jensen-Shannon distance, that is a key measure for two analyses presented in this study.
The components A [α] ij (i, j = 1, 2, ..., N; being N the number of ROIs in this study) of the adjacency matrix A [α] -encoding layer α-are obtained from the components of the multilayer adjacency tensor as ij > 0 if there is correlation between ROIs i and j in the frequency band represented by α. The Von Neumann entropy (Braunstein et al., 2006;Passerini and Severini, 2010) of the corresponding complex network is defined by where ij , and S is the diagonal matrix of the strengths of the nodes. From the eigen-decomposition of the Laplacian, it is possible to show that the entropy can be calculated by where . The similarity of two layers can be calculated in terms of differences in their entropy. Given two rescaled Laplacian matrices L [α] and L [β] , it is possible to quantify to which extent layer α is different from layer β by their Kullback-Liebler divergence, defined by encoding the information gained about L [β] when the expectation is based only on L [α] . This divergence is not a metric and a more suitable dissimilarity measure is the Jensen-Shannon divergence, defined by where . It can be shown that D JS , usually called Jensen-Shannon distance, takes values in [0, 1], satisfies all the properties of a metric distance and provides a very powerful measure of dissimilarity between layers.

Random Forest Classification
Machine learning has been used to train a classifier to distinguish between control and schizophrenic individuals. We used the random forest classifier (Breiman, 2001), well-known for its robustness and for facilitating the interpretation of results. We have fixed to 5 the maximum number of terminal nodes trees the forest can have and to 2 the number of variables randomly sampled as candidates at each split. We have verified that forests consisting of 700 trees where enough to reach stable results within this setup.
Given the importance of interconnections weight for our analysis and, at the same time, the lack of knowledge about its value, we used random forest to learn also which value of D would be more suitable for calculations.
We have performed a first exploratory classification using a leave-one-out approach to maximize the amount of data used for training the classifier. The result of each classification, corresponding to exactly one different individual (without replacement) left out, was accompanied by the importance assigned by the classifier to each ROI in terms of mean decrease in its Gini index. Therefore, for each individual and each value of D, we have ranked the ROIs according to this measure and, eventually, summed up the ranks corresponding to all classifications.
The result of the exploratory classification was an overall ranking suggesting which ROIs, in general, have been more crucial than others in the classification process. Therefore, we performed a second classification round by using only the top ROIs according to the above ranking. We first varied the number of kept features and the value of D, to find the values with best classification performances in terms of accuracy (see Supplementary Figure 1). The numerical analysis indicated that the best classification is achieved for interconnections weight close to 24.7708 and about 30 top ROIs: that value of D and that sub-set of ROIs have been used for analysis reported in the text.
Using a similar approach, we have compared the best performance obtained from the full multiplex functional network (12 layers) against multiplex functional networks obtained by keeping layers in the typical band (Supplementary Figure 2) and against classifier trained by including phenotypic data (Supplementary Figure 3). In all cases, the classification obtained using the full multiplex functional network was equal or better than the other ones.
FIGURE 1 | Schematic illustration of brain multiplex functional network construction. (A) We measure the brain activity with a set of 264 ROIs (here, we only draw five ROIs, for simplicity), and estimate the coherence spectrum of signals between any pair of ROIs. (B) Averaged coherence values are calculated in 12 frequency bands (here we only show four bands, for simplicity), to quantify the strength of frequency-specific functional connectivity. The statistical significance of each connection is calculated (see Methods) and connections with Z-score smaller than 3 are discarded. (C) The remaining connections are used to build adjacency matrices, weighted by Z-scores, that constitute the layers of the multiplex functional network once interconnected. (D) Resulting single-layer and multiplex networks obtained from this procedure.

ROIs PageRank Centrality
PageRank is a measure of node's centrality originally introduced by Google founders to rank Web pages according to their importance in the Internet (Brin and Page, 1998;Ermann et al., 2015). The algorithm consider a random walker exploring the network with the following rules: 85% of times the walker jumps from the current node to one chosen with uniform probability from the neighborhood, whereas 15% of times the walker is allowed to jump to any node of the network, with uniform probability. The stationary probability of finding the walker in a specific node is then used to rank the importance of nodes in the network, the rationale being that central nodes have high number of incoming links from other important nodes.
The natural extension of the PageRank algorithm to the context of multiplex networks has been recently introduced (De Domenico et al., 2015c) and proved to perform better than its single-layer counterpart in some applications. Let us indicate with R iα jβ the transition tensor, governing the dynamics of a random walker jumping to neighboring nodes with rate 0.85 and teleporting to any other node in the network with rate 0.15. This rank-4 tensor is given by where T iα jβ governs the standard moves of a classical random walker from a node i in layer α to one of its neighbors j in layer β, L is the total number of layers and u iα jβ is the rank-4 tensor with all components equal to 1. The steady-state solution of the master equation obtained in the limit t −→ ∞, provides the PageRank centrality for interconnected multiplex networks. To compute the overall PageRank of a node, accounting for the whole interconnected topology, we can safely sum up the stationary probabilities π ⋆ jβ over the layers, to obtain the components of the centrality profile vectorπ ⋆ j = L β=1 π ⋆ jβ used in our analysis. It is worth remarking that the interconnection weight used for this purpose is D = 24.7708, the one yielding the highest classification accuracy.

Building Aggregated and Multiplex Functional Connectivity Networks
We use a publicly available COBRE data set of resting state fMRI, consisting of 71 patients affected by Schizophrenia and 74 healthy controls (age: 18-65). The set of 264 regions of interest (ROIs) introduced by Power et al. (2011) is used to extract the mean signal within each ROI, for each individual separately. After estimating coherence between all pairs of ROIs, the frequency-specific connectivity matrices are obtained by averaging coherence within 12 frequency bands, defined by decomposing the frequency range from 0.01 to 0.25 Hz into intervals with equal widths of 0.02 Hz. The upper bound of this frequency range corresponds to the Nyquist frequency of fMRI signals, while the lower bound is obtained by following conventional way to eliminate long term drift (Cordes et al., 2002). Weighted adjacency matrices, defining the functional network for each frequency component separately, are yielded by discarding from frequency-specific connectivity matrices those connections with non-significant amount of correlation (see Materials and Methods). The resulting multiplex network is obtained, for each individual separately in control and patient groups, by interconnecting the layers encoding functional connectivity in each frequency band (Figures 1A-C). We also define two single-layer networks, obtained by averaging coherence signals within 0.01-0.25 Hz and 0.01-0.1 Hz frequency ranges ( Figure 1D). We refer to such conventional networks as full-band and typical-band single-layer networks, respectively, both representing averaged and filtered versions of the full multiplex functional networks.

Frequency-Dependent Structural Analysis and Small-Worldness
For each subject and for each layer of the corresponding multiplex brain network, we calculate some well-known structural descriptors to better characterize the networks.
We show in Figure 2 the distribution of average degree and average strength-characterizing ROI's mean number and weight, respectively, of functional connections in each frequency band-and ROI's assortative mixing (Newman, 2002)-characterizing the tendency of ROIs to connect to other ROIs with similar or dissimilar connectivity. For both healthy and schizophrenic brains we observe similar distributions and FIGURE 3 | Frequency-dependent clustering. As in Figure 2, showing the distribution of edge density (left panels), average weighted clustering coefficient (central panels) and modularity (right panels). positive assortativity, meaning that, on average, nodes with similar degree tend to be functionally connected each other. Figure 3 shows the distribution of edge density-defined by the ratio between the number of links in the network and the maximum number of possible connections-, average weighted clustering coefficient (Barrat et al., 2004)characterizing the tendency of nodes to form weighted triads-and modularity (Blondel et al., 2008)-quantifying the mesoscale organization of nodes into functional modules. We find that, in both groups, layers are moderately dense and modular, with highly clustered ROIs.
Finally, we are interested in quantifying to which extent the characteristic path length of each layer (Figure 4, left panel)defined by the average over all shortest paths connecting ROIs each other-and the average local clustering (Figure 4, middle panel) are different from their expectation when using an Erdos-Renyi graph, with the same number of nodes and edges, to model the network topology. We find that regardless of the frequency band and the group, the functional brain networks are characterized by short average path length and high clustering, property typical of small-world topologies (Watts and Strogatz, 1998). To better quantify this finding, we calculate the small-world index (Humphries et al., 2006;Sporns, 2011) defined by S(f ) = (C(f )/C rand (f ))/(l(f )/ℓ rand (f )), whereC(f ) is the average local clustering coefficient andl(f ) is the average path length, with C rand (f ) and ℓ rand (f ) being their random expectations. These quantities are calculated for each layer separately and depend on the corresponding frequency band. The results, shown in the right panel of Figure 4 confirm that all layers are, on average, characterized by a small-world topology, being the median small-world index close to 2.5.

Structural Reducibility of the Multiplex Functional Connectivity Network
First, we verify if the multiplex network is a valid and suitable model of the underlying brain connectivity. For this purpose, we analyze the structural reducibility of a multiplex network (De Domenico et al., 2015b), allowing to identify layers carrying redundant topological information. The method incorporates redundant layers into other ones to reduce the overall structure, while still maximizing the distinguishability between the multiplex network model and the corresponding fully aggregated graph, obtained by summing up the connectivity FIGURE 4 | Frequency-dependent small-worldness. As in Figure 2, showing the distribution of average path length (left panels), average local clustering coefficient (central panels) and small-world index (right panels). of all layers (Figure 5A). The difference between connectivity in different layers is quantified by Jensen-Shannon divergence (see Materials and Methods), a powerful information-theoretical measure of (dis)similarity. A quality function controls the reduction process and its global maxima identify optimal structural reduction strategies.
We perform structural reducibility for each individual separately and calculate the corresponding quality functions, for control and patients groups ( Figure 5B). In both cases, we found that the maximum value of the quality function is attained when no reduction is performed at all, providing evidence that the topological information carried by each functional network, corresponding to a different frequency component, should not be disregarded from structural analyses. It is worth noting that the behavior of the quality function alone does not allow to distinguish between the two groups of individuals.
To gain insights about (dis)similarities between different layers of the multiplex functional network in the two groups, we use the quantum Jensen-Shannon distance (see Material and Methods) calculated during the structural reducibility analysis. The distance matrix, whose entries provide the Jensen-Shannon distance between any pairs of layers, is first built for each individual separately, and group average µ and standard deviation σ are calculated. The signal-to-noise ratio (SNR) defined by their ratio is successively calculated for each pair of layers and for each group, separately (Figure 5C), as well as the relative difference between the two group-averaged values. We observed differences of up to 30% in absolute value between the two groups, for specific pairs of layers. Dissimilarities between layers within the typical-band were higher in healthy individuals than in schizophrenic patients. On one hand, functional connectivity in healthy subjects is rather volatile and, in general, exhibits topological differences across individuals (Sasai et al., 2014) that we did not observe in patients, suggesting the possibility that schizophrenia might alter brain's integrated operation to reduce such a functional diversity. On the other hand, an abnormal amount of dissimilarity between functional networks corresponding to other frequency bands (such as those within relatively higher ranges, e.g., 0.09-0.19 Hz) was observed in patients but not in healthy individuals. it allows to identify frequency bands providing redundant topological information and to verify the validity of the multiplex model with respect to conventional single-layer models. Global maxima in the quality function identify optimal structural reductions. (B) The median quality function is shown for healthy control (solid) and schizophrenic patients (dashed), with shaded areas indicating the standard deviation around each value. (C) Signal-to-noise ratio (SNR; see text for further details) for Jensen-Shannon distance calculated for each pair of layers, color-coded for both groups, and corresponding relative difference between the two groups.
These results suggest that the dependence on frequency of patients' functional connectivity is different from that of healthy individuals and we might use such dissimilarity patterns as a fingerprint of brain's functional organization for each group.

Identifying Schizophrenic Patients by the Centrality Profile of Their Brain Functional Connectivity
The importance of a region with respect to the overall brain functional connectivity can be quantified by centrality descriptors Sporns et al., 2007;Bassett et al., 2008;Lynall et al., 2010;Power et al., 2013;Rubinov and Bullmore, 2013;van den Heuvel and Sporns, 2013;Thompson and Fransson, 2015). Here, we propose to use PageRank (Brin and Page, 1998;Ermann et al., 2015) as a measure of centrality, which is based on the rationale that nodes linked by influential nodes are more central than those linked by un-influential nodes. It has been used in several applications, from ranking relevant Web pages in the World Wide Web (Brin and Page, 1998) to identifying important nodes in the human functional connectome (Zuo et al., 2012) and, more recently, in a variety of multiplex networks (De Domenico et al., 2015c). Once centrality scores are calculated for each node, the set of all their values constitutes the centrality profile of the underlying functional network. The Spearman's correlation coefficient between the centrality profiles corresponding to the multiplex, full-band and typical-band functional networks is calculated (Figures 6A,B). While very strong correlations are observed for centrality profiles calculated from single-layer networks, no significant correlation with multiplex centrality profiles were found.
These results suggest the appealing possibility to use the multiplex centrality profiles to gain new insights about brain functional connectivity. To this aim, we interpret the centrality profiles as characteristic features of each individual (control or patient) and we use the well-known and robust random forest method (Breiman, 2001) to train a classifier distinguishing between healthy and schizophrenic individuals (see Material and Methods for further details). At the very beginning, we trained the classifier by using all 264 centrality scores available for each individual and found a classification accuracy of about 60-65%, regardless for the type of centrality profile used (i.e., multiplex, full-band and typical-band). One of the main advantage of random forest classification is that it also ranks the features based on their classification power, i.e., on the degree of discrimination they have. We capitalize on this precious information to perform FIGURE 6 | Comparing centrality profiles of multiplex and conventional functional networks. Spearman's correlation coefficient between the centrality profiles obtained from multiplex, full-band and typical-band functional networks for (A) healthy and (B) patient groups. a second round of classification, this time using only top-ranked features instead of the full set. We varied between 10 and 260 the number of top features used to discriminate between healthy and schizophrenic individuals. The comparison between the results obtained from different centrality profiles are shown in Figure 7A (see Supplementary Figures 1-3 for further details) for top features varying between 10 and 50, being the accuracy of the classification rapidly decreasing for increasing size of centrality profiles. Remarkably, multiplex centrality profiles allow a more accurate discrimination of the two populations, confirming the hypothesis that multiplex functional networks are able to encode richer information than their single-layer aggregations, providing a suitable framework to better discriminate healthy brains from schizophrenic ones. This result is robust against the selection of the number of features used to discriminate, with the multiplex approach significantly outperforming the other ones.
To gain further insights, we focus on the regions corresponding to the top 30 ROIs of the multiplex centrality profile, where we attained the maximum discrimination between control and patient groups. The spatial distribution of the corresponding brain's regions are shown in Figure 7B. Anatomical information about these ROIs is summarized in Supplementary Table 1.

Characterizing Regions Distinctive of Schizophrenic Brain Functional Activity
Capitalizing on results from group discrimination using centrality profiles, we investigate more in detail the role of hub regions. In particular, our interest is twofold. On one hand, we wonder if the most central regions obtained from multiplex and conventional functional networks are the same (it is worth remarking that previous correlation analysis of centrality profiles does not provide this information, because differences might be due to low-ranked regions, for instance). On the other hand we want to clarify which hub regions are found only in healthy individuals, which ones are found only in schizophrenic individuals and which ones are found in both groups.
Hubs were identified as ROIs ranked in top 5% in terms of group-averaged region centrality. Figure 8 shows the spatial distributions in the brain of such hubs, for each group, while the corresponding anatomical information is reported in Supplementary Table 2. In all cases we found hub regions peculiar for each group and hubs regions that are common to both groups. While significant differences are not observed between networks built from conventional approaches, hubs from multiplex analysis constitute a distinct set.
In both conventional networks, healthy-specific hubs are located in cuneus, precuneus, transvers temporal, superior temporal, and superior parietal cortices, whereas schizophrenicspecific hubs are located in superior frontal, middle frontal, and precentral cortices as well as thalamus. Hubs shared by both groups are identified along the midline of the brain, in particular the medial superior frontal, precuneus and cingulate cortices. In the multiplex network, healthy-specific hubs controls are located in anterior cingulate, superior frontal, insula and superior temporal cortices, whereas those pertaining to schizophrenic  patients are distributed over frontal, parietal and occipital cortices. Hub regions shared by both groups groups are localized in frontal, occipital cortices and cerebellum. Notably, no hub region has been identified in the precuneus cortex, a region well known to function as a hub in healthy individuals (van den Heuvel and Sporns, 2013).
Finally, we compared the top 30 hubs found in each group against the top 30 distinctive ROIs found from the classification procedure previously described. In the case of random forest classification, we ranked ROIs by their average rank to identify the most discriminating ones; let us indicate by R class the set of the top 30 ROIs according to this ranking. Similarly, we first rank the ROIs by their PageRank versatility and then calculate their group-averaged ranks; let us indicate by R H pr and R S pr the sets of the top 30 ROIs found in healthy and schizophrenic groups, respectively, according to this ranking. In the case of healthy subjects, 3 hubs in R H pr are found to be also discriminating features; the Jaccard index of the two sets is 0.053. Remarkably, in the case of schizophrenic subjects, 13 hubs in R S pr are found to be also discriminating features; the Jaccard index of the two sets is 0.277. The ROIs corresponding to the two cases are reported in Supplementary Table 4. This finding confirms that highly central regions in schizophrenic brains are different from highly central regions in healthy brains, and that this result can be used in practical applications to better identify patients affected by schizophrenia.

DISCUSSION
Resting state functional connectivity has been widely investigated with fMRI in the past two decades. Since the first study conducted by Biswal et al. (1995), functional connectivity has been defined as an inter-regional temporal correlation of fMRI signals that are preprocessed with band-pass filters, removing frequency components below 0.01 and above 0.1 Hz. In fact, the power spectrum of spontaneous fluctuations of fMRI signals roughly follows a 1/f power-law scaling (He, 2011), where powers in the higher frequency range are relatively weaker than lower ones, suggesting the hypothesis that only the lower frequency range substantially contributes to brain's function. However, recent studies have reported that conventionally excluded frequency bands might provide additional insights on brain activity Liao et al., 2013;Sasai et al., 2014;Chen and Glover, 2015;Thompson and Fransson, 2015). As a consequence, brain functional networks exhibit a peculiar architecture, consisting of a few regions acting as hubs, strongly dependent on the frequency components of brain activity that contribute to interregional interactions. However, a rigorous method to identify such hubs in networks built from different frequency bands simultaneously is a challenging problem remaining largely unexplored.
Our results, based on multiplex modeling and analysis of the brain activity, provide convincing evidence that characterization of brain functional networks can not prescind from considering the whole information observed from different frequency bands, simultaneously. This crucial finding allows to exploit new theoretical and computational tools for the analysis of brain activity and opens a new direction toward a deeper understanding of brain function and its operated integration. As a first hint of the power of the new methodology, we have shown that multiplex characterization of brain regions, in terms of network centrality, allows to find new areas of the brain that have never been classified as relevant in brain's functional integration (or the opposite). This is the case of ROIs in the precuneus cortex, a well-known region of highly central functional hubs (van den Heuvel and Sporns, 2013), that are not found by our multiplex network analysis, reflecting the importance of considering the whole information simultaneously, rather than aggregating or neglecting part of it (De Domenico et al., 2015c).
We wondered if this result could be exploited for practical applications, where the choice of specific frequency bands might play a crucial role. We focused our attention on characterizing brain disorders in schizophrenic patients, a research topic of great interest that has been largely explored (Bassett et al., 2008;Lynall et al., 2010;van den Heuvel et al., 2013), although individual diagnosis based on brain imaging remains still undeveloped (Rubinov and Bullmore, 2013). With the aid of the MRI technique, it has been recently shown that regions affected by schizophrenia are distributed across the brain (Glahn et al., 2008;Ellison-Wright and Bullmore, 2009), impelling researchers to move from the conventional perspective where the causes of disorders are localized in specific areas, to a wider perspective with emphasis on abnormality in brain structural and functional connectivity (van den Heuvel and Fornito, 2014). Studies on structural connectivity provided evidence that schizophrenic brains exhibit abnormal network architecture, characterized by reduced hierarchical organization, the loss of frontal hubs with emergence of non-frontal hubs (Bassett et al., 2008;Lynall et al., 2010) and degraded rich-club organization (van den Heuvel et al., 2013). Methods not based on networks were able to provide satisfactory performance in discriminating schizophrenic patients from the analysis of their brain activity (Yang et al., 2010;Chyzhyk et al., 2015), although they are often based on very complicated machine learning algorithms and make use of heterogenous data sources, thus not improving our understanding of brain function. Here, we have found that multiplex centrality profile of brain regions allow to discriminate between control and schizophrenic groups of individuals more accurately than centrality profiles calculated from networks obtained by using conventional approaches, such as aggregating and/or disregarding the measured activity. Nevertheless, the discrimination accuracy is comparable to other methods, with the additional advantage of providing a framework facilitating the interpretation of results, without relying on external data sources or phenotypic information. In fact, we were able to identify many regions distinctive of schizophrenic brains, some of them localized where abnormality has been previously suggested (Honea et al., 2005;Rubinov and Bullmore, 2013). The analysis of dissimilarities between networks corresponding to different layers of the multiplex functional network, confirmed significant differences between healthy and schizophrenic individuals in specific frequency ranges, including the higher ones. This finding demonstrates that brain activity in higher frequencies provides unique information about functional interaction in the brain, even if their amplitudes are under-represented in the power spectrum.
Nevertheless, the present study presents some limitations. First, the sampling rate of fMRI signals (0.5 Hz) is lower than the values recently used to investigate frequency-specificity of functional connectivity (Wu et al., 2013;Gohel and Biswal, 2015). We think that a future study would benefit from taking into account a wider frequency range, as the one provided by a higher sampling rate. Second, neural mechanisms generating frequency components of fMRI and their interactions still remain unclear. Electrophysiological signals, as well as fMRI signals, include many different frequency components showing distinct network topology (Siegel et al., 2012). The fact that electrophysiological studies have repeatedly shown cross-frequency coupling as a mechanism of interactions between different frequency layers (Canolty et al., 2006;Tort et al., 2009;Axmacher et al., 2010;Belluscio et al., 2012), we conjecture that inter-layer interactions between functional networks built in the present study may reflect the corresponding mechanism in the case of fMRI signals. The differences between healthy and schizophrenic functional brains found in our study are related to some layers, suggesting that pathological abnormality in schizophrenia may occur on neural mechanisms with specific frequency-dependent fingerprints. Further studies, for instance with simultaneous recording of electroencephalography and fMRI, could allow us to examine the electrophysiological background of fMRI frequency components and possible mechanisms of interactions between different components.
The proposed methodology suggests a guideline for future studies designed to consider brain's inter-regional interactions at different frequencies, encouraging the application of other multiplex network measures to functional networks obtained, for instance, from variable brain states.

AUTHOR CONTRIBUTIONS
MD and SS contributed equally to this work. MD and SS analyzed the data and performed the analysis. MD, SS, and AA designed the study and wrote the paper.