## Abstract

Network science has been extensively developed to characterize the structural properties of complex systems, including brain networks inferred from neuroimaging data. As a result of the inference process, networks estimated from experimentally obtained biological data represent one instance of a larger number of realizations with similar intrinsic topology. A modelling approach is therefore needed to support statistical inference on the bottom-up local connectivity mechanisms influencing the formation of the estimated brain networks. Here, we adopted a statistical model based on exponential random graph models (ERGMs) to reproduce brain networks, or connectomes, estimated by spectral coherence between high-density electroencephalographic (EEG) signals. ERGMs are made up by different local graph metrics, whereas the parameters weight the respective contribution in explaining the observed network. We validated this approach in a dataset of *N* = 108 healthy subjects during eyes-open (EO) and eyes-closed (EC) resting-state conditions. Results showed that the tendency to form triangles and stars, reflecting clustering and node centrality, better explained the global properties of the EEG connectomes than other combinations of graph metrics. In particular, the synthetic networks generated by this model configuration replicated the characteristic differences found in real brain networks, with EO eliciting significantly higher segregation in the *alpha* frequency band (8–13 Hz) than EC. Furthermore, the fitted ERGM parameter values provided complementary information showing that clustering connections are significantly more represented from EC to EO in the *alpha* range, but also in the *beta* band (14–29 Hz), which is known to play a crucial role in cortical processing of visual input and externally oriented attention. Taken together, these findings support the current view of the functional segregation and integration of the brain in terms of modules and hubs, and provide a statistical approach to extract new information on the (re)organizational mechanisms in healthy and diseased brains.

## 1. Introduction

The study of the human brain at rest provides precious information that is predictive of intrinsic functioning, cognition, as well as pathology [1]. In the last decade, graph theoretic approaches have described the topological structure of resting-state connectomes derived from different neuroimaging techniques, such as functional magnetic resonance imaging (fMRI) or magneto- (MEG) and electroencephalography (EEG).

These estimated connectomes, or brain networks, tend to exhibit similar organizational properties, including small-worldness, cost-efficiency, modularity and node centrality [2], as well as characteristic dependence from the anatomical backbone connectivity [3–5] and genetic factors [6]. Furthermore, they potentially show clinical relevance, as demonstrated by the recent development of network-based diagnostics of consciousness [7,8], Alzheimer's disease [9], stroke recovery [10] and schizophrenia [11]. In this sense, quantifying the topological properties of intrinsic functional connectomes by means of graph theory has enriched our understanding of the structure of functional brain connectivity maps [2,12–14]. Nevertheless, these results refer to a descriptive analysis of the observed brain network, which is only one instance of several alternatives with similar structural features. This is especially true for functional networks inferred from empirically obtained data, where the edges (or links) are noisy estimates of the true connectivity and thresholding is often adopted to filter the relevant interactions between the system units [15–17].

Statistical models are, therefore, needed to reflect the uncertainty associated with a given observation, to permit inference about the relative occurrence of specific local structures and to relate local-level processes to global-level properties [18]. A first approach consists in generating synthetic random networks that preserve some observed properties, such as the degree distribution or the random walk distribution, and then contrasting the values of the graph indices obtained in these synthetic networks with those extracted from the estimated connectomes [13]. While these methods often provide appropriate null models, and can improve the identification of relevant network properties [19–21], they do not inform on the organizational mechanisms modelling the whole network formation [22,23]. Alternative approaches consider probabilistic growth models such as those based on spatial distances between nodes [24]. Interesting results have been achieved in identifying some basic connectivity rules reproducing both structural and functional brain networks [25,26]. However, these methods suffer from the rough approximation (e.g. Euclidean) of the actual spatial distance between nodes, and, moreover, they do not indicate whether the identified local mechanisms are either necessary or sufficient as descriptors of the global network structure.

To support inference on the processes influencing the formation of network structure, statistical models have been conceived to consider the set of all possible alternative networks weighted on their similarity to the observed one [18]. Among others, exponential random graph models (ERGMs) represent a flexible category that allows the simultaneous assessment of the role of specific graph features in the formation of the entire network. These models were first proposed as an extension of the triad model defined in [27] to characterize Markov graphs [28,29] and have been widely developed to understand how simple interaction rules, such as transitivity, could give rise to the complex network of social contacts [30–38].

Recently, the use of ERGMs has been proved to successfully model imaging connectomes derived respectively from spontaneous fMRI activity [39] and diffusion tensor imaging (DTI) [40]. Despite its potential, the use of ERGMs in network neuroscience is still in its infancy and more evidence is needed to better elucidate its applicability to connectomes inferred from other types of neuroimaging data and across different experimental conditions. In addition, many methodological issues remain unanswered, such as the relationships between the graph metrics included in the ERGM and the graph indices used to describe the topology of the observed connectomes.

To address the above issues, we proposed and evaluated several ERGM configurations based on the combination of different local connectivity structures (i.e. graph metrics). Specifically, we modelled brain networks estimated from high-density EEG signals in a group of healthy individuals during eyes-open (EO) and eyes-closed (EC) resting states. Our goal was to identify the best ERGM configuration reproducing EEG-derived connectomes in terms of functional integration and segregation, and to evaluate the ability of the estimated ERGM parameters in providing new information discriminating between EO and EC conditions.

## 2. Material and methods

### 2.1. Electroencephalographic data and brain network construction

We used high-density EEG signals freely available from the online PhysioNet BCI database [41,42]. EEG data consisted of 1 min resting state with EO and 1 min resting state with EC recorded from 56 electrodes in 108 healthy subjects. EEG signals were recorded with an original sampling rate of 160 Hz. All the EEG signals were referenced to the mean signal gathered from electrodes on the ear lobes. We subsequently downsampled the EEG signals to 100 Hz after applying a proper anti-aliasing low-pass filter. The electrode positions on the scalp followed the standard 10–10 montage.

We used the spectral coherence [43] to measure functional connectivity (FC) between EEG signals of sensors *i* and *j* at a specific frequency band *f* as follows:
2.1where *S*_{ij} is the cross-spectrum between *i* and *j*, and *S*_{ii} and *S*_{jj} are the autospectra of *i* and *j*, respectively. Specifically, we computed cross- and auto-spectra by means of Welch's averaged modified periodogram with a sliding Hanning window of 1 s and 0.5 s of overlap. The number of fast Fourier transform points was set to 100 for a frequency resolution of 1 Hz. As a result, we obtained for each subject a connectivity matrix *W*(*f*) of size 56 × 56 where the entry *w*_{ij}(*f*) contains the value of the spectral coherence between the EEG signals of sensors *i* and *j* at the frequency *f*.

We then averaged the connectivity matrices within the characteristic frequency bands *theta* (4–7 Hz), *alpha* (8–13 Hz), *beta* (14–29 Hz) and *gamma* (30–40 Hz). These matrices constituted our raw brain networks whose nodes corresponded to the EEG sensors (*n* = 56) and links corresponded to the *w*_{ij} values. Finally, we thresholded the values in the connectivity matrices to retain the strongest links in each brain network. Specifically, we adopted an objective criterion, i.e. the efficiency cost optimization (ECO), to filter and binarize a number of links such that the final average node degree *k* = 3 [44]. We also considered *k* = 1, 2, 4, 5 to evaluate the main brain network properties around the representative threshold *k* = 3. The resulting sparse brain networks, or graphs, were represented by adjacency matrices *A*, where each entry indicates the presence *a*_{ij} = 1 or the absence *a*_{ij} = 0 of a link between nodes *i* and *j*.

### 2.2. Graph indices

We evaluated the global structure of brain networks by measuring graph indices at large-scale topological scales. We focused on well-known properties of brain networks such as optimal balance between integration and segregation of information [2,45,46]. Integration is the tendency of the network to favour distributed connectivity among remote brain areas; conversely, segregation is the tendency of the network to maintain connectivity within specialized groups of brain areas [47].

In graph theory, integration has been typically quantified by the global-efficiency *E*_{g} and by the characteristic path length *L*,
2.2where *d*_{ij} is the distance, or the length of the shortest path, between nodes *i* and *j* [48,49].

Segregation is typically measured by means of the local-efficiency *E*_{l} and by the clustering coefficient *C*:
2.3where *G*_{i} is the subgraph formed by the nodes connected to *i*; *t*_{i} is the number of triangles around node *i*; and *k*_{i} is the degree of node *i* [48,49].

In addition, we evaluated the strength of division of a network into modules by measuring the modularity *Q*:
2.4where is the number of edges, *m*_{i} is the module containing node *i* and *δ*_{mi,mj} = 1 if *m*_{i} = *m*_{j} and 0 otherwise. We used the Walktrap algorithm to generate a sequence of community partitions [50] and we selected the one that maximized *Q* according to the standard algorithm proposed in [51]. Modularity can be seen as a compact measure of the integration and segregation of a network, as it measures the propensity to form dense connections between nodes within modules (i.e. segregation) but sparse connections between nodes in different modules (i.e. inverse of integration).

### 2.3. Exponential random graph model

Let *G* be a graph in a set of possible network realizations, *g* = [*g*_{1}, *g*_{2}, … , *g*_{r}] be a vector of graph statistics, or metrics, and *g** = [*g**_{1}, *g**_{2}, … , *g**_{r}] be the values of these metrics measured over *G*. Then, we can statistically model *G* by defining a probability distribution *P*(*G*) over such that the following conditions are satisfied:
2.5and
2.6 where 〈*g*_{i}〉 is the expected value of the *i*th graph metric over .

By maximizing the Gibbs entropy of *P*(*G*) constrained to the above conditions, the probability distribution reads as:
2.7 where is the graph Hamiltonian, *θ*_{i} is the *i*th model parameter to be estimated and is the so-called partition function [52]. The estimated value of a parameter *θ*_{i} indicates the change in the (log-odds) likelihood of an edge for a unit change in graph metric *g*_{i}. If the estimated value of *θ*_{i} is large and positive, the associated graph metric *g*_{i} plays an important role in explaining the topology of *G* more than would be expected by chance. Note that here chance corresponds to randomly choosing a network from the space . If instead the estimated value of *θ*_{i} is negative and large, then *g*_{i} still plays an important role in explaining the topology of *G* but it is less prevalent than expected by chance [53].

In general, the fact that the space can be very large even for relatively small *n*, as well as the inclusion of graph metrics that are not simple linear combinations of *G*_{ij}, in practice make it impossible to derive analytically the model parameters vector ** θ** = [

*θ*

_{1},

*θ*

_{2}, … ,

*θ*

_{r}] [27,31].

Numerical methods, such as Markov chain Monte Carlo (MCMC) approximations of the maximum-likelihood estimators (MLEs) of the model parameters vector ** θ**, are typically adopted to circumvent this issue [54].

#### 2.3.1. Model construction and implementation

We considered graph metrics reflecting the basic properties of complex systems such as hub propensity and transitivity in the network [46,55,56]. Specifically, we focused on *k*-stars to model highly connected nodes (hubs) and *k*-triangles to model transitivity, where *k* refers to the order of the structures as illustrated in figure 1.

In general, this leads to a large number of model parameters to be estimated, i.e. *n* − 1 for *k*-stars and *n* − 2 for *k*-triangles. To avoid consequent degeneracy issues in the ERGM estimation, we adopted a compact specification for these metrics that combines them in an alternating geometric sequence [31,57].

Because *k*-stars are related to the node degree distribution *D* [33], we used the *geometrically weighted degree distribution* GW_{K} as a graph metric to characterize hub propensity:
2.8where *τ* > 0 is a *ratio parameter* to penalize nodes with extremely high node degrees.

Similarly, because *k*-triangles are related to the *shared pattern distribution* *S*, we used the *geometrically weighted edgewise shared partner distribution* to characterize transitivity:
2.9where the element *S*_{i} is the number of dyads that are directly connected and that have exactly *i* neighbours in common.

In addition, complementary metrics have been defined based on the *shared partner distribution*:

GW

_{N}:*geometrically weighted non-edgewise shared partner distribution*given by equation (2.9), with*S*_{i}considering exclusively dyads that are not connected.GW

_{D}:*geometrically weighted dyadwise shared partner distribution*given by equation (2.9), with*S*_{i}considering any dyad, connected or not.

The above specifications yield particular ERGMs that belong to the so-called curved exponential family [33] and that have been extensively used in social science [32,58,59].

We constructed different ERGM configurations by including these graph metrics as illustrated in table 1. For the sake of simplicity, we only considered combinations of two graph metrics at most, except in one case where we also included the number of edges as a further metric [39,40].

We tested the different configurations by fitting the ERGM to brain networks in each single subject (*N* = 108), frequency band (*theta*, *alpha*, *beta*, *gamma*) and condition (EO, EC). To fit ERGMs, we used an MCMC algorithm (Gibbs sampler) that samples networks from an exponential graph distribution. Specifically, we set the initial values of the model parameters ** θ^{0}** by means of a maximum pseudo-likelihood estimation (MPLE) [54,60]. Then, we adopted Fisher's scoring method to update the model parameters

**until they converged to the approximated MLEs [31]. As we used curved ERGMs, the ratio parameters**

*θ**τ*were not fixed but were estimated.

Eventually, for each fitted ERGM configuration we generated 100 synthetic networks in order to obtain appropriate confidence intervals.

#### 2.3.2. Goodness of fit

First, we used the Akaike information criterion (AIC) to evaluate the relative quality of the ERGMs' fit by taking into account the maximum value of the likelihood function and the number of model parameters [61].

We also adopted a different approach to assess the absolute quality of the fit by comparing the synthetic networks generated by the estimated ERGMs and the observed brain networks. Specifically, we defined the following score based on the integration and segregation properties of networks:
2.10where *η*_{Eg},*η*_{El} are the relative errors between the mean values of the global/local efficiency of the simulated networks and the value of the observed brain network. By selecting the maximum absolute error, we were considering the worst case, similar to what was proposed in [26]. Based on the above criteria, we selected the best model, which minimizes the AIC and *δ* mean values. To validate the model adequacy (equation (2.6)), we computed the *Z*-scores between the graph metrics' values of brain networks and synthetic networks.

Furthermore, we cross-validated the best model configuration by evaluating the synthetic networks' fit to graph indices that were explicitly not included in either the ERGM or the model selection criteria. We computed Pearson's correlation coefficient between the values of the characteristic path length (*L*), clustering coefficient (*C*) and modularity (*Q*) extracted from the observed brain networks and the mean values obtained from the corresponding simulated networks. In addition, we used the Mirkin index (MI) [62] to evaluate the similarity between the community partitions of the observed networks and the consensus partitions of the corresponding synthetic networks.

### 2.4. Statistical group analysis

We assessed the statistical differences between the values of the graph indices extracted from the brain networks in the EO and EC resting-state conditions. We also computed between-condition differences using the synthetic networks fitted by the best ERGM. In this case, we considered the mean values of the graph indices in order to have one value corresponding to one brain network. Eventually, we computed the statistical differences between the values of the best ERGM parameters in the EO and EC conditions in order to assess their potential to provide complementary information to that provided by standard graph analysis. For each comparison, we used a non-parametric permutation *t*-test and we fixed a statistical threshold of *α* = 0.001 and 100 000 permutations.

## 3. Results

### 3.1. Characteristic functional segregation of electroencephalographic resting-state networks

The group analysis revealed a significant increase in the local-efficiency in EO, compared with that in EC, for the *alpha* band (*T* = 3.529, *p* = 0.0007, figure 2). We also reported a significant increment (*T* = 3.557, *p* = 0.0007) for the modularity in the *alpha* band, while no other statistically significant differences were observed in the other frequency bands, graph indices or metrics (electronic supplementary material, table S3).

These differences were obtained for brain networks thresholded with an average node degree *k* = 3 according to the ECO criterion [44]. We reported a similar increase in functional segregation (local-efficiency) in the *alpha* band for *k* = 5 (electronic supplementary material, figure S1). More details on the analysis for *k* = 5 can be found in the electronic supplementary material (Supp_text.pdf).

In terms of existing relationships between graph indices and ERGM metrics, we could not establish univocal associations between *E*_{g} and *E*_{l} values and the metrics' values used in the ERGMs (electronic supplementary material, table S1). This was especially true for the global-efficiency, which exhibited significantly high correlations with all the other graph metrics (Spearman's |*R*| > 0.43, *p* < 10^{−39}).

### 3.2. Triangles and stars as fundamental constituents of functional brain networks

All the proposed ERGM configurations exhibited a relatively good fitting in terms of AIC, except for *M*_{11} (electronic supplementary material, figure S2). Notably, the latter was the only configuration where the number of edges was considered as a model parameter and not as a constraint. *M*_{1} gave the lowest *δ*(*E*_{g}, *E*_{l}) scores compared with the other configurations in both the EO and EC conditions (figure 3). Notably, the configurations giving lower *δ*(*E*_{g}, *E*_{l}) scores included, directly or indirectly, the metric GW_{E}, with the exception of *M*_{11}.

We selected *M*_{1} as a potentially good candidate to model EEG-derived brain networks. According to this model configuration the mass probability density reads . The group-median values of the estimated parameters (*θ*_{1} and *θ*_{2}) were all positive and larger than 1 in each band and condition (table 2). This means that the likelihood of an edge existing in a simulated network is larger if that edge is part of a triangle (GW_{E}) or of a star (GW_{K}), and that these connectivity structures are statistically relevant for the brain network formation.

Overall, the GW_{E} and GW_{K} values of the synthetic networks generated by *M*_{1} were not significantly different from those of the observed brain networks (figure 4). This was true in every subject for GW_{E} (*Z* < 2.58, *p* > 0.01) and in at least 94% of the subjects for GW_{E} (*Z* < 2.58, *p* > 0.01). Furthermore, the values of the characteristic path length (*L*), clustering coefficient (*C*) and modularity (*Q*) extracted from synthetic networks were significantly correlated (Pearson's *R* > 0.44, *p* < 10^{−6}) with those of the brain networks in each frequency band (figure 5; electronic supplementary material, table S2). In addition, synthetic networks exhibited a similar community partition to individual brain networks, as revealed by the low MI values (MI < 0.21) (electronic supplementary material, figure S3). These results confirmed that *M*_{1} adequately models the obtained EEG brain networks.

### 3.3. Simulating network differences between absence and presence of visual input

Figure 6 illustrates the brain networks for a representative subject in the *alpha* band along with the corresponding synthetic networks generated by *M*_{1}. In both the EO and EC conditions, simulated networks and brain networks share similar topological structures characterized by diffused regularity and more concentrated connectivity in parietal and occipital regions.

The group analysis over the synthetic networks revealed the ability of *M*_{1} to capture not only the individual properties of brain networks but also the main observed difference between the EC and EO resting states, reflecting, respectively, the absence and presence of visual input. Similarly to observed brain networks, we obtained, for simulated networks, a marginally significant increase in the local-efficiency from EC to EO, in the *alpha* band (*T* = 3.168, *p* = 0.002). No other significant differences were reported in any other band or graph index/metric (electronic supplementary material, table S3).

Finally, by looking at the values of the estimated parameters, we observed that *θ*_{1} values were significantly larger in EO than in EC for both the *alpha* (*T* = 3.746, *p* = 0.0002) and *beta* (*Z* = 1.514, *p* = 0.0009) frequency bands, while no significant differences were found for *θ*_{2} values (table 2).

## 4. Discussion

In recent years, the use of statistical methods to infer the structure of complex systems has gained increasing interest [39,64–66]. Beyond the descriptive characterization of networks, statistical network models aim to statistically assess the local connectivity processes involved in the global structure formation [18]. This is a crucial advance with respect to standard descriptive approaches because imaging connectomes, as with other biological networks, is often inferred from experimentally obtained data and therefore the estimated edges can suffer from statistical noise and uncertainty [67].

In our study, we used ERGMs to identify the local connectivity structures that statistically form the intrinsic synchronization of large-scale electrophysiological activities. This model formulation has the advantage of statistically inferring the probability of edge formation accounting for highly dependent configurations, such as transitivity structures, something that is lacking in, for example, the Bernoulli model. Furthermore, it is possible to include, in theory, graph metrics measuring global and local properties and discriminating node and edges attributes, such as homophily effects. In addition, it generalizes well-known network models such as the stochastic block model, where a block structure is imposed by including the count of edges between groups of nodes as a model metric [68].

Here, the results showed that the tendency to form triangles (GW_{E}) and stars (GW_{K}) was sufficient to statistically reproduce the main properties of the EEG brain networks, such as functional integration and segregation, measured by means of global-efficiency *E*_{g} and local-efficiency *E*_{l} (electronic supplementary material, table S3). Our findings partially deviate from previous studies, which have used ERGMs to model fMRI and DTI brain networks, where GW_{E} and the geometrically weighted non-edgewise shared partner GW_{N} were selected under the assumption that these could be related, respectively, to local- and global-efficiency [39,40]. However, here we showed that a univocal relationship between the ERGM graph indices and the metrics used to describe the EEG connectomes could not be statistically established (electronic supplementary material, table S1). While the propensity to form triangles (GW_{E}) can lead to cohesive clustering in the network (*E*_{l}), the propensity to form redundant paths of length 2 (i.e. GW_{N}) is not clearly related to the formation of short paths between nodes (*E*_{g}) [57]. Thus, while in general a good fit can be achieved by including GW_{N} in the ERGM, the subsequent interpretation in terms of brain functional integration appears less straightforward. Here, we showed that GW_{E} together with the tendency to form stars (GW_{K}) gave the best fit in terms of local- and global efficiency. Triangles and stars, giving rise to clustering and hubs, are fundamental building blocks of complex systems reflecting important mechanisms such as transitivity [48] and preferential attachment [69]. Notably, the existence of highly connected nodes is compatible with the presence of short paths (e.g. in a star graph the characteristic path length *L* = 2). This supports the recent view of brain functional integration where segregated modules exchange information through central hubs and not necessarily through the shortest paths [70,71].

In the cross-validation phase, the selected model configuration captured other important brain network properties as measured by the clustering coefficient *C*, the characteristic path length *L* and the modularity *Q* (figure 5). In terms of the differences between conditions, the simulated networks gave a marginally significant increase (*p* = 0.002) in *E*_{l} in the *alpha* band during EO as compared with EC, while, differently from observed brain networks, no significant differences were reported for the modularity *Q* (electronic supplementary material, table S3). The latter could be, in part, ascribed to the absence of specific metrics in the ERGM accounting for modularity. In this respect, stochastic block models, which explicitly force modular structures, could represent an interesting alternative to explore in the future [72,73]. Here, the increased *alpha* local-efficiency suggests a modulation of augmented specialized information processing, from EC to EO, that is consistent with typical global power reduction and increased regional activity [74]. Possible neural mechanisms explaining this effect have been associated with the automatic gathering of non-specific information resulting from more interactions within the visual system [75] and with shifts from interoceptive towards exteroceptive states [76–78].

As a crucial result, we provided complementary information by inspecting the fitted ERGM parameters. The positive *θ*_{1} > 1 and *θ*_{2} > 1 values indicated that both GW_{E} and GW_{K} are fundamental connectivity features that emerge in brain networks more than expected by chance (table 2). However, only *θ*_{1} values showed a significant difference (EO > EC) in the *alpha* band, as well as in the *beta* band (table 2), suggesting that the tendency to form triangles, rather than the tendency to form stars, is a discriminating feature of EO and EC modes. More concentrated EEG activity among parieto-occipital areas has been largely documented in the *alpha* and also in the *beta* bands, the latter reflecting either cortical processing of visual input or externally oriented attention [74,79]. Notably, the role of the *beta* band could not be found when analysing either brain networks or synthetic networks (figure 2; electronic supplementary material, table S3) and we speculate that this result specifically stems from the inherent ability of ERGMs to account for potential interaction between different graph metrics [57].

### 4.1. Methodological considerations

We estimated EEG connectomes by means of spectral coherence. While this measure is known to suffer from possible volume conduction effects [80], it has also been demonstrated that, probably due to this effect, it has the advantage of generating connectivity matrices that are highly consistent within and between subjects [81]. In addition, spectral coherence is still one of the most used measures to infer FC in the electrophysiological literature on resting states because of its simplicity and relatively intuitive interpretation. Thus, constructing EEG connectomes by means of spectral coherence allowed us to better contextualize the results obtained with ERGM from a neurophysiological perspective. Future studies will have to assess if and how different connectivity estimators affect the choice of the model parameters.

We used a density-based thresholding procedure to filter information in the EEG raw networks by retaining and binarizing the strongest edges. Despite the consequent information loss, thresholding is often adopted to mitigate the uncertainty of the weakest edges, reduce the false positives and facilitate the interpretation of the inferred network topology [13,17].

Selecting a binarizing threshold does have an impact on the topological structure of brain networks [82]. Based on the optimization of fundamental properties of complex systems, i.e. efficiency and economy, the adopted thresholding criterion (ECO) leads to sparse networks, with an average node degree *k* = 3, causing possible nodes to be disconnected. However, it has demonstrated empirically that the size of the resulting largest component typically contains more than 60% of the brain nodes, thus ensuring a sparse but meaningful network structure [44]. In a separate analysis, we verified that the validity of the model and the characteristic between-condition differences observed in the *alpha* band were also globally preserved for *k* = 5 (electronic supplementary material, Supp_text.pdf).

## Authors' contributions

C.O. and F.D.V.F. both designed the study; C.O. performed the analysis and wrote the paper; and F.D.V.F. wrote the paper.

## Competing interests

We declare we have no competing interests.

## Funding

This work has been partially supported by the Agence Nationale de la Recherche (French programme) ANR-10-IAIHU-06 and ANR-15-NEUC-0006-02.

## Acknowledgements

We are grateful to M. Chavez and J. Guillon for their useful comments and suggestions.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3695518.

- Received November 23, 2016.
- Accepted February 10, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.