Frequency-based brain networks: From a multiplex framework to a full multilayer description

During the last few years, network science has undergone a conceptual revolution with the extension of well-established techniques of network analysis to multilayer brain networks (Boccaletti et al., 2014; De Domenico et al., 2013; Kivelä et al., 2014), which provide a convenient way to simultaneously encode different types of interactions, subsystems, and other complications in networks. Consequently, it has been necessary to revisit our intuitive understanding of both structural and dynamical properties of networks—including structural phase transitions (Radicchi & Arenas, 2013), diffusion and other spreading processes (De Domenico, Granell, Porter, & Arenas, 2016; Gómez et al., 2013; Salehi et al., 2015), percolation and robustness (Buldyrev, Parshani, Paul, Stanley, & Havlin, 2010; Gao, Buldyrev, Stanley, Xu, & Havlin, 2013), synchronization (Aguirre, Sevilla-Escoboza, Gutiérrez, Papo, & Buldú, 2014), and others—to the new possibilities in multilayer descriptions, leading in many cases to counterintuitive results.

The study of brain networks is currently undergoing a process of adaptation of classical single-layer (“monolayer”) concepts and analyses to a more general multilayer description. (For reviews, see Betzel & Bassett, 2016; De Domenico, 2017; Vaiana & Muldoon, 2018; also see Figure 1 of Betzel & Bassett, 2016 and Figures 1–2 of De Domenico, 2017 for relevant schematics.) Some studies have considered integration of data from structural and functional brain imaging into a multilayer network to account for both anatomical and dynamical information. In an early study using monolayer networks, Honey, Koötter, Breakspear, and Sporns (2007) showed that a dynamical model simulated on the anatomical network of a macaque neocortex can successfully identify the positions of the anatomical hubs when signals are averaged appropriately. More recently, Stam et al. (2016) analyzed how anatomical networks support activity, leading to specific functional networks (either undirected ones or directed ones), demonstrating that a dynamical model close to a critical transition is able to unveil interactions between structural and functional networks.

It is also possible to combine anatomical and functional interactions of just a few nodes instead of an entire brain network. For example, Battiston, Nicosia, Chavez, and Latora (2017) examined network motifs (i.e., overrepresented network substructures that consist of a few nodes; Milo et al., 2002) that combine anatomical connections (in one layer) and functional relations between cortical regions (in a second layer), linking data obtained, respectively, from diffusion-weighted magnetic resonance imaging (DW–MRI) and functional magnetic resonance imaging (fMRI).

Efforts to combine anatomical and functional networks into a single multilayer network face the challenge of how to normalize the weights of edges that arise from different origins. To tackle this issue, Simas, Chavez, Rodriguez, and Díaz-Guilera (2015) proposed translating functional and anatomical networks into a common embedding space and then comparing them in that space. They constructed functional networks (each with N nodes) from the fMRIs of q healthy individuals. They then used the functional networks of the q = 20 individuals to construct a single functional multiplex brain network, a special type of multilayer network in which corresponding entities (brain regions) in different layers (individuals) can be connected to each other via interlayer edges, but other types of interlayer edges cannot occur (Kivelä et al., 2014). They followed a similar procedure to construct an anatomical multiplex network using data obtained from DW–MRI. They then independently projected each of the multiplex networks into a common embedding space using a series of algebraic operations that allow one to calculate an “algebraic aggregation” of all layers into a single layer (see Simas et al., 2015, for details). Using such a projection, it is possible to quantify the differences between anatomical and functional networks. Simas et al. also calculated an “averaged aggregation” for each of the functional and anatomical multiplex networks by averaging the weights of the corresponding edges over all layers. They compared the two types of aggregation, and they were thereby able to identify certain brain regions (e.g., ones related to visual, auditory, and self-awareness processes) with significant differences between the functional and anatomical networks for both types of aggregation. However, only the algebraic aggregation was able to detect differences between the functional and anatomical networks in other regions (specifically, the thalamus, the amygdala, the postcentral gyrus, and the posterior cingulate), suggesting that the averaged aggregation disregards significant information (Simas et al., 2015).

One possible alternative for reducing the complexity of analysis of brain networks is to concentrate only on “functional” (dynamical, in fact) interactions between brain regions and to define multilayer functional networks as a concatenation of a series of layers, each of which captures the interplay between brain regions during some time window. This approach, in which a layer in a multilayer network represents connection similarities over some time window, was taken in papers such as Bassett et al. (2011, 2013) and Braun et al. (2014) to analyze the temporal evolution of network modules and examine dynamical reconfiguration and “flexibility” of functional networks.

Another alternative is to construct functional multilayer networks whose layers correspond to the well-known frequency bands at which a brain operates (Buzsáki, 2006). As demonstrated by Brookes et al. (2016), it is possible to construct frequency-based multilayer networks from magnetoencephalographic (MEG) recordings by (i) band-pass filtering the raw MEG signals, (ii) obtaining the envelope of the amplitude at each frequency band, and (iii) measuring the pairwise correlations between each envelope (for whichever frequency it accounts). Using such an approach, Brookes et al. (2016) constructed frequency-based multilayer networks, in which each layer includes the interactions in a given frequency band, and showed that the corresponding supra-adjacency matrices (which encode a linear-algebraic representation of connections in a multilayer network) convey statistically significant differences when comparing a control group with a group of schizophrenia sufferers.

Recently, De Domenico, Sasai, and Arenas (2016) took the important step of analyzing the spectral properties of matrices associated with frequency-based multiplex networks. The spectrum of an $Ñ×Ñ$ matrix consists of the set ${λ1,λ2,λ3,…,λÑ}$ of its eigenvalues, and it encodes valuable information about the structural properties of the corresponding network. In turn, these eigenvalues (as we will explain later) are related to various dynamical properties, such as network synchronizability, robustness, and diffusion (Newman, 2010). De Domenico, Sasai, and Arenas (2016) compared a group of schizophrenic patients with a control group by using fMRI data, and they found that the second-smallest eigenvalue (i.e., the algebraic connectivity or “Fiedler value”; Van Mieghem, 2011) λ₂ of the combinatorial supra-Laplacian matrix¹ associated with a multiplex network is a better discriminator between the two groups than what one can obtain by studying either unfiltered or single-band functional networks (i.e., by using monolayer networks). De Domenico et al. also calculated centrality measures (i.e., measures of the importance of network components; Newman, 2010) on the frequency-based multiplex networks to demonstrate the existence of hubs that had not been classified previously as important brain regions for functional integration. Hubs of the control group were located in anterior cingulate, superior frontal, insula, and superior temporal cortices; however, hubs for schizophrenic patients were distributed over frontal, parietal, and occipital cortices. These results revealed that frequency-based multiplex networks include relevant information about the functional organization of brain networks that is not captured by using a classical monolayer approach.

Multiple recent papers have demonstrated the benefits of using a multiplex description for analyzing the functional networks of patients suffering from Alzheimer’s disease (AD). For example, Yu et al. (2017) constructed frequency-based multiplex networks from MEG data and demonstrated that centrality measures calculated with layers analyzed independently are unable to detect significant differences between a control group and the AD group. However, when centralities are evaluated on a frequency-based multiplex network, one can find statistically significant differences in the hippocampus, posterior default-mode network, and occipital areas. Guillon et al. (2017) also used frequency-based multiplex networks to differentiate between controls and people with AD. In this case, the authors proposed the use of a multi-participation coefficient (MPC) to enhance classification of which individuals are suffering from AD. Their MPC consists of an adaptation of the classical (i.e., monolayer) participation coefficient (Guimerà & Amaral, 2005) to networks composed of several layers (see also Yu et al., 2017). A potential benefit of using this MPC is that it does not depend on interlayer edge weights, which thus do not need to be calculated. As shown in Guillon et al. (2017), using an MPC yields better classification accuracy and sensitivity than using only monolayer network diagnostics.

In this paper, we investigate how to translate the dynamics of different brain regions into a frequency-based (i.e., functional) multilayer network, in which individual layers account for coordination within a given frequency band. We focus specifically on the consequences of analyzing a multiplex network versus a more general multilayer one. The former allows interlayer connections only between the same brain region (i.e., node-layer) in different network layers, so coupling between oscillations in different frequency bands occurs only between the same brain region, whereas the latter allows one to model coordination between any brain region at any frequency band. We use resting-state MEG recordings because of their high temporal resolution (on the order of milliseconds), which makes it possible to analyze a broad spectrum of frequency bands (de Pasquale et al., 2010; van Diessen et al., 2015). In our case, a set of MEG signals consist of N time series, each of which comes from a sensor that captures the activity above a different cortical region. We then filter signals at four frequency bands (theta, alpha, beta, and gamma) and construct, for each individual, a four-layer functional multilayer network from the dynamical coordination within and between frequency bands. The existence of interlayer edges in frequency-based functional networks relies on the phenomenon of cross-frequency coupling (CFC), which is responsible for integrating brain activity at different spatial and temporal scales (Canolty et al., 2006). The quantification of CFC is a hard task, because the interplay between frequency bands is very intricate (Aru et al., 2015). To address this issue, one can set the weights of the interlayer edges to a value p, which one can estimate by using an optimization procedure. For example, it is possible to construct multiplex networks of two groups of individuals with different profiles and determine a value of p that best distinguishes between the two groups (see, e.g., De Domenico, Sasai, & Arenas, 2016). However, we adopt a different strategy: We obtain the weights of the interlayer edges directly from time series by calculating mutual information (MI; MacKay, 2003). We thereby capture heterogeneity (in the weights) of both intralayer edges and interlayer edges, and we investigate the influence of such heterogeneity on the spectral properties of frequency-based functional networks.

Using both synthetic network models and data from laboratory experiments, we investigate the effects that heterogeneity of interlayer edge weights have on the spectral properties of both multiplex and more general multilayer networks. Specifically, we focus on the algebraic connectivity λ₂, which is closely related to both structural and dynamical properties of networks (Masuda, Porter, & Lambiotte, 2017; Newman, 2010; Van Mieghem, 2011). On one hand, algebraic connectivity is an indicator of modular structure in networks (Fortunato & Hric, 2016). In the framework of multilayer networks, one can interpret the value of λ₂ and how it changes as a function of interlayer coupling strength as a way to quantify structural integration and segregation of different network layers (Radicchi & Arenas, 2013); also see the discussions in DeFord and Pauls (2017). On the other hand, 1/λ₂ is proportional to the time required to reach equilibrium in a linear diffusion process (Gómez et al., 2013). Additionally, the time t_sync to reach synchronization of an ensemble of phase oscillators that are linearly and diffusively coupled is also proportional to 1/λ₂, and it is known that t_sync and 1/λ₂ are positively correlated in some situations with nonlinear coupling (Almendral & Díaz-Guilera, 2007). A recent investigation of the consequences of modifying interlayer edge weights in multiplex networks has illustrated that enhancing interlayer coupling tends to shorten the time to achieve interlayer synchronization in a Kuramoto model of coupled oscillators (Allen-Perkins, Albuquerque de Assis, Pastor, & Andrade, 2017).

In the framework of functional brain networks, algebraic connectivity has been used as an indicator of AD, such as in de Haan et al. (2012), who obtained statistically significant differences for λ₂ of functional networks obtained from MEG in a comparison of a group of patients suffering from AD with a control group. Phillips, McGlaughlin, Ruth, Jager, and Soldan (2015) calculated the value of λ₂ in a group of individuals with mild cognitive impairment and AD, although they did not report significant differences between them. Computing λ₂ is also necessary for calculating the most standard type of synchronizability parameter, which is the ratio λ_N/λ₂, where λ_N is the largest eigenvalue of the combinatorial Laplacian matrix (Boccaletti, Latora, Moreno, Chavez, & Hwang, 2006). In a series of studies, the synchronizability of brain networks was calculated for different frequency bands (Bassett, Meyer-Lindenberg, Achard, Duke, & Bullmore, 2006), during epileptic seizures (Schindler, Bialonski, Horstmann, Elger, & Lehnertz, 2008), and for schizophrenic individuals (Siebenhühner, Weiss, Coppola, Weinberger, & Bassett, 2013). We aim to improve the interpretation of algebraic connectivity for functional brain networks, and we thus investigate (i) how the fact that a considerable fraction of all possible interlayer edges are not present in multiplex networks leads to a deviation from the theoretical values expected for λ₂ and (ii) how these deviations are related to the mean weight of the interlayer edges. We thereby scrutinize the consequences of using a multiplex formalism, in which only CFC inside the same brain region is allowed, instead of employing a fully multilayer approach (i.e., without any restrictions on the type of coupling that one considers).

In Figure 1, we illustrate the process of constructing frequency-based multiplex and multilayer brain networks. Our starting point is a data set of MEG recordings of a group of q = 89 individuals during resting state (see Materials and Methods for details), but one can use other experimental paradigms—including different brain-imaging techniques, such as electroencephalography (EEG) or fMRI—to construct multilayer networks with the same procedure. Specifically, we record MEG activity at N cortical regions, with 235 ≤ N ≤ 246 (see the Supplementary Information for details; Buldú & Porter, 2017), and we then clean the data to remove artifacts and obtain corresponding unfiltered signals. We thereby analyze the signal recorded by each sensor instead of carrying out a source reconstruction. We then band-pass filter each signal to obtain four different filtered time series for each brain region. We use the four classical frequency bands: theta [3–8] Hz, alpha [8–12] Hz, beta [12–30] Hz, and gamma [30–100] Hz. The number l of layers of the multilayer network is the number of different frequencies that we examine (and 2 ≤ l ≤ 4 in our work), and the N nodes in each layer are associated with the dynamics of the N sensors filtered at a given frequency band. We number the nodes so that nodes n, n + N, n + 2N, …, n + lN (with n ∈ {1, … N}) correspond to the signals of the same brain region n at the l different frequency bands (i.e., layers).

We then quantify coordination between each pair of nodes of a multilayer network, regardless of which layers they are in, using mutual information (see Materials and Methods). Calculating mutual information (MI) between time series of the same frequency band yields intralayer connections between brain regions (see edge type “1” in the bottom-left plot of Figure 1A for an example), so each layer corresponds to a specific frequency band. Edges between the signals of the same sensor at different frequency bands result in interlayer connections between layers (see edge type “2”). Such “diagonal” interlayer edges are the only type of interlayer edges that are allowed in multiplex networks (Kivelä et al., 2014). Finally, cross-frequency coupling between different brain regions yields the other (“nondiagonal”) interlayer edges in a full multilayer network (see edge type “3”). As we show in Figure 1B, we thereby obtain a supra-adjacency matrix, where blocks along the diagonal account for intralayer connections (layers alpha and beta in the two-layer example) and blocks off of the diagonal, marked as P^αβ and P^βα, encode the interlayer edges. Because MI_ij = MI_ji (see Equation 4 in Materials and Methods) for time series of nodes i and j, the supra-adjacency matrix is symmetric, so P^αβ = [P^βα]^T.

Importantly, although we have chosen to use MI, there are a diversity of similarity measures for capturing amplitude–amplitude and phase–amplitude correlations between different frequency bands (see Aru et al., 2015 for a review of cross-frequency coupling measures), and each measure has its own advantages and drawbacks. Nevertheless, as we will see, the same methodological implications exist no matter which specific measure one uses to evaluate coordination between brain sites.

In our discussion, we focus on the analysis of a two-layer network with alpha and beta layers. We will discuss the consequences of considering alternative frequency bands and numbers of layers in the last section of the paper (see Conclusions and Discussion).

Our starting point is to compare the results from four different constructions of functional networks:

1. Unfiltered functional networks. We obtain these networks from the original (unfiltered) signals of each brain region—that is, without decomposing the signals into different frequency bands—so these are monolayer networks.

2. Aggregated networks. We obtain these networks from componentwise addition of the weights of the alpha and beta layers to form monolayer networks. (Note that this is a uniform aggregation.)

3. Multiplex networks. Each layer corresponds to a specific frequency band (as explained above), and interlayer edges are allowed only between node-layers that correspond to the same brain region.

4. Full multilayer networks. These include the same layers as their multiplex counterparts, but they allow all possible interlayer edges.

In Figure 2, we show the probability density functions (PDFs) for a group of q = 89 individuals of the values of λ₂, the standard deviation of the interlayer edge weights (to quantify their heterogeneity), and the percentage of missing interlayer edges (also see Table 1).

In column one of Figure 2, we observe that the (monolayer) unfiltered functional network has a similar mean and standard deviation of λ₂ as the multilayer network (also see Table 1). However, the aggregated networks tend to have larger values of λ₂, which makes sense, as we construct an aggregated network by adding the weights of the two layers (alpha and beta), and the total “strength” of the network (i.e., the sum of all of its edge weights) is close to double the strength of each layer. Importantly, the mean λ₂ for the multiplex networks is two orders of magnitude smaller than the mean λ₂ of the multilayer networks. We expect this discrepancy, because we construct a multiplex network by deleting all interlayer edges of a full multilayer network, except for edges (so-called “diagonal” edges) that link the same “physical” nodes (i.e., the same brain regions). Thus, the total strength of the interlayer matrix P^αβ is considerably smaller in multiplex networks than in corresponding full multilayer networks. Because λ₂ is an indicator of the amount of interconnections between communities in a network (Newman, 2010), one expects such drastic edge removals to yield a smaller value of λ₂, as one can construe layers as communities with a small number of edges between them (only N of the N² possible interlayer edges of a full multilayer network).

In columns two and three of Figure 2, we quantify the heterogeneity and the number of missing interlayer edges of the four different functional networks. In column two, we plot the PDFs of the standard deviation of all edges (for the unfiltered and aggregated networks) and interlayer edges (for the multiplex and multilayer networks). In all cases, we observe that the functional networks that we construct from experimental data have non-negligible heterogeneity (also see Table 1). Again, aggregated networks have values that are roughly twice as large as those of the other kinds of networks because of the (rough) doubling of the mean strength. In column three, we show that there is a large percentage of missing edges. This arises from the fact that we construe edges whose weights are smaller than those obtained with an appropriate surrogate time series as not being statistically significant, so we set their values to 0 (see Materials and Methods for details). The (monolayer) aggregated networks have the smallest mean percentage, followed by the (monolayer) unfiltered networks, and then the two types of multilayer networks. For the multiplex and full multilayer networks, the percentage of missing edges, which is larger than 40% in both cases, refers to the number of all possible interlayer edges. Note that missing edges are unavoidable in functional brain networks, because not all brain regions communicate with each other through direct connections². Moreover, the different amounts of coordination between pairs of brain regions lead to functional networks with heterogeneous weights.

Given our prior observations, a crucial question arises: What are the consequences of heterogenous and missing interlayer edges, both of which are intrinsic features of brain-imaging data, on multiplex and full multilayer functional networks? More specifically, how do they affect the value of λ₂, and what are the ensuing structural and dynamical implications of these differences?

With the aim of answering the above questions, we perform a series of numerical computations in which we compare the theoretical values of λ₂ in multiplex and full multilayer networks from homogeneous interlayer-edge distributions with ones from networks with heterogeneous and missing interlayer edges.

In Figure 3, we show the consequences of heterogeneity in the distribution of the interlayer-edge weights of the multiplex network M^αβ. In this example, the multiplex network has an interlayer connectivity matrix $Pαβ=p→I$⁠, where $I$ is an N_α × N_α (equivalently, N_β × N_β, as N_α = N_β in this example) identity matrix and $p→=ph→$ is a row vector that controls the weights of interlayer edges, where p modulates the vector’s magnitude and the vector $h→$ encodes the heterogeneity of the interlayer edges. We set the elements of $h→$ to follow a uniform distribution over the interval [h_min, h_max]. These elements have a mean value of $h-=1$ and a standard deviation of σ. We set $h-=1$ and construct networks with interlayer-edge-weight heterogeneities that range from σ = 0 (blue circles) to σ ≈ 0.581 (green circles). We then examine the interplay between the weights of the interlayer edges and the heterogeneity by increasing the value of p. Note that σ = 0 corresponds to what we call a homogeneous multiplex network, which has uniformly-weighted interlayer edges (i.e., $p→=p$ for all interlayer edges). We obtain the results in Figure 3 by calculating a mean over 100 realizations of two-layer networks with the G(N, p_con) Erdős–Rényi (ER) model (with p_con = 0.25) in each layer and N = N_α = N_β = 250 nodes (Newman, 2010).

As explained in Radicchi and Arenas (2013), modifying the weight parameter p of the interlayer edges has important consequences for the value of λ₂ for homogeneous multiplex networks. The existence of two regimes of qualitatively distinct dynamics, separated at a transition point p*, was discussed in Radicchi and Arenas (2013) (and in various subsequent papers): When p ≪ p*, the algebraic connectivity λ₂ follows the linear relation λ₂ = 2p; for p ≫ p*, however, the value of λ₂ approximates that of the (uniformly) aggregated network [i.e., $λ2,agg=12λ2(Lα+Lβ)$]. In Figure 3, the dashed and solid lines indicate the theoretical predictions of λ₂ in the homogeneous case for small (dashed) and large (solid) values of p. We obtain a mean value of p* ≈ 2.870 by locating the intersection between the eigenvalues λ₂ and λ₃ of the supra-Laplacian matrix $L$^αβ (Sahneh, Scoglio, & Van Mieghem, 2015).

In Figure 3A, we illustrate that the heterogeneity of the (weights of the) interlayer edges leads to non-negligible differences in λ₂. Furthermore, values of p near p* have the maximum discrepancy between heterogeneous cases (colored circles) and the homogeneous case (black circles). In studies of multiplex networks, it is simplest to take interlayer edge weights to be homogeneous (i.e., given by a constant value p), especially when it is not clear how to estimate their values (Kivelä et al., 2014). When possible, however, it is desirable to be more sophisticated, and it is possible to quantify the coupling strength between different frequency bands in various ways when studying frequency-based functional brain networks (Aru et al., 2015). The choice of incorporating versus ignoring heterogeneity of inter-frequency edge weights then leads to unavoidable differences in the estimation of λ₂, especially near the transition point p*. It is also important to think about the different interlayer edge weights that one obtains from different measurements of coupling strength (see Materials and Methods).

Discrepancies in the values of λ₂ are typically even larger when some interlayer edges are missing. In Figure 3B, we show the effects of removing some percentage of the interlayer edges in a “holey” multiplex network (i.e., a multiplex network with an interlayer connectivity matrix P^αβ whose elements are either 1 or 0 along the diagonal and are 0 in all other entries). We observe that increasing the number of missing interlayer edges causes the multiplex networks to have drastically reduced values of λ₂. When p ≪ p*, the value of λ₂ always grows with a slope that is smaller than 2p; for p ≫ p*, however, the value of λ₂ never reaches the value of λ₂ for the aggregated network. Finally, values of p close to p* again have the largest discrepancies with respect to the homogeneous case.

We now identify the qualitative differences (and some of their consequences) between multiplex networks and full multilayer networks. The latter have an interlayer connectivity matrix P^αβ = pC, where C has elements c_ij that account for the weight between each pair {i, j} of nodes i ∈ α and j ∈ β, and the parameter p allows one to modulate the mean weight of the interlayer edges. In Figure 4A, we connect the same layers as in our previous example, but now we use a homogeneous interlayer connectivity matrix P^αβ with weights p_ij = p (i.e., c_ij = 1 for all i and j). Increasing p from 0 leads to the transition point p* ≈ 0.013, which we obtain from the intersection of the eigenvalues λ₂ and λ₃ of the combinatorial supra-Laplacian $L$^αβ. We observe for p < p* that the value of λ₂ follows the linear function λ₂ = 2p〈c〉 (see the dashed line in Figure 4A), where 〈c〉 = 250 is the mean weighted degree of the matrix C, as demonstrated in Radicchi and Arenas (2013). Interestingly, after the transition point p = p*, one can describe the value of λ₂ by the function λ₂ = λ_2,min(α,β) + λ₂(pC), where λ_2,min(α,β) is the value of the smaller λ₂ of the two isolated layers. As we show in the inset of Figure 4A, when we introduce a certain amount of heterogeneity into the interlayer connectivity matrix, we observe slight differences from homogeneous case. Specifically, for heterogeneous cases, these differences are larger for values of p above the mean 〈p*〉 of the transition point.

As we saw in our analysis of multiplex networks, deleting interlayer edges from a full multilayer network increases the discrepancy in λ₂ values between networks with homogeneous interlayer edges and networks with heterogeneous ones. As we show in Figure 4B, the value of λ₂ decreases with the number of missing edges, which one should expect, because having a smaller number of interlayer edges implies that one needs larger edge weights to maintain the same amount of interlayer coupling. Interestingly, the deviation from the theoretical predictions is significant even for p < 〈p*〉. Additionally, for p > 〈p*〉, we observe (as expected) a change in the slope of λ₂, but the theoretical predictions given by λ₂ = λ_2,min(α,β) + λ₂(pC) (solid line) begin to fail, leading to a discrepancy that increases with the percentage of missing edges.

Because full multilayer networks have up to N² interlayer edges, whereas multiplex networks can have only N of them, the former tend to have interlayer connectivity matrices with larger strengths $SP=∑ijpij$⁠. In Figure 5, we show the algebraic connectivity λ₂ of the combinatorial supra-Laplacian matrices $L$^αβ for a series of networks with identical layers and interlayer strength S_P, but with a different number of interlayer edges, ranging from a multiplex network (with N interlayer edges) to a full multilayer network with no nonzero entries (i.e., with N² interlayer edges). The dashed lines correspond to the analytical solutions for the full multilayer network (black dashed line; λ₂ = 2p〈c〉, with 〈c〉 = 250) and the multiplex network (red dashed line; λ₂ = 2p) when p < p*. The solid curves are the corresponding theoretical solutions when p > 〈p*〉 for the full multilayer network (black line; λ₂ = λ_2,min(α,β) + λ₂(pC)) and the multiplex network (red curve; $λ2,agg=12λ2(Lα+Lβ))$⁠. We observe the effect that adding interlayer edges to multiplex networks has on the value of λ₂ and the associated transition from a multiplex network to a full multilayer architecture. Interestingly, the different numbers of interlayer edges in the two types of networks leads to a difference in the position of 〈p*〉 (which one can infer by looking at the change of slope of λ₂) that can reach several orders of magnitude.

The inset of Figure 5 illustrates the same results normalized by the total strength of the interlayer connectivity matrices. That is, we show S_P/(S_α + S_β), where $Sα=∑Mijα$ and $Sβ=∑Mijβ$⁠, respectively, are the strengths of layers alpha and beta. This allows us to compare networks with the same value of S_P, regardless of whether they are close to a multiplex architecture or to a fully multilayer architecture with no nonzero entries. For p < 〈p*〉, we obtain similar values of λ₂ for all network architectures. It is only for p > 〈p*〉 that the particular structure of the interlayer connectivity matrix begins to play a role in the value of λ₂. We observe that differences start to arise at S_P/(S_α + S_β) ≈ 0.1, which is a relatively small value.

Although we used ER intralayer networks in our above calculations, we obtain similar results for other network models. In particular, the results are qualitatively the same when we construct the intralayer networks using a Barabási–Albert (BA) model (Barabási & Albert, 1999). See the Supplementary Information for details (Buldú & Porter, 2017).

Now that we have analyzed the effects of edge-weight heterogeneity and the number of missing interlayer edges, let’s revisit the multiplex and multilayer networks that we constructed from the MEG recordings. As we saw in Figure 2, both edge-weight heterogeneity and missing interlayer edges occur in our experimental data, and it is thus desirable to investigate how close the experimental networks that we are analyzing are to a transition point p* and how this proximity (or lack thereof) influences the expected value of λ₂.

In Figure 6, we show the values of λ₂ that we obtain for four different network reconstructions based on the MEG data: a homogeneous multiplex network (black circles), a heterogeneous multiplex network (red circles), a homogeneous full multilayer network (orange circles), and a heterogeneous full multilayer network (blue circles). Each network has two layers—one for the alpha band and one for the beta band—and each point corresponds to the mean over the group of 89 individuals. We obtain the heterogeneous versions of the multiplex and multilayer networks by calculating MI values between the brain regions in the frequency bands, as we described previously (also see Materials and Methods). Second, we construct the homogeneous versions of both the multiplex and multilayer networks by assigning the same weight 〈c〉 to all interlayer edges, where 〈c〉 is the mean of the weights of the interlayer edges in their heterogeneous counterparts. Note that, in this case, the homogeneous multiplex and multilayer networks do not correspond to real functional networks; instead, they are reference networks that we use to quantify the consequences on λ₂ of the intrinsic edge heterogeneity and missing edges in real functional networks.

To assess how close the real networks are to the transition point p*, we multiply the values of the interlayer edges by a parameter p, which we increase from p = 0 to a value above 1. We then calculate the strength S_P of the corresponding interlayer connectivity matrix P^αβ to allow a comparison between the multiplex network and the multilayer networks, and we plot the value of λ₂ versus S_P. The red dashed line in Figure 6 corresponds to the theoretical predictions for p < p* (i.e., λ₂ = 2p〈c〉, where 〈c〉 is the mean weighted degree of the nodes in the interlayer connectivity matrix P^αβ). The black and blue dashed lines are, respectively, the value of λ₂ for the aggregated network divided by 2 (i.e., $λ2,agg=12λ2(Lα+Lβ))$ and the value of λ₂ for the unfiltered (monolayer) functional network. The vertical solid line corresponds to the case p = 1 for the mean of the heterogeneous full multilayers networks (i.e., the networks that we obtain by taking into account all interlayer correlations without modifying their weights). Interestingly, this network occupies the region in which the evolution of λ₂ with respect to p changes slope from λ₂ = 2p〈c〉 to λ₂ = λ_2,min(α,β) + λ₂(pC), which suggests that the frequency-based multilayer networks that we are analyzing are close to the transition point p*. As we have seen, it is near this point where the value of λ₂ is most influenced by the effects of edge heterogeneity and missing interlayer edges. It is also worth noting that the four network representations have rather different values of algebraic connectivity (and hence, we expect, have different spectral properties more generally), except when p ≪ p* (see the inset of Figure 6). Although it looks from the plot that the algebraic connectivities of the heterogeneous multilayer networks (blue circles) and those of the homogeneous multiplex networks (black circles) may converge to the same value for large S_P, they eventually cross each other when we further increase S_P (not shown).

In Figure 7, we show (top panel) the PDF of the values of p* for the multilayer networks of each of the 89 individuals and (bottom panel) the percentage of deviation of λ₂ with respect to the value λ_2,agg of the aggregated network. We observe that the peak of the PDF for p* is near p = 1. That is, the multilayer networks that we constructed from empirical data are close to the transition point. Figure 7 also illustrates discrepancies in the values of λ₂ between the aggregated and multilayer networks, indicating that it is necessary to differentiate between the two cases when interpreting the value of λ₂, especially when comparing results from different studies.

Finally, we investigate how the combination of layers from different frequency bands affects the value of λ₂. In our analysis thus far, we have focused on a two-layer network that we constructed using the alpha and beta frequency bands, because they are often associated with brain activity during the resting state. Nevertheless, because the signal has been filtered into four different frequency bands, there are eight possible combinations of two layers. In Figure 8, we show the relation between λ₂ for all possible combinations of two-layer networks versus that for the full four-layer multilayer network. Observe the strong correlation of the two-layer networks that include the gamma layer (especially the one that consists of the theta and gamma layers) with the full multilayer network. One can explain such a correlation by inspecting the total strength of each frequency band. In Table 2, we separate the intralayer and interlayer strengths to facilitate interpretation of the results. We observe that the gamma band is the least-active layer, as it is the one with the smallest intralayer strength. Nevertheless, it has the largest interlayer strength (i.e., the sum of the weights of its interlayer edges with all other layers is the largest), so it is the layer that appears to interact most strongly with the other layers. Because (as we have seen) the full multilayer functional networks are close to the mean 〈p*〉 of the transition point, the weight of the intralayer connections has a strong influence on the value of λ₂. Therefore, the two-layer theta–gamma network, which includes the layers with the largest interlayer strengths, is the one with the strongest correlation of λ₂ with the full multilayer network (see Figure 8C).

Using network analysis as a tool for analyzing brain-imaging data, and (more specifically) implementing and studying a multilayer description of brain activity, has both advantages and drawbacks that must be investigated carefully. As we have discussed, it is possible to encode such imaging data either as a multiplex network or as a more general type of multilayer network, but different choices lead to different results, which one must then interpret from a neuroscientific perspective. In our paper, we have analyzed the implications of such a choice on spectral information—and specifically on the algebraic connectivity λ₂, which has been related to structural, diffusion, and synchronization properties of networks (Almendral & Díaz-Guilera, 2007; Gómez et al., 2013; Masuda et al., 2017; Newman, 2010; Radicchi & Arenas, 2013; Van Mieghem, 2011). We have seen how heterogeneity (of the weights) of the interlayer edges of multiplex networks leads to deviations of the theoretical predictions that one obtains when all interlayer edges have an equal weight p, and we observed that the deviation is even larger when interlayer edges are missing. The importance of these results, which entail large differences in qualitative dynamics, is underscored by the fact that both edge heterogeneity and missing interlayer edges are common features of brain-imaging data.

It is also important to understand that although a multiplex description of brain networks is an important and useful approach to integrate multivariate information, and it is often an extremely natural approach—such as for integrating anatomical and functional networks, examining the temporal evolution of a functional network, and so on—it is just an initial step towards developing increasingly complete models to better analyze the spatial and temporal complexity of brain networks (Papo, Zanin, & Buldú, 2014). In this quest, a natural step—although certainly not the final step—is to use a full multilayer description for subsequent analyses of frequency-based brain networks.³ In this type of network, it is more biologically plausible to represent brain activity as a full multilayer network than as a multiplex network. The reasons are twofold: (i) a brain region does not necessarily coordinate with itself at different frequencies; and (ii) there may also exist CFCs between different brain regions. Consequently, interlayer connections do not necessarily follow the multiplex paradigm, because some of the interlayer edges probably should not be present, whereas others should be included to account for interlayer coupling between different nodes (i.e., different brain regions), thereby necessitating a fully general multilayer approach.⁴ Nevertheless, the use of frequency-based multiplex networks has been an effective strategy in previous work for successfully distinguishing between healthy individuals and individuals who are suffering from schizophrenia (De Domenico, Sasai, & Arenas, 2016), so there are times when a multiplex approach works well.

As we have illustrated with our computations, because the number of interlayer edges scales with N in multiplex networks but with N² in the multilayer networks that we examined, there are important quantitative differences in the values of the algebraic connectivity in the two types of networks. Nevertheless, the effects of edge-weight heterogeneity and missing interlayer edges on full multilayer networks are qualitatively similar to those in multiplex networks. Specifically, the analytical predictions that we obtained with homogeneous multilayer networks deviate from those of heterogeneous cases, showing an important discrepancy for values of the coupling-strength parameter p that are larger than the transition point p*. The effect of missing interlayer edges is particularly dramatic for p > p*.

Our analysis of experimental imaging data confirms our results with synthetic networks. The small number of interlayer edges in a multiplex network leaves the network in a region in which layers behave as if they are structurally independent, and a deeper analysis reveals that it is necessary to increase the weights of the interlayer edges from the experimental data by several orders of magnitude to reach the transition point p*. Consequently, using a multiplex representation leads to a network with structurally independent layers, unless one uses extremely large weights (i.e., p ≫ 1) for the interlayer edges. However, the weights of the interlayer edges cannot be increased artificially without limits, because eventually they lose their biophysical meaning in attempts to compare them to intralayer edge weights.

Interestingly, our results changed drastically when we used a full multilayer description, in which a brain region whose dynamics are in a given frequency band can now couple with other brain regions in a different frequency band. When we considered all possible interlayer edges, the majority of the 89 individuals in the experiment are close to the transition point p*, even though the percentage of weight in the interlayer edges is always about 10% of that of the intralayer ones (see Table 2). As we showed in our analysis of synthetic data, it is around this point that edge-weight heterogeneity and missing interlayer edges begin to yield important discrepancies from the theoretical results of homogeneous multilayer networks (see Figure 4A).

It is important to note that the percentage of total weight of interlayer edges (relative to the total weight of all edges) depends strongly on how one measures coordination between brain regions. When one uses diagnostics other than mutual information, one calculates different values for the percentages, which can yield multilayer networks either above or below the transition point p*. Thus, no matter what measure one uses, it is mandatory to first calculate the percentage of total edge weight that arises from interlayer edges and to interpret the value of λ₂ with respect to the value of p*.

The fact that, in frequency-based multilayer networks, the weights of interlayer edges strongly influence the value of λ₂ despite being much smaller than the weights of intralayer edges (see Table 2) highlights the importance of adequately evaluating CFCs, which traditionally have been disregarded when representing brain activity in terms of functional networks. As has been discussed prominently in neuroscience, including in critiques of connectomics (Kopell, Gritton, Whittington, & Kramer, 2014), the dynamics matter, and investigations of dynamics must include incorporation and analysis of coupling between different frequency bands. The methodology that is used to quantify interactions between brain regions at different frequencies leads to different values of λ₂; and there are also other important dynamical issues, such as phase–amplitude correlations, that we have not investigated in this paper. (See Bastos & Schoffelen, 2016 for a review of how common reference, volume conduction, field spread, or common input affect the quantification of coordinated activity between brain regions.) These caveats notwithstanding, our analysis illustrates a viable approach for examining the effects of such phenomena on λ₂ (and hence on spectral structure) in multilayer functional networks, and we expect that similar qualitative phenomena will arise in both multiplex and full multilayer networks that are constructed using other choices (e.g., different measures than MI) from the ones that we chose in order to provide a concrete illustration.

We initially studied two-layer networks of the alpha and beta frequency bands, because (i) it is known that these bands incorporate a large portion of the power spectrum of brain activity during resting state and (ii) these frequency bands typically exhibit stronger synchronization between brain regions (e.g., see Brookes et al., 2016) than the others. The alpha and beta bands are thus the most commonly studied frequency bands in resting-state studies, and we followed this tradition. Nevertheless, in our particular case, our comparison between two-layer and four-layer networks illustrates that the gamma band is the one that most influences the spectral properties of the full four-layer multilayer network (see Figure 8). As we showed in Table 2, the theta and gamma layers are the ones whose interlayer edges give the largest contributions, thereby leading to the strongest correlation between a two-layer network (with theta and gamma layers) and the complete (four-layer) multilayer network. This fact highlights the importance of the well-known phase–amplitude correlations between the theta and gamma frequency bands (Aru et al., 2015; Canolty et al., 2006), as the former acts as a carrier of fast-amplitude fluctuations in the latter. Consequently, theta–gamma coupling may be of fundamental importance for understanding the multilayer nature of functional brain networks.

It also worth mentioning that the facts that we are constructing frequency-based multilayer networks from (i) data at the sensor level and (ii) using MI as a measure to quantify mutual interdependency between brain regions leads to unavoidable errors in the quantification of edge weights, due to signal mixing and spurious correlations from common sources. However, to date, there does not exist an error-free methodology to construct functional brain networks, and other alternatives (such as source reconstruction or using different measures to evaluate interdependency between brain regions) have their own drawbacks. For a detailed discussion, see the subsection Methodological Considerations (in Materials and Methods).

Finally, although our calculations with experimental data used resting-state MEG recordings, we expect to observe similar behavior in frequency-based multilayer networks of different origins—whether obtained from any of a large variety of different cognitive or motor tasks, with different brain-imaging techniques, or even if they come from a completely different system (such as functional climate networks; Tsonis & Roebber, 2004). However, despite the generality of our results, it is important to examine richer models of frequency-based brain network, such as ones that include spatial constraints and temporal evolution.

The data were made available by the Human Connectome Project (HCP); see http://www.humanconnectome.org/ and Larson-Prior et al. (2013) for details. The experimental data consist of magnetoencephalographic (MEG) recordings of a group of 89 individuals, during resting state, for a period of approximately 2 min. During the scan, subjects were supine and maintained fixation on a projected red crosshair on a dark background. Brain activity was scanned on a whole-head MAGNES 3600 (4D Neuroimaging, San Diego, CA, USA) system housed in a magnetically-shielded room, and it included up to 248 magnetometer channels. The root-mean-squared noise of the magnetometers is about 5 fT/sqrt (Hz) on average in the white-noise range (above 2 Hz). The data were recorded at a sampling rate of f_s ≈ 508.63 Hz. Five current coils attached to a subject, in combination with structural-imaging data and head-surface tracings, were used to localize the brain in geometric relation to the magnetometers and to monitor and partially correct for head movement during MEG acquisition. Artifacts, bad channels, and bad segments were identified and removed from the MEG recordings, which were processed with a pipeline based on independent component analysis to identify and clean environmental and subject artifacts (Larson-Prior et al., 2013). After this process, the number N of channels considered for each individual was in the range 235 ≤ N ≤ 246 (with a mean of 243.42), because some channels were used as references and others were disregarded.

An important issue is the applicability of our results to functional brain networks that are constructed using different approaches, such as by using source reconstruction or with different ways of evaluating coordination between brain regions. In our study, we constructed functional brain networks from magnetoencephalographic recordings, which are measured at the sensor level, instead of using the actual magnetic field that is generated in each brain region (i.e., at a source). (One can infer such a field using source-reconstruction methods.) Working with time series at the sensor level necessarily entails signal mixing (Schoffelen & Gross, 2009), whereas the use of source reconstruction has what is known as the “inverse problem” (of inferring the actual magnetic field that is created by the brain regions), and it thus requires the introduction of a priori assumptions in the model that one uses for source reconstruction. This issue has led to a diversity of algorithms to obtain source-reconstructed time series and, in many cases, different algorithms yield qualitatively different time series at the source level (see, e.g., Schoffelen & Gross, 2009; Belardinelli, Ortiz, Barnes, Noppeney, & Preissl, 2012). Consequently, there is an open debate about what constitutes the most appropriate methodology to construct functional networks using source reconstruction (Palva & Palva, 2012; van Diessen et al., 2015).

Regardless of how one constructs a functional brain network (and whether one uses sensors or sources), our analysis has the same qualitative implications. Specifically, estimating the value of λ₂ from matrices associated with function brain networks is affected strongly by edge-weight heterogeneity and missing interlayer edges. One also observes the same qualitative differences between a multiplex construction and a full multilayer one, as this phenomenon depends on the number of interlayer edges and not on specific details of the construction of functional brain networks.

Similar reasoning applies to the construction of frequency-based brain networks using different ways of quantifying coordination between brain regions. We used MI because it is able to successfully capture both linear and nonlinear correlations, it is an adequate method for quantifying interdependencies between signals that are split into frequency bands (David, Cosmelli, & Friston, 2004), and it has been used widely in the past to construct functional brain networks (see, e.g., Bassett et al., 2009; Becker et al., 2012; Chai, Walther, Beck, & Fei-Fei, 2009; Deuker et al., 2009; Jin, Seol, Kim, & Chung, 2011). On the down side, MI cannot avoid zero-lag correlations that originate from common sources. (This flaw also occurs in many other commonly-employed similarity measures, such as partial correlation, coherence, and phase-locking value (Bastos & Schoffelen, 2016).) The systematic deletion of zero-lag correlations entails a very strict assumption, but because not all zero-lag correlations are due to common sources or volume conduction (Christodoulakis et al., 2013; Porz, Kiel, & Lehnertz, 2014; Vicente, Gollo, Mirasso, Fischer, & Pipa, 2008), there does not exist a consensus on how to eradicate them (see, e.g., Brookes, Woolrich, & Barnes, 2012; Stam, Nolte, & Daffertshofer, 2007). It is thus important for future work to conduct a systematic comparison between the most widespread coordination measures and their consequences on the spectral properties of functional networks. As we stated previously, the effect that edge-weight heterogeneity and missing interlayer edges have on functional multiplex and multilayer networks is a general phenomenon, as it does not rely on the specific similarity measure that one employs to evaluate coordination between brain regions, as long as there is sufficient heterogeneity in the weights of the interlayer edges. As shown in Brookes et al. (2016), such heterogeneity is expected in frequency-based functional brain networks, given the complicated interactions between different brain regions at different frequencies.

Javier M. Buldú: Conceptualization; Data curation; Formal analysis; Methodology; Writing — original draft; Writing — revisions. Mason Alexander Porter: Conceptualization; Methodology; Writing — original draft; Writing — revisions.

We thank John Allen for fruitful conversations. We also thank the referees for their helpful suggestions.

J.M.B. acknowledges financial support from Spanish MINECO (project FIS2013-41057) and from Salvador de Madariaga Program (PRX15/00107), which allowed him to visit University of Oxford in summer 2016.

Javier M. Buldú, Ministerio de Economía y Competitividad (http://dx.doi.org/10.13039/501100003329), Award IDs: FIS2013-41057 and FIS2017-84151-P. Javier M. Buldú, Ministerio de Economía y Competitividad (http://dx.doi.org/10.13039/501100003329), Award ID: PRX15/00107.

Multilayer brain network:

A brain network with more than one layer. Layers are connected to each other through interlayer edges, which link node-layers from different layers.

Node-layer:

Within each layer, a node-layer represents the dynamics recorded by a given magnetometer that is filtered at a specific frequency band. Thus, each magnetometer has an associated node-layer on each layer.

Layer:

A portion of a multilayer network with a particular set of connections (called “intralayer edges”) between nodes. In our paper, each layer represents a specific frequency band of a functional brain network.

Multiplex brain network:

A particular kind of multilayer brain network in which interlayer edges occur only between counterpart nodes (which, in this paper, represent the same brain region) in different layers.

Magnetoencephalography (MEG):

Brain-imaging technique that measures the magnetic field of the brain using a set of magnetometers placed on the head of an individual.

Algebraic connectivity:

The second-smallest eigenvalue (λ₂) of the combinatorial Laplacian (or supra-Laplacian) matrix.

Combinatorial supra-Laplacian matrix:

The combinatorial Laplacian matrix of a supra-adjacency matrix.

Cross-frequency coupling (CFC):

In the brain, CFC is the mechanism through which dynamics at a certain frequency affect, or are affected by, the dynamics at other frequencies.

Combinatorial Laplacian matrix:

Given the adjacency matrix A of a network and the degree matrix D, which consists of the degree of the nodes along the diagonal and 0s in all other entries, the combinatorial Laplacian matrix is defined as L = D − A.

Mutual information (MI):

A measure to evaluate statistical interdependencies between a pair variables, X and Y, based on comparing a joint distribution p(X, Y) with the product p(X)p(Y) of marginal distributions.

Supra-adjacency matrix:

An adjacency matrix associated with a multilayer network; it codifies how much different node-layers interact with each other.