Community detection in multi-frequency EEG networks

Functional connectivity networks of the human brain are commonly studied using tools from complex network theory. Existing methods focus on functional connectivity within a single frequency band. However, it is well-known that higher order brain functions rely on the integration of information across oscillations at different frequencies. Therefore, there is a need to study these cross-frequency interactions. In this paper, we use multilayer networks to model functional connectivity across multiple frequencies, where each layer corresponds to a different frequency band. We then introduce the multilayer modularity metric to develop a multilayer community detection algorithm. The proposed approach is applied to electroencephalogram (EEG) data collected during a study of error monitoring in the human brain. The differences between the community structures within and across different frequency bands for two response types, i.e. error and correct, are studied. The results indicate that following an error response, the brain organizes itself to form communities across frequencies, in particular between theta and gamma bands while a similar cross-frequency community formation is not observed following the correct response.


I. INTRODUCTION
The human brain consists of various units which interact with each other through structural and functional links.Advances in functional and structural neuroimaging technology allow the brain to be modeled as a network using graph theoretic tools.In this modeling, the nodes correspond to the different brain units and the edges represent structural or functional connections among these units [1].In order to characterize the topology and dynamics of brain networks, various descriptive and inferential network measures are utilized; such as centrality, degree distribution and small-worldness [2]- [6] with respect to disease, task, learning, cognitive control, attention and memory [1], [4], [6]- [11].
Current network models have been mostly limited to examining a single network instance either of a subject, a frequency band or a task.However, most neurophysiological recordings, such as the electroencephalogram (EEG), allows one to capture brain dynamics across multiple temporal and spatial scales.Reducing this rich information into a single network disregards the high amount of dependency that exists between networks of different subjects, frequency bands and time points.Thus, a principled mathematical framework to accurately study this multiplicity of brain connectivity is needed.
Recently, multilayer networks [12]- [14] have been proposed as a mathematical framework to study multiple networks simultaneously.Multilayer networks consist of multiple layers, each of which carry information from a different network while the interactions between layers represent the dependency between these networks.Due to their ability to represent and study multi-dimensional and multi-scale data, multilayer networks gained attention in network neuroscience [15]- [17].To this end, measures that are developed to characterize the layered structure of multilayer networks are employed to study brain connectivity, allowing for a richer analysis than prior single network based studies [18].
In this paper, we aim to characterize the topological organization of multilayer brain networks by developing a multilayer community detection algorithm.The proposed approach is outlined in Fig. 1.The main contributions of the proposed framework are: • First, we introduce a signal processing based approach for constructing multi-frequency networks from EEG data.In this approach, the intra-and inter-layer edges are quantified by novel time-frequency based phase synchrony and phase amplitude coupling (PAC) measures, respectively.Thus, the constructed network is a general multilayer network with inter-layer edges allowed between all brain regions.
• Second, modularity function is defined for multilayer networks based on a multilayer null model that preserves the layer-wise node strengths while randomizing the remaining characteristics of the network.The corresponding modularity maximization algorithm is capable of detecting communities that are both common across a subset of layers and private to each layer.
• Third, once the multilayer community structure for each subject is detected, a consensus clustering approach based on multiview modeling is used to obtain group community structure of a set of subjects.In the multiview network model, each layer corresponds to a subject's coclustering matrix constructed from community structures obtained from multiple runs of modularity maximization.Multilayer spectral clustering is then applied to the re- sulting multiview network.This approach considers all of the subjects simultaneously and detects a common community structure across subjects summarizing their network topology.
• Finally, this is the first work of its kind that considers a general multilayer network perspective for EEG data and uses the community structure to characterize the brain activity across multiple frequencies, simultaneously.
In particular, the group level differences between the two response types during Flanker task, i.e., error and correct, are evaluated from a multi-frequency network perspective.

II. RELATED WORK A. Multilayer Brain Networks
Recently, the study of brain networks has undergone a process of adapting classical single-layer functional connectivity network concepts to a more general multilayer description [15], [17]- [20].The meaning of layer can vary depending on context; including multiple modalities, subjects, time points, and frequency bands.Some examples include multi-modal networks where each layer represents a different modality.For example, [21] study a two layer network that combines structural and functional networks of a subject, which are constructed from the subject's diffusion tensor imaging (DTI) data and funtional MRI (fMRI), respectively.Two layers are connected with inter-layer edges, which link a brain region in the functional network to itself in the structural network.Another example is the temporal network where each layer represents connectivity over some time period [5], [7].The inter-layer edges are usually allowed only between consecutive time periods and between nodes that represent the same brain regions.Temporal brain networks are employed to analyze how the community structure and the network topology of functional networks change over time.
More recently, functional multilayer networks are constructed, where layers correspond to well-known frequency bands at which the brain operates [20], [22]- [25].Compared to single-layer networks, which do not distinguish the contributions coming from different frequency bands and consider a single frequency band, multilayer networks allow the study of connectivity across multiple frequency bands simultaneously.For instance, recent work has considered multi-frequency functional connectivity networks for resting state fMRI [26]- [29].These studies have shown that the network organization such as hubs may differ across different frequency bands [28].More recently, multi-frequency networks were considered in neurophysiological studies [20], [22], [30], [31].Brookes et al. [22] use magnetoencephalographic (MEG) recordings to construct functional multilayer networks, where each layer corresponds to the links within a given frequency band, and the inter-layer edges correspond to the cross-frequency coupling across frequency bands resulting in a general multilayer network.It is shown that there is statistically significant difference between supra-adjacency matrices of the multilayer networks of control and schizophrenia subjects.Similar multi-frequency networks are constructed from resting state MEG recordings of subjects affected by Alzheimer's disease.Analyzes of these networks shows the abnormal distribution of regional connectivity across frequency bands in unhealthy subjects compared to healthy individuals, revealing an abnormal loss of interfrequency centrality in memory-related association areas [32].More recently, a multi-frequency network model was considered for EEG data [31], [33], where a multilayer network is constructed by concatenating the functional connectivity matrices for 40 different frequency bins.Generalized Louvain algorithm for multiview networks is then employed to detect the community structure of these multi-frequency networks.However, in this work, the inter-layer connections were only allowed between nodes corresponding to the same brain areas, implying that the network is not a general multilayer network.

B. Community Detection in Multilayer Networks
Community detection is one of the main network analysis tools and has been applied to several complex systems, including biological and social [34].In single-layer networks, communities are defined as groups of nodes that are more strongly connected among themselves than they are to the rest of the nodes, while in multilayer networks, these groups of nodes are shared across multiple layers.Traditional community detection methods usually need modifications to achieve a good performance on a multilayer network.
In [36], the definition of modularity is extended for multilayer networks.The proposed modularity index is then maximized using Girvan-Newman and Louvain algorithms.In [38], [40], the normalized cut measure is extended to multiview and multilayer networks by constructing a block Laplacian matrix where each block corresponds to a layer.Then they apply spectral clustering to this matrix for finding the communities.In [48], a Weighted Stochastic Block Model (WSBM) is proposed to detect communities on each layer and also the ones that are shared across layers of a heterogeneous weighted network .However, this method only detects common communities if they are shared by all of the layers.In [44], a method based on random walks, Multiplex-Infomap, is proposed.Although this method identifies intra-and inter-layer communities in multilayer networks, it usually groups physical nodes across layers in the same community.In [46], the authors propose an approach based on nonnegative matrix factorization that detects intra-and inter-layer communities in general multilayer networks.However, this method is only applicable to networks with two layers.In [45], a community detection method based on Label Propagation Algorithm (LPA) is proposed for multilayer networks.LPA follows the intuition that a label would get trapped inside a group of densely connected nodes and end up having the same label indicating they belong to the same community.This method simultaneously identifies the communities and the subset of layers where these communities belong to.Although these methods are developed networks with multiple layers, only a few of them are applicable to general multilayer networks where connections may exist between different physical nodes across layers [36], [40], [46].

III. BACKGROUND A. Multilayer Networks
An undirected single-layer network G = (V, E, A) is represented by a node set V with |V | = N , an edge set E ⊆ V × V , where e uv are associated with weights w uv and adjacency matrix A ∈ R N ×N with A uv = A vu = w uv if e uv ∈ E and 0, otherwise.The strength of node u is the total weight of the edges incident to it, i.e. s u = n v=1 A uv .An undirected multilayer network is a quadruplet M = (V, L, V, E) where V is the set of physical entities, L is the set of layers with |L| = L [12].V ⊆ V × L with |V | = N is the set of nodes, which are representations of physical entities in layers and E ⊆ V × V is the edge set.In this paper, nodes of a multilayer network are represented as u h , where u ∈ V and h ∈ L.An edge between u h and v k is denoted by e hk uv and associated with the weight w hk uv .V can be partitioned into layers, that is , where E h is the intra-layer edges of layer h and E hk is the inter-layer edges between nodes in layers h and k .Using this notation, one can define intralayer graphs G h = (V h , E h ) and bipartite inter-layer graphs where A h is the adjacency matrix of G h and A hk is the incidence matrix of the bipartite graph, G hk .Layer-wise strength of a node u h is the sum of weights corresponding to the edges connected to the nodes in layer k , i.e., s k Multilayer graphs can be categorized based on the structure of their inter-layer graphs.In this paper, we consider two types of multilayer graphs.The first one is multiview graphs, where inter-layer edges are only allowed between nodes representing the same physical entities.The second type is general multilayer graphs, which do not have any restrictions on the structure of G hk 's.In the following, general multilayer graphs will be referred to as multilayer graphs for simplicity.

B. Modularity
Given an undirected single-layer graph G = (V, E, A), community detection aims to partition the node set This is usually achieved by optimizing a quality function that quantifies the goodness of a given partition [34], [49].Among various objective functions, the most commonly used one is modularity [50], which quantifies the quality of a partition by comparing the intracommunity edge density to that expected under a null model and is calculated as follows: where P ij is the expected edge weight between nodes i and j under the null model, g i is the community of node i, and δ gigj = 1 if g i = g j and 0, otherwise.Different null models can be used to define P ij depending on the graph under study.The most commonly used null models are the configuration and Erdős-Rényi null models [51].Modularity as defined in (2) suffers from the resolution limit [52], i.e., it cannot detect small communities.Therefore, it has been extended to include a resolution parameter γ [51]: By tuning γ, one changes the resolution of the modularity function such that larger γ values can detect smaller communities.

A. Data Acquisition
The EEG data was acquired during a cognitive controlrelated error processing task where the subjects performed a letter version of the speeded reaction Flanker task [53].The experimental protocol of this study was approved by the Institutional Review Board (IRB) of the Michigan State University (IRB: LEGACY13-144).The data collection was conducted by following the regulations approved by this protocol.Prior to data acquisition, all subjects signed an informed and written consent form.During recording, each subject was presented with a string of five letters at each trial.Letters could be congruent (e.g., SSSSS) or incongruent stimuli (e.g., SSTSS) and the subject was instructed to respond to the center letter with a standard mouse.The trials started with a flanking stimulus (e.g., SS SS) of 35ms followed by the target stimuli (e.g., SSSSS/SSTSS) displayed for about 100 ms.The total display time is 135 ms, followed by a 1200 to 1700 ms inter-trial break between the trials.These trials capture the Error-Related Negativity (ERN) after an error response and the Correct-Related Negativity (CRN) after a correct response.As earlier studies suggested a rise in synchronization related to ERN for the 25-75 ms time window [54], all the analysis was conducted for the 25-75 ms time period of the data.For each subject, 480 total trials (each of 1-second in duration) were recorded, where the number of error trials varied from 20 to 61 across the subjects.For a fair comparison between ERN/CRN, the same number of correct trials were selected randomly.The EEG signals were recorded with a BioSemi ActiveTwo system using a cap with 64 Ag-AgCl electrodes placed at standard locations of the International 10-20 system.The sampling rate of the data was 512 Hz.After cleaning the artifacts, volume conduction was minimized using the Current Source Density (CSD) Toolbox [55].A multilayer network with four layers is constructed for each subject and each response type (error vs. correct) where layers correspond to the four EEG frequency bands: θ (4-7 Hz), α (8-12 Hz), β (13-30 Hz), γ (31-100 Hz).In this paper, we consider data from 20 participants.

B. Intra-layer Edges
The intra-layer edges are computed using a reduced interference Rihaczek (RID-Rihaczek) time-frequency distributionbased phase synchrony index known as RID Rihaczek timefrequency phase synchrony (RID-TFPS) [56], [57].Compared to existing measures like wavelet transform, RID-TFPS provides uniformly high time and frequency resolution and has been reported as a robust phase synchrony index for noisy signals [57].For an arbitrary signal x(t), the RID-Rihaczek distribution is obtained as [57]: where e −(θτ ) 2 /σ refers to the Choi-Williams kernel and e jθτ /2 defines the kernel function for the Rihaczek distribution [58].C(t, f ) belongs to Cohen's class of distribution and thus satisfies the marginals and preserves the energy.This complex time-frequency distribution is utilized to calculate the phase difference φ u,v (t, f ), between two signals x u and x v as: The phase synchrony between two brain locations measures the consistency of the phase differences φ u,v (t, f ) across trials and is quantified using the Phase Locking Value (PLV) index.PLV measures the degree of phase locking by computing intertrial variability of the phase differences as follows [59]: where K is the total number of trials and φ k u,v (t, f ) is the phase difference between x k u and x k v as given in (5) for trial k.After computing the pairwise PLV values between all possible brain locations, the average pairwise synchrony within a predefined time window of interest, W = [t 1 , t 2 ], and a chosen frequency band is used as intra-layer edge weights, i.e., where N is the number of brain locations and and |h| is the bandwidth of the particular frequency band h.

C. Inter-layer Edges
The inter-layer edges are calculated using a form of crossfrequency coupling known as PAC, which computes the modulation of the amplitude/power of a high frequency rhythm along the phase of a slower frequency rhythm [60], [61].In this study, PAC between two brain locations is estimated by applying a RID-Rihaczek time-frequency-based PAC approach [62], [63].To quantify PAC, first we obtain analytic signals C u (t, f ) and C v (t, f ) at node u and v, respectively, by using RID-Rihaczek distribution as defined in (4).These analytic signals are utilized to extract the instantaneous high frequency amplitude envelope a u fa (t) at node u, and the instantaneous low frequency phase φ v fp (t) at node v, where f p and f a are respective frequencies within the hth and k th frequency bands.
For node u, a u fa (t) is obtained from the frequency constrained time marginal of C u (t, f ) as: where f a1 and f a2 is the bandwidth around the chosen high frequency.Similarly, the low frequency phase components at node v is obtained from C v (t, f ), as: Once the amplitude and phase components are extracted, PAC is estimated by distributing a u fa (t) and φ v fp (t) to a composite vector in the complex plane at each time point and measuring the amplitude normalized modulation index (MI) proposed in [64] as follows: As the amplitude of the signal normalizes the computed MI, the MI value is scaled between 0 and 1 [65].Thus, the interlayer edges between node u and v are constructed as

A. Multilayer Modularity
Modularity function compares intra-community edges of the observed network to those of a null model to quantify the quality of a community structure.The null model is a random graph model with some properties, e.g.edge density, of the observed network.The properties to be preserved are determined such that they are assumed to be unrelated to the community [66].For example in the configuration null model, the degree of each node is preserved so that the identified community structure would not be affected by the heterogeneity of the degree distribution.This assumption is based on the fact that nodes with a high degree tend to connect with each other merely because they have high number of connections and not necessarily because they are within the same community [67].To prevent this tendency to bias community detection, the null model preserves the node degrees.On the other hand, Erdős-Rényi null model does not make such an assumption and allows the identified community structure to be influenced by the degree distribution.
Based on this insight on the role of null models, we extend the definition of modularity function to multilayer networks by considering which properties of the observed multilayer network we want to preserve in the null model.In Fig. 2a, the supra-adjacency matrices of the multilayer EEG networks are shown for one subject for both error and correct responses.From this figure, it can be seen that the edge weights are heterogeneous across layers for both response types.For instance, the intra-layer edges for θ band and interlayer edges between θ and γ bands are very strong, while intra-layer connections for α, β and γ bands are weaker for error response.Any community detection method applied to this network without taking this heterogeneity into account would likely partition the nodes based on the layer label rather than the true community membership.In order to prevent this, the null model used in the definition of the modularity function should preserve the heterogeneity of edge weights across layers as follows.
Multilayer Configuration Null Model: Let M be a multilayer network with m hh defined as the total weight of the intralayer edges in layer h and m hk defined as the total weight of the inter-layer edges between layers h and k .The multilayer configuration null model preserves layer-wise node strengths while randomizing the remaining characteristics of M. The expected edge weight between u h and v k is then: where δ hk = 1 if h = k and 0, otherwise.The multilayer modularity is then defined as follows: where γ hk 's are the resolution parameters corresponding to intra-and inter-layer graphs and ω hk 's are the scaling parameters that weigh the importance of different intra-and inter-layer connections.In the following, we use a common resolution parameter across layers, i.e. γ hk = γ.Scaling parameters are set as ω hk = ω if h = k and 1, otherwise.( 11) can be optimized with greedy algorithms, such as the Louvain algorithm [68], developed for maximizing the singlelayer modularity function defined in (2).In this work, we used the Leiden algorithm, which is an extension of the Louvain algorithm with better performance [69].
Since the proposed multilayer configuration null model preserves the heterogeneity of edge weights across layers, the corresponding modularity matrix should not have the heterogeneity observed for the original supra-adjacency matrices in Fig. 2a.Fig. 2c.shows the modularity matrices, where the heterogeneity across layers is reduced.On the other hand, Fig. 2b.shows the modularity matrices computed with the singlelayer configuration null model.In this case, the heterogeneity across layers is very similar to that of the supra-adjacency matrix.These results indicate that multilayer configuration null model defined in (10) preserves the heterogeneity of edge weights across layers as desired.
While recent work [36] has considered a similar definition of modularity for multilayer networks, in that work, the resolution parameter is set to 1 and the inter-layer scale parameters are determined from m hk .However, the resolution limit of modularity function necessitates one to consider different resolution parameters to obtain better community structures [70].Moreover, as our results show, considering different values for inter-layer scale can reveal important aspects of multilayer community structure.Therefore, in this work, we do not fix these parameters and try to learn their optimal values as described in Results section.

B. Group Community Structure
Once the community structures of the multilayer networks for a group of subjects are detected, it is often desirable to find a group community structure, which summarizes the shared communities across subjects.There are two widely used approaches in network neuroscience to find group community structure [19].The first approach is based on selecting a representative subject and using this subject's community structure as the group structure [71].Representative subject is selected as the one whose community structure has the greatest average similarity to the other subjects' community structures.Since this approach ignores the community structure of all but one subject, it generally does not yield satisfactory results.The second approach is consensus clustering [72], which combines information from multiple community structures to derive the group community structure.This approach constructs a coclustering matrix A where A uv is the number of times nodes u and v are in the same community across all subjects and applies community detection to A to find the group structure.However, this method is highly dependent on the quality of the detected community structure for each subject.Since modularity maximization is an NP-hard problem [73], modularity maximization algorithms yield locally optimal results.Furthermore, modularity function has a strong degeneracy [74], i.e., multiple community structures with similar modularity values may be significantly different from each other.Therefore, there can be a diverse set of informative community structures for each subject and relying on a single community structure from each subject can be problematic for consensus clustering [75].
In order to address these issues, in this paper we propose a group community structure detection method based on multiview graphs.Given L subjects, for each subject we maximize the modularity function with the optimal γ and ω values 100 times to obtain 100 community structures.From these community structures, we construct a co-clustering matrix A h for each subject h ∈ {1, 2, . . ., L} as described above.This process results in L co-clustering matrices, which can be modeled as the layers of a multiview graph, where each layer is an undirected, weighted graph corresponding to each subject.The group community structure can be found from this multiview graph.
Multiview community detection methods generally fall into two categories [35].The first category of methods find a single common community structure using information across all layers.The second category of methods finds a community structure for each layer while regularizing these structures to have some similarity.Since our goal is to find the group community structure of L subjects, we focus on the first category.In particular, we use Spectral Clustering on Multi-Layer graphs (SC-ML) [76], which finds the multiview community structure by applying spectral clustering to a modified Laplacian defined as: where L h is the the normalized graph Laplacian defined as , where D h is the diagonal matrix of node strengths and U h is the low-rank subspace embedding of layer h.In this work, we set α = 0.5, following the guidelines in [76].

A. Resolution Parameter and Inter-layer Scale Selection
We propose a statistical testing approach comparing the modularity value of the observed multilayer network to that of surrogate networks to determine the resolution and interlayer scale parameters in (11).Such methods have been previously used in network neuroscience literature for selecting the resolution parameter in single-layer networks and coupling parameter in multiview graphs [19], [70], [77].Given an observed multilayer network M and c surrogate multilayer networks, we perform community detection on surrogate multilayer networks for a given pair of γ and ω values, and calculate the modularity values of the detected community structures.Let Q surr be the average of these c modularity values.Next, we perform modularity maximization for M c times and compute the average of the modularity values for the c community structures, Q obs .This process is repeated for different pairs of γ and ω values, and the pair with the largest difference, Q obs − Q surr , is selected as the optimal parameter values.
The surrogate networks are generally constructed from the observed graph by randomly swapping the edges while preserving some properties, e.g.node strength, of the observed graph [78], [79].Since the multilayer EEG networks are fully connected and weighted, we focus on randomization techniques presented in [79], where two approaches are proposed for fully connected single-layer networks.The first technique swaps edges while preserving the edge weights.The second approach randomly modifies the edge weights while preserving node strengths.Both edge weights and node strengths cannot be preserved at the same time when randomizing fully connected weighted networks, as this "in most cases only leaves one possible surrogate network, namely, the original network" [79].We employed the first technique, as the second one may result in surrogate networks with unrealistic edge weights for EEG networks, e.g.edge weights larger than 1.Thus, surrogate multilayer networks are generated as follows: we select two edges e hk uv and e lm st and swap their edge weights.Edges are selected such that h = l and k = m, which ensures that the heterogeneity of edge weights across layers is preserved.
Figs. 3a. and 3b.show the average of Q obs − Q surr across subjects for error and correct responses, respectively.For both response types, optimal γ is found to be close to 0.99, while optimal ω values are observed to be more diverse across subjects, ranging between 0.0-0.2 for error and between 0.0-0.1 for correct.In Fig. 3c, we plotted the histogram of the optimal ω values across subjects.This figure shows that the optimal ω values are non-zero for all subjects except one for the error response.On the other hand, for correct response, the optimal ω values are close to 0 for 7 subjects, while most of the remaining subjects have optimal ω values close to 0. This result indicates that inter-layer connections are not important for the community structure of correct response, as ω = 0 means inter-layer part of the modularity function is equal to zero.For error response, inter-layer edges are influential in community formation indicating the importance of crossfrequency coupling.This result is in line with prior work showing increased cross-frequency coupling between highfrequency oscillations related to motor activity and to visual processing in the gamma band and low frequency cognitive control signals which are activated after an error response [62].It has been suggested that low frequency network oscillations in prefrontal cortex, e.g.θ band, guide the expression of motorrelated activity in action planning [80].

B. Consistency of Community Structures for Error and Correct
After obtaining the community structure for each subject and both response types at the respective optimal γ and ω values, we assess the consistency of community structures across subjects within each response type.We employed the multiview graph constructed in Section V-B, where we quantified the consistency of community structures across subjects by computing the similarity between the layers.In particular, we used Jensen-Shannon (JS) distance, which quantifies the similarity of layers based on their spectral properties [81].To measure the similarity of layers h and k in the multiview graph, where each layer represents the co-clustering matrix of a subject, we let L h and L k be the respective combinatorial Laplacian matrices of layers h and k and define L h = L h /tr(L h ).JS distance is the square-root of JS divergence, which is defined as [82]: where JS distance takes values between 0 and 1.We measured JS distance between each pair of subjects within each response type and plotted the average distance of each subject to other subjects in Fig. 4.This plot shows that the average divergence for each subject with respect to the other subjects is lower for error response compared to the correct response.This indicates that there is more group level consistency in terms of topological organization for the error response.This is in line with prior work [54] that shows that the organization of the functional connectivity networks for correct response is similar to pre-stimulus networks.Thus, there is more variation across subjects for the correct response compared to response-evoked networks following an error response.

C. Group Community Structure for Error and Correct Responses
Once the optimal community structures are obtained for each subject and for each response type, the group community structure is extracted by using SC-ML described above.The number of communities is determined as the average of the number of communities detected for each subject.These values are 5 and 9 for the error and correct responses, respectively.Fig. 5 illustrates the group community structure for error and correct responses for the multi-frequency networks.For error response, it can be seen that there is consistency between the community structure of the low frequency bands, e.g., θ and α.We also note that during ERN there is a community comprised of the frontal-central nodes in the θ-band.This community structure is consistent with prior work that indicates the role of medial prefrontal cortex (mPFC) during ERN in the θband [54].However, the proposed approach shows that this is an across-frequency community as the frontal-central regions from the low-frequency bands are in the same community as the parietal-occipital regions from the γ-band.While this is a new finding in terms of network organization, our recent work indicates that the phase of the θ band oscillations from the frontal-central regions modulate the amplitude of the γ band oscillations in the parietal-occipital regions following an error response [83].Prior studies also hypothesized that error-related negativity initiates the medial frontal based topdown control mechanisms to improve the performance after an error response [84].Thus, the communities detected are consistent with previous literature reflecting higher thetagamma coupling in the medial frontal cortex and relating this with error-related negativity.
The community structure for the correct response is mostly within-layer indicating the lack of coupling across different frequency bands.Our prior work comparing PAC between response types supports this observation as there is significantly higher cross-frequency coupling during error monitoring [62].
The similarity of the group structures for error and correct responses was quantified using three different metrics: Normalized Mutual Information (NMI), Adjusted Rand Index (ARI) and Variation of Information (VI).NMI and ARI are normalized and a value close to 1 indicates similarity, while for VI small values indicate similarity.The similarity is computed for each frequency band.From Table I, it can be seen that the community structure in the γ band has the highest similarity between the two response types while the community structures in the θ band are not similar at all.This is due to the fact that during error monitoring, there is increased phase synchronization in the θ band.Thus, the community structure following an error response is different from the correct response in this frequency band.

VII. CONCLUSIONS
This paper introduced a multilayer model of functional connectivity of the brain.In particular, we provided a datadriven way to construct multi-frequency connectivity networks where layers correspond to different frequency bands.The resulting networks capture both within frequency connectivity and cross-frequency coupling in a single framework.We then introduced a new definition of modularity for multilayer networks such that the null model preserves the heterogeneity of edge weights across layers.The community detection algorithm resulting from the maximization of this new multilayer modularity function is applied to EEG data collected during error monitoring.The results indicate that following an error response, the brain organizes itself to form communities across frequencies, in particular between theta and gamma bands.This cross-frequency community formation is not observed for the correct response which indicates that the cross-frequency coupling is primarily associated with cognitive control.Moreover, we observed that the community structures detected for the error response were more consistent across subjects compared to the community structures for correct response.
Future work will consider extension of this multilayer model to higher dimensions, e.g.multi-aspect multilayer brain networks such as temporal multi-frequency connectivity networks.Compared to current work where subjects' community structure is found separately and then combined through multiview community detection, future work can use multi-aspect multilayer networks constructed from subjects' multilayer networks.This approach will allow simultaneous detection of communities of subjects similar to [85].Future work will also consider different null models in the definition of modularity such as the constant Potts model, which is shown to be resolution limit free [86].Finally, in this work we aimed to find the optimal resolution and inter-layer scale parameter; future work can focus on a multi-scale approach where the aim is to combine community structures from different resolutions and inter-layer scales [87].

Figure 1 .
Figure 1.Flowchart of the proposed approach for community detection of multi-frequency EEG networks.Bottom two panels illustrate multilayer network construction (left) and community detection for each subject (right).

Figure 2 .
Figure2.Supra-adjacency and two types of modularity matrices for multilayer EEG networks for a single subject.Top row corresponds to the error response, while the bottom row corresponds to the correct response.a) Supra-adjacency matrices, b) Modularity matrices obtained when the multilayer network is considered as a single-layer network, c) Modularity matrices obtained using the modularity function defined in(11).Rows and columns of the matrices are color coded to indicate the layers corresponding to the nodes.

Figure 3 .
Figure 3. Selection of the resolution (γ) and inter-layer scale (ω) parameters: a) and b) show the average of Q obs − Qsurr across 20 subjects for error and correct responses, respectively.c) shows the histogram of optimal ω values for error (top) and correct (bottom) responses across subjects.

Figure 4 .
Figure 4. Consistency of the community structure for error and correct responses as measured by JS distance.

Figure 5 .
Figure 5. Multilayer community structures for error (a)) and correct (b)) responses.Each electrode is shown with a circle that includes 4 colored circles, representing communities of electrodes across 4 frequency bands.c) shows the correspondence of inner circles to frequency bands.