# Effect of channel density, inverse solutions and connectivity measures on EEG resting-state networks reconstruction: A simulation study

Sahar Allouch, Aya Kabbara, Joan Duprez, Mohamad Khalil, Julien Modolo, Mahmoud Hassan

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sahar Allouch, Aya Kabbara, Joan Duprez, Mohamad Khalil, Julien Modolo, et al., Effect of channel density, inverse solutions and connectivity measures on EEG resting-state networks reconstruction: A simulation study. NeuroImage, 2023, 271, pp.120006. 10.1016/j.neuroimage.2023.120006. hal-04060943

HAL Id: hal-04060943

https://hal.science/hal-04060943

Submitted on 29 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_1.png" /></figure>

# Effect of channel density, inverse solutions and connectivity measures on EEG resting-state networks reconstruction: A simulation study

<figure><img src="image_2.png" /></figure>

Sahar Allouch<sup>a,b,*</sup>, Aya Kabbara<sup>c,d</sup>, Joan Duprez<sup>a</sup>, Mohamad Khalil<sup>b,e</sup>, Julien Modolo<sup>a,1</sup>, Mahmoud Hassan<sup>c,f,1</sup>

## ARTICLE INFO

**Keywords:**

EEG resting-state networks

Channel density

Inverse solution

Functional connectivity

Neural mass model

Analytical variability

## ABSTRACT

Along with the study of brain activity evoked by external stimuli, the past two decades witnessed an increased interest in characterizing the spontaneous brain activity occurring during resting conditions. The identification of connectivity patterns in this so-called "resting-state" has been the subject of a great number of electrophysiology-based studies, using the Electro/Magneto-Encephalography (EEG/MEG) source connectivity method. However, no consensus has been reached yet regarding a unified (if possible) analysis pipeline, and several involved parameters and methods require cautious tuning. This is particularly challenging when different analytical choices induce significant discrepancies in results and drawn conclusions, thereby hindering the reproducibility of neuroimaging research. Hence, our objective in this study was to shed light on the effect of analytical variability on outcome consistency by evaluating the implications of parameters involved in the EEG source connectivity analysis on the accuracy of resting-state networks (RSNs) reconstruction. We simulated, using neural mass models, EEG data corresponding to two RSNs, namely the default mode network (DMN) and dorsal attentional network (DAN). We investigated the impact of five channel densities (19, 32, 64, 128, 256), three inverse solutions (weighted minimum norm estimate (wMNE), exact low-resolution brain electromagnetic tomography (eLORETA), and linearly constrained minimum variance (LCMV) beamforming) and four functional connectivity measures (phase-locking value (PLV), phase-lag index (PLI), and amplitude envelope correlation (AEC) with and without source leakage correction), on the correspondence between reconstructed and reference networks. We showed that, with different analytical choices related to the number of electrodes, source reconstruction algorithm, and functional connectivity measure, high variability is present in the results. More specifically, our results show that a higher number of EEG channels significantly increased the accuracy of the reconstructed networks. Additionally, our results showed significant variability in the performance of the tested inverse solutions and connectivity measures. Such methodological variability and absence of analysis standardization represent a critical issue for neuroimaging studies that should be prioritized. We believe that this work could be useful for the field of electrophysiology connectomics, by increasing awareness regarding the challenge of variability in methodological approaches and its implications on reported results.

## 1. Introduction

Over the past two decades, there has been a considerable growth in the number of studies investigating human brain activity at rest (Raichle et al., 2001; Raichle and Snyder 2007; van den Heuvel et al., 2008; Damoiseaux et al., 2006). Characterizing synchronous activity

across spatially distributed regions, using either functional magnetic resonance imaging (fMRI) or magneto/electro-encephalography (M/EEG) has revealed consistent patterns of brain connectivity in the absence of a goal-directed task, with a preserved consistency across subjects as well as across different neuroimaging modalities (Damoiseaux et al., 2006; Brookes et al., 2011; Francesco de Pasquale et al., 2012; Kabbara et al.,

<footnote><sup>a</sup> Univ Rennes, INSERM, LTSI - UMR 1099, Rennes F-35000, France</footnote>

<footnote><sup>b</sup> Azm Center for Research in Biotechnology and Its Applications, EDST, Tripoli, Lebanon</footnote>

<footnote><sup>c</sup> MINDIG, Rennes F-35000, France</footnote>

<footnote><sup>d</sup> LASeR - Lebanese Association for Scientific Research, Tripoli, Lebanon</footnote>

<footnote><sup>e</sup> CRSI research center, Faculty of Engineering, Lebanese University, Beirut, Lebanon</footnote>

<footnote><sup>f</sup> School of Science and Engineering, Reykjavik University, Reykjavik, Iceland</footnote>

<footnote>* Corresponding author at: Univ Rennes, INSERM, LTSI - UMR 1099, Rennes F-35000, France.
E-mail address: saharallouch@gmail.com (S. Allouch).</footnote>

<footnote n="1">¹ These authors contributed equally to this work.</footnote>
<figure><img src="image_3.png" /></figure>

# Effect of channel density, inverse solutions and connectivity measures on EEG resting-state networks reconstruction: A simulation study

<figure><img src="image_4.png" /></figure>

Sahar Allouch<sup>a,b,*</sup>, Aya Kabbara<sup>c,d</sup>, Joan Duprez<sup>a</sup>, Mohamad Khalil<sup>b,e</sup>, Julien Modolo<sup>a,1</sup>, Mahmoud Hassan<sup>c,f,1</sup>

## ARTICLE INFO

**Keywords:**

EEG resting-state networks

Channel density

Inverse solution

Functional connectivity

Neural mass model

Analytical variability

## ABSTRACT

Along with the study of brain activity evoked by external stimuli, the past two decades witnessed an increased interest in characterizing the spontaneous brain activity occurring during resting conditions. The identification of connectivity patterns in this so-called "resting-state" has been the subject of a great number of electrophysiology-based studies, using the Electro/Magneto-Encephalography (EEG/MEG) source connectivity method. However, no consensus has been reached yet regarding a unified (if possible) analysis pipeline, and several involved parameters and methods require cautious tuning. This is particularly challenging when different analytical choices induce significant discrepancies in results and drawn conclusions, thereby hindering the reproducibility of neuroimaging research. Hence, our objective in this study was to shed light on the effect of analytical variability on outcome consistency by evaluating the implications of parameters involved in the EEG source connectivity analysis on the accuracy of resting-state networks (RSNs) reconstruction. We simulated, using neural mass models, EEG data corresponding to two RSNs, namely the default mode network (DMN) and dorsal attentional network (DAN). We investigated the impact of five channel densities (19, 32, 64, 128, 256), three inverse solutions (weighted minimum norm estimate (wMNE), exact low-resolution brain electromagnetic tomography (eLORETA), and linearly constrained minimum variance (LCMV) beamforming) and four functional connectivity measures (phase-locking value (PLV), phase-lag index (PLI), and amplitude envelope correlation (AEC) with and without source leakage correction), on the correspondence between reconstructed and reference networks. We showed that, with different analytical choices related to the number of electrodes, source reconstruction algorithm, and functional connectivity measure, high variability is present in the results. More specifically, our results show that a higher number of EEG channels significantly increased the accuracy of the reconstructed networks. Additionally, our results showed significant variability in the performance of the tested inverse solutions and connectivity measures. Such methodological variability and absence of analysis standardization represent a critical issue for neuroimaging studies that should be prioritized. We believe that this work could be useful for the field of electrophysiology connectomics, by increasing awareness regarding the challenge of variability in methodological approaches and its implications on reported results.

## 1. Introduction

Over the past two decades, there has been a considerable growth in the number of studies investigating human brain activity at rest (Raichle et al., 2001; Raichle and Snyder 2007; van den Heuvel et al., 2008; Damoiseaux et al., 2006). Characterizing synchronous activity

across spatially distributed regions, using either functional magnetic resonance imaging (fMRI) or magneto/electro-encephalography (M/EEG) has revealed consistent patterns of brain connectivity in the absence of a goal-directed task, with a preserved consistency across subjects as well as across different neuroimaging modalities (Damoiseaux et al., 2006; Brookes et al., 2011; Francesco de Pasquale et al., 2012; Kabbara et al.,

<footnote><sup>a</sup> Univ Rennes, INSERM, LTSI - UMR 1099, Rennes F-35000, France</footnote>

<footnote><sup>b</sup> Azm Center for Research in Biotechnology and Its Applications, EDST, Tripoli, Lebanon</footnote>

<footnote><sup>c</sup> MINDIG, Rennes F-35000, France</footnote>

<footnote><sup>d</sup> LASeR - Lebanese Association for Scientific Research, Tripoli, Lebanon</footnote>

<footnote><sup>e</sup> CRSI research center, Faculty of Engineering, Lebanese University, Beirut, Lebanon</footnote>

<footnote><sup>f</sup> School of Science and Engineering, Reykjavik University, Reykjavik, Iceland</footnote>

<footnote>* Corresponding author at: Univ Rennes, INSERM, LTSI - UMR 1099, Rennes F-35000, France.
E-mail address: saharallouch@gmail.com (S. Allouch).</footnote>

<footnote n="1">¹ These authors contributed equally to this work.</footnote>
2021). Interestingly, accurate identification of the so-called resting-state networks (RSNs) has been found to be primordial for both cognitive (Keller et al., 2015; Jockwitz et al., 2017; Shen et al., 2018; Alavash et al., 2015) and clinical studies (Sheline and Raichle 2013; Kabbara et al., 2018; Hassan et al., 2017; Gratton et al., 2018; Jimenez et al., 2019).

In this context, the EEG source connectivity method enables inferring functional brain networks at the cortical level with high temporal resolution (Hassan and Wendling, 2018). However, due to the complexity of the analysis workflow, many methodological questions related to the EEG source connectivity pipeline remain unanswered. In fact, each step of the analysis involves different choices that might significantly affect the resulting functional network. Hence, the effect of several factors needs to be investigated and quantified for more consistent and reliable use of this method.

One major factor in terms of impact on the EEG source connectivity analysis is the number of EEG electrodes (electrodes density). In fact, several studies investigated the effect of electrodes density on EEG source localization in simulations in the context of epilepsy (Song et al., 2015; Sohrabpour et al., 2015; Goran Lantz et al., 2003). More specifically, it has been shown that EEG electrodes density has a direct influence on the localization error: as intuitively expected, a higher number of electrodes is associated with a significant decrease in localization error (Goran Lantz et al., 2003; Song et al., 2015; Sohrabpour et al., 2015). In both studies by Goran Lantz et al. (2003) and Sohrabpour et al. (2015), a dramatic decrease in localization error occurred when increasing the number of electrodes from 32 to 64.

Another critical influencing factor in the EEG source connectivity pipeline is the algorithm chosen to solve the inverse problem, which is ill-posed due to its non-uniqueness and the instability of its solution (see Grech et al. (2008) for a review). Several studies quantifying the performance of different inverse methods, in simulated and experimental EEG/MEG data, concluded that the choice of the inverse method significantly influences source estimation results (Anzolin et al., 2019; Mahjoory et al., 2017; Hedrich et al., 2017; Bradley et al., 2016; Grova et al., 2006; Halder et al., 2019; Tait et al., 2021; Allouch et al., 2022). However, no consistent conclusions have been made regarding one method that would stand apart from the others in terms of performance, which can also be related to the analyzed conditions.

The choice of the functional connectivity metric is also a critical step. A wide range of measures are used in the field, and each differs in the aspect of the data being investigated (amplitude- vs phase-based measures / directional vs non-directional connectivity, prone/robust to source leakage), (see Friston (2011), Pereda et al. (2005), Cao et al. (2022) for a review), resulting in significant variability of performance and interpretations (Colclough et al., 2016; Wang et al., 2014; Wendling et al., 2009; Hassan et al., 2017; Allouch et al., 2022).

In this context, simulation studies are of utmost interest, since they provide an accessible ground truth that enables an objective evaluation of the methods/techniques under investigation. Here, we used a model for cortical dynamics, namely neural mass models (NMMs), to simulate resting state networks (RSNs), and more specifically the default mode network (DMN) and dorsal attentional network (DAN). Then, we quantified the effect of three key factors involved in the EEG source connectivity analysis:

(1) EEG channel density (from 19,32,64,128 to 256 channels).

(2) Inverse solution method used to reconstruct EEG sources. We selected three of the widely used algorithms in the EEG/MEG community, namely (i) weighted minimum norm estimate (wMNE) (Fuchs et al., 1999; Lin et al., 2006), (ii) exact low-resolution brain electromagnetic tomography (eLORETA) (Pascual-Marqui 2007), and (iii) linearly constrained minimum variance (LCMV) beamformer (Van Veen et al., 1997).

(3) Functional connectivity metric. We selected two phase-based and two amplitude-based metrics: the phase-locking value (PLV) (Lachaux et al., 2000), phase-lag index (PLI) (Stam et al., 2007), and amplitude envelope correlation (AEC) with and without source leakage correction (Hipp et al., 2012; Colclough et al., 2016, 2015).

Finally, cortical networks obtained for each combination (number of electrodes × inverse solution × connectivity measure) were quantitatively compared to the reference simulated networks.

## 2. Materials and methods

A schematic diagram of the analysis pipeline is presented in Fig. 1.

### 2.1. Simulations

The simulated cortical networks (DMN and DAN) each included six regions based on the Desikan-Killiany atlas (Desikan et al., 2006) in terms of region parcellation. The DMN consisted of the right and left posterior cingulate cortex (PCC), medial orbitofrontal (MOF) gyrus, and inferior parietal lobe (IPL). Regarding the DAN, this network consisted of the right and left inferior parietal lobe (IPL), caudal middle frontal gyrus (cMFG), and superior parietal lobe (SPL). The choice of those regions was based on their frequent occurrence in previous resting-state studies (Allen et al., 2018; Allen et al., 2014; Damoiseaux et al., 2006; Fox and Raichle, 2007; Greicius et al., 2003; Baker et al., 2014; Shirer et al., 2012; Kabbara et al., 2017; Kabbara et al., 2021). Cortical-level activity was generated using a flexible neural mass model framework, named COALIA (Bensaid et al., 2019). This multi-population neural mass model enables the simulation of brain-scale electrophysiological activity while accounting for the macro- (between regions) and micro-circuitry (within a single region) of the brain, with one neural mass representing the local field potential of one Desikan-Killiany atlas region [for details, readers may refer to Bensaid et al. (2019)]. Activity in the alpha band ([8 – 12] Hz) was attributed to reference RSNs regions, while background activity was assigned to remaining cortical regions. To account for variability between simulated data segments, each “virtual subject” had different connectivity matrices provided to the model, while each epoch for the same subject had a different randomly generated input noise describing the average density of afferent action potentials from neighboring and distant populations (mean = 90 Hz, standard deviation = 30 Hz). More specifically, for each subject, a different fractional anisotropy matrix of the HCP dataset was used (Van Eschen et al., 2013), and the weights corresponding to an RSN-connection were modified and set to a value of (1 ± 20%). A corresponding scaling of the matrices followed in accordance with COALIA’s requisites and the type of each input matrix (inhibitory/excitatory). A total of 50 “virtual subjects”, 4 epochs per subject (i.e., 200 data segments) were simulated; with a duration of 40 s each and sampling rate of 2048 Hz. An example of simulated cortical signals is shown in Supplementary Materials, Fig. S1.

### 2.2. Forward problem

Scalp EEG signals were estimated from simulated cortical activity by solving the forward problem as follows:

$$X(t) = (x_1(t) \dots x_M(t))^T = G(s_1(t) \dots s_P(t))^T = G.S(t) \quad (1)$$

where $X(t)$ represents scalp EEG signals, $S(t)$ simulated cortical time series, and $G$ the $(M \times P)$ Gain (lead field) matrix. More specifically, $G$ quantifies the contribution of each cortical source to the generation of scalp signals by taking into account the geometrical and electrical characteristics of the head. Here, the gain matrix was computed using a realistic head model (Colin27 MRI template) using the boundary element method (BEM), implemented within the OpenMEEG package (Gramfort et al., 2010) in the Fieldtrip toolbox (Oostenveld et al., 2011).
2021). Interestingly, accurate identification of the so-called resting-state networks (RSNs) has been found to be primordial for both cognitive (Keller et al., 2015; Jockwitz et al., 2017; Shen et al., 2018; Alavash et al., 2015) and clinical studies (Sheline and Raichle 2013; Kabbara et al., 2018; Hassan et al., 2017; Gratton et al., 2018; Jimenez et al., 2019).

In this context, the EEG source connectivity method enables inferring functional brain networks at the cortical level with high temporal resolution (Hassan and Wendling, 2018). However, due to the complexity of the analysis workflow, many methodological questions related to the EEG source connectivity pipeline remain unanswered. In fact, each step of the analysis involves different choices that might significantly affect the resulting functional network. Hence, the effect of several factors needs to be investigated and quantified for more consistent and reliable use of this method.

One major factor in terms of impact on the EEG source connectivity analysis is the number of EEG electrodes (electrodes density). In fact, several studies investigated the effect of electrodes density on EEG source localization in simulations in the context of epilepsy (Song et al., 2015; Sohrabpour et al., 2015; Goran Lantz et al., 2003). More specifically, it has been shown that EEG electrodes density has a direct influence on the localization error: as intuitively expected, a higher number of electrodes is associated with a significant decrease in localization error (Goran Lantz et al., 2003; Song et al., 2015; Sohrabpour et al., 2015). In both studies by Goran Lantz et al. (2003) and Sohrabpour et al. (2015), a dramatic decrease in localization error occurred when increasing the number of electrodes from 32 to 64.

Another critical influencing factor in the EEG source connectivity pipeline is the algorithm chosen to solve the inverse problem, which is ill-posed due to its non-uniqueness and the instability of its solution (see Grech et al. (2008) for a review). Several studies quantifying the performance of different inverse methods, in simulated and experimental EEG/MEG data, concluded that the choice of the inverse method significantly influences source estimation results (Anzolin et al., 2019; Mahjoory et al., 2017; Hedrich et al., 2017; Bradley et al., 2016; Grova et al., 2006; Halder et al., 2019; Tait et al., 2021; Allouch et al., 2022). However, no consistent conclusions have been made regarding one method that would stand apart from the others in terms of performance, which can also be related to the analyzed conditions.

The choice of the functional connectivity metric is also a critical step. A wide range of measures are used in the field, and each differs in the aspect of the data being investigated (amplitude- vs phase-based measures / directional vs non-directional connectivity, prone/robust to source leakage), (see Friston (2011), Pereda et al. (2005), Cao et al. (2022) for a review), resulting in significant variability of performance and interpretations (Colclough et al., 2016; Wang et al., 2014; Wendling et al., 2009; Hassan et al., 2017; Allouch et al., 2022).

In this context, simulation studies are of utmost interest, since they provide an accessible ground truth that enables an objective evaluation of the methods/techniques under investigation. Here, we used a model for cortical dynamics, namely neural mass models (NMMs), to simulate resting state networks (RSNs), and more specifically the default mode network (DMN) and dorsal attentional network (DAN). Then, we quantified the effect of three key factors involved in the EEG source connectivity analysis:

(1) EEG channel density (from 19,32,64,128 to 256 channels).

(2) Inverse solution method used to reconstruct EEG sources. We selected three of the widely used algorithms in the EEG/MEG community, namely (i) weighted minimum norm estimate (wMNE) (Fuchs et al., 1999; Lin et al., 2006), (ii) exact low-resolution brain electromagnetic tomography (eLORETA) (Pascual-Marqui 2007), and (iii) linearly constrained minimum variance (LCMV) beamformer (Van Veen et al., 1997).

(3) Functional connectivity metric. We selected two phase-based and two amplitude-based metrics: the phase-locking value (PLV) (Lachaux et al., 2000), phase-lag index (PLI) (Stam et al., 2007), and amplitude envelope correlation (AEC) with and without source leakage correction (Hipp et al., 2012; Colclough et al., 2016, 2015).

Finally, cortical networks obtained for each combination (number of electrodes × inverse solution × connectivity measure) were quantitatively compared to the reference simulated networks.

## 2. Materials and methods

A schematic diagram of the analysis pipeline is presented in Fig. 1.

### 2.1. Simulations

The simulated cortical networks (DMN and DAN) each included six regions based on the Desikan-Killiany atlas (Desikan et al., 2006) in terms of region parcellation. The DMN consisted of the right and left posterior cingulate cortex (PCC), medial orbitofrontal (MOF) gyrus, and inferior parietal lobe (IPL). Regarding the DAN, this network consisted of the right and left inferior parietal lobe (IPL), caudal middle frontal gyrus (cMFG), and superior parietal lobe (SPL). The choice of those regions was based on their frequent occurrence in previous resting-state studies (Allen et al., 2018; Allen et al., 2014; Damoiseaux et al., 2006; Fox and Raichle, 2007; Greicius et al., 2003; Baker et al., 2014; Shirer et al., 2012; Kabbara et al., 2017; Kabbara et al., 2021). Cortical-level activity was generated using a flexible neural mass model framework, named COALIA (Bensaid et al., 2019). This multi-population neural mass model enables the simulation of brain-scale electrophysiological activity while accounting for the macro- (between regions) and micro-circuitry (within a single region) of the brain, with one neural mass representing the local field potential of one Desikan-Killiany atlas region [for details, readers may refer to Bensaid et al. (2019)]. Activity in the alpha band ([8 – 12] Hz) was attributed to reference RSNs regions, while background activity was assigned to remaining cortical regions. To account for variability between simulated data segments, each “virtual subject” had different connectivity matrices provided to the model, while each epoch for the same subject had a different randomly generated input noise describing the average density of afferent action potentials from neighboring and distant populations (mean = 90 Hz, standard deviation = 30 Hz). More specifically, for each subject, a different fractional anisotropy matrix of the HCP dataset was used (Van Eschen et al., 2013), and the weights corresponding to an RSN-connection were modified and set to a value of (1 ± 20%). A corresponding scaling of the matrices followed in accordance with COALIA’s requisites and the type of each input matrix (inhibitory/excitatory). A total of 50 “virtual subjects”, 4 epochs per subject (i.e., 200 data segments) were simulated; with a duration of 40 s each and sampling rate of 2048 Hz. An example of simulated cortical signals is shown in Supplementary Materials, Fig. S1.

### 2.2. Forward problem

Scalp EEG signals were estimated from simulated cortical activity by solving the forward problem as follows:

$$X(t) = (x_1(t) \dots x_M(t))^T = G(s_1(t) \dots s_P(t))^T = G.S(t) \quad (1)$$

where $X(t)$ represents scalp EEG signals, $S(t)$ simulated cortical time series, and $G$ the $(M \times P)$ Gain (lead field) matrix. More specifically, $G$ quantifies the contribution of each cortical source to the generation of scalp signals by taking into account the geometrical and electrical characteristics of the head. Here, the gain matrix was computed using a realistic head model (Colin27 MRI template) using the boundary element method (BEM), implemented within the OpenMEEG package (Gramfort et al., 2010) in the Fieldtrip toolbox (Oostenveld et al., 2011).
<figure><img src="image_3.png" /><figcaption>Fig. 1. Processing pipeline used in the present study. (A) Generation of reference, ground-truth cortical-level electrophysiological signals using a realistic model of neuronal activity. (B) Projection of source-level activity on scalp-level EEG sensors by using an anatomically accurate head model. (C) Network reconstruction using different methods of source reconstruction based on scalp-level signals, and of statistical coupling between signals (functional connectivity metrics). (D) Quantification of the level of matching between estimated and reference networks. DMN - default mode network. DAN - Dorsal attentional network. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

The sensor space was defined based on the GSN HydroCel electrodes configuration (EGI, Electrical geodesic Inc) with 256, 128, 64, and 32 channels, as well as the international 10–20 system with 19 channels. The lead field matrix used for solving the forward problem, as well as the inverse problem, was constrained to 66 lead field vectors representing the contribution of sources located at the centroid of the regions of interest defined by the Desikan-Killiany atlas (Desikan et al., 2006) (right and left insula were excluded, leaving 66 regions of interest).

## 2.3. Inverse problem

The first step to reconstruct cortical networks was to estimate the dynamics of cortical sources from scalp EEG data, i.e., determining the position, orientation, and magnitude of dipolar sources $\hat{S}(t)$. Cortical sources were located at the centroids of Desikan-Killiany regions and oriented normally to the cortical sheet. Thus, the inverse problem was reduced to computing the magnitude of dipolar sources $\hat{S}(t)$:

$$ \hat{S}(t) = W \cdot X(t) \qquad (2) $$

Several algorithms have been proposed to solve this problem, and estimate W based on different assumptions related to the spa-

tiotemporal properties of sources and regularization constraints (see Awan et al. (2019), Grech et al. (2008), Baillet et al. (2001) for a review). Inverse solutions can be classified into two major families: minimum norm estimates and beamformers. The former reconstructs all sources simultaneously by minimizing the difference between the data $X(t)$ and predicted data $G$. $\hat{S}(t)$, while beamformers take an adaptive spatial-filtering approach in which each source is scanned independently. In this study, we focused on three commonly used solutions in EEG source reconstruction, as detailed below.

### 2.3.1. Weighted minimum norm estimate (wMNE)

wMNE (Lin et al., 2006; Fuchs et al., 1999) is a derivative of the minimum norm estimate (MNE) (Hämäläinen and Ilmoniemi, 1994), which proposes a solution that fits the measurements with a least square error. However, wMNE compensates further for the tendency of MNE to favor weak and surface sources:

$$ W_{wMNE} = BG^T (GB^T + \lambda C)^{-1} \qquad (3) $$

where $\lambda$ is the regularization parameter and $C$ is the noise covariance matrix (set to the identity matrix in our case). The matrix $B$ is a diagonal matrix built from matrix $G$ with non-zero terms inversely proportional
<figure><img src="image_5.png" /><figcaption>Fig. 1. Processing pipeline used in the present study. (A) Generation of reference, ground-truth cortical-level electrophysiological signals using a realistic model of neuronal activity. (B) Projection of source-level activity on scalp-level EEG sensors by using an anatomically accurate head model. (C) Network reconstruction using different methods of source reconstruction based on scalp-level signals, and of statistical coupling between signals (functional connectivity metrics). (D) Quantification of the level of matching between estimated and reference networks. DMN - default mode network. DAN - Dorsal attentional network. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

The sensor space was defined based on the GSN HydroCel electrodes configuration (EGI, Electrical geodesic Inc) with 256, 128, 64, and 32 channels, as well as the international 10–20 system with 19 channels. The lead field matrix used for solving the forward problem, as well as the inverse problem, was constrained to 66 lead field vectors representing the contribution of sources located at the centroid of the regions of interest defined by the Desikan-Killiany atlas (Desikan et al., 2006) (right and left insula were excluded, leaving 66 regions of interest).

## 2.3. Inverse problem

The first step to reconstruct cortical networks was to estimate the dynamics of cortical sources from scalp EEG data, i.e., determining the position, orientation, and magnitude of dipolar sources $\hat{S}(t)$. Cortical sources were located at the centroids of Desikan-Killiany regions and oriented normally to the cortical sheet. Thus, the inverse problem was reduced to computing the magnitude of dipolar sources $\hat{S}(t)$:

$$ \hat{S}(t) = W \cdot X(t) \qquad (2) $$

Several algorithms have been proposed to solve this problem, and estimate $W$ based on different assumptions related to the spa-

tiotemporal properties of sources and regularization constraints (see Awan et al. (2019), Grech et al. (2008), Baillet et al. (2001) for a review). Inverse solutions can be classified into two major families: minimum norm estimates and beamformers. The former reconstructs all sources simultaneously by minimizing the difference between the data $X(t)$ and predicted data $G$. $\hat{S}(t)$, while beamformers take an adaptive spatial-filtering approach in which each source is scanned independently. In this study, we focused on three commonly used solutions in EEG source reconstruction, as detailed below.

### 2.3.1. Weighted minimum norm estimate (wMNE)

wMNE (Lin et al., 2006; Fuchs et al., 1999) is a derivative of the minimum norm estimate (MNE) (Hämäläinen and Ilmoniemi, 1994), which proposes a solution that fits the measurements with a least square error. However, wMNE compensates further for the tendency of MNE to favor weak and surface sources:

$$ W_{wMNE} = BG^T (GBG^T + \lambda C)^{-1} \qquad (3) $$

where $\lambda$ is the regularization parameter and $C$ is the noise covariance matrix (set to the identity matrix in our case). The matrix $B$ is a diagonal matrix built from matrix $G$ with non-zero terms inversely proportional
to the norm of lead field vectors. This matrix adjusts the properties of the solution by reducing the bias inherent to the standard MNE solution:

$$B_{ij} = \begin{cases} (G_i^T G_i)^{\frac{1}{2}} & \text{if } i = j; \\ 0 & \text{if } i \neq j; \end{cases} \qquad (4)$$

Regarding the regularization parameter $\lambda$, we used the recommended default value ($\frac{1}{SNR}$; SNR = 3) included in the Brainstorm toolbox (Tadel et al., 2011).

### 2.3.2. Exact low-resolution brain electromagnetic tomography (eLORETA)

eLORETA belongs to the family of weighted minimum norm inverse solutions. However, in addition to compensating for depth bias, it also has exact-zero error localization in the presence of measurement and structured biological noise (Pascual-Marqui, 2007):

$$B_{ij} = \begin{cases} (G_i^T (G_i B G_i^T + \lambda C)^{-1})^{\frac{1}{2}} & \text{if } i = j; \\ 0 & \text{if } i \neq j; \end{cases} \qquad (5)$$

Regarding the regularization parameter $\lambda$, we used the default value (0.05) of the Fieldtrip toolbox (Oostenveld et al., 2011).

### 2.3.3. Linearly constrained minimum-variance (LCMV) beamformer

The LCMV beamformer (Van Veen et al., 1997) takes an adaptive spatial-filtering approach and estimates the activity for a source at a given location while simultaneously suppressing contributions from all other sources and noise captured in the data covariance.

$$W_{LCMV} = ((G^T \cdot (C + \lambda \cdot I)^{-1}) \cdot G)^{-1} \cdot (G^T \cdot (C + \lambda \cdot I)^{-1}) \qquad (6)$$

The regularization parameter $\lambda$ was set to 0.05.

## 2.4. Functional connectivity

Following the reconstruction of cortical dynamics, functional connectivity was assessed to estimate cortical networks. Two distinct approaches can be adopted to compute functional connectivity: phase-based and amplitude-based techniques. We tested here four different connectivity measures:

### 2.4.1. Phase-locking value (PLV)

For two signals $x(t)$ and $y(t)$, the phase-locking value (Lachaux et al., 2000) is defined as:

$$PLV = \left| E \left\{ e^{i(\phi_x(t) - \phi_y(t))} \right\} \right| \qquad (7)$$

where $E\{\cdot\}$ is the expected value operator, and $\phi(t)$ is the instantaneous phase derived from the Hilbert transform. PLV was computed over consecutive non-overlapping sliding windows, with the length of the window $\delta$ set to $\frac{1}{\text{central frequency}}$ as recommended in (Lachaux et al., 2000), and where 6 is the number of cycles in a given frequency band. Thus, $\delta$ equals 600 ms in the considered alpha band ([8 - 12] Hz) where the central frequency was equal to 10 Hz. PLV values were then averaged over all sliding windows.

### 2.4.2. Phase-lag index (PLI)

The PLI originally proposed in Stam et al. (2007) is a measure of the asymmetry of the distribution of phase differences between two signals. It aims at overcoming source leakage by discarding phase differences centered around 0 and $\pi$, i.e., removing zero-lag connections. For two signals $x(t)$ and $y(t)$, PLI is defined as follows:

$$PLI = \left| E\{\text{sign}(\phi_x(t) - \phi_y(t))\} \right| \qquad (8)$$

where $E\{\cdot\}$ is the expected value operator, and $\phi(t)$ the instantaneous phase derived from the Hilbert transform. Similar to PLV computation, PLI was assessed over consecutive non-overlapping sliding windows (600 ms). PLI values were then averaged over all sliding windows.

### 2.4.3. Amplitude envelope correlation (AEC)

AEC was computed as the Pearson correlation between signals' envelopes derived from the Hilbert transform (Hipp et al., 2012; Brookes et al., 2011). Similar to PLV and PLI computation, a sliding window approach was adopted. Based on (O'Neill et al., 2017), the window length was set to 6 s with an overlap of 0.5 s.

### 2.4.4. Amplitude envelope correlation with source leakage correction (AEC*)

Zero-lag signal overlaps caused by spatial leakage were removed using a multivariate symmetric orthogonalization approach (Colclough et al., 2015). Briefly, the closest orthonormal matrix to the uncorrected regional time courses is first computed; then the magnitudes of the orthogonalized vectors are adjusted iteratively to minimize the least-squares distances between corrected and uncorrected signals. Orthogonalization was applied over the entire data segment. Following the orthogonalization procedure, AEC was computed as previously described in Section 2.4.3.

## 2.5. Results quantification and statistical testing

We used three metrics, detailed below, to assess the performance of different channel densities and the tested inverse solutions and connectivity measures.

### 2.5.1. Pearson correlation

First, the Pearson correlation between the reference and reconstructed connectivity matrices was computed as a global measure of the similarity between networks. Matrices were not thresholded nor binned, and the Pearson correlation was computed between the vector representation of the upper triangular elements.

### 2.5.2. Closeness accuracy

Second, closeness accuracy (CA), defined as follows, was used for a node-wise comparison between the reference and reconstructed networks:

$$CA = \frac{1}{1 + AD} \qquad (9)$$

where AD is the average distance (in cm) between the reference and reconstructed networks given by:

$$AD = \frac{\sum_k d(N_k, N_v)}{M} \quad k \in [1, M]; v \in [1, W] \qquad (10)$$

where $d(N_k, N_v)$ denotes the Euclidean distance between the node $N_k$ in the reconstructed network, and the nearest node $N_v$ in the reference simulated network. M and W represent the total number of nodes detected in reconstructed and reference networks, respectively. Before computing the CA, networks were thresholded by keeping the edges with the highest 0.7% weight values, the choice of this proportion corresponding to the number of edges in the simulated RSNs (30 edges).

### 2.5.3. Edge contribution

Finally, we investigated the contribution of individual edges to the correlation values obtained between the reference and reconstructed networks (Finn et al., 2015; Colclough et al., 2016). The set of edges of reference and reconstructed (non-thresholded) networks were first z-score normalized (mean = 0, std = 1). Second, the edge contribution $\varphi$ was calculated as the element-wise product between the two normalized edge vectors (upper triangular elements). Suppose that $[X_i^{\text{Ref}}]$ and $[X_i^{\text{Rec}}]$ are the set of edges of the reference and reconstructed networks for an epoch $i$, after z-score normalization, then:

$$\varphi_i(e) = X_i^{\text{Ref}}(e) * X_i^{\text{Rec}}(e), e \in [1, E] \qquad (11)$$

where $i$ is the epoch index, $e$ the edge index, and $E$ the total number of edges. The average of $\varphi_i$ across all subjects and all epochs is denoted
to the norm of lead field vectors. This matrix adjusts the properties of the solution by reducing the bias inherent to the standard MNE solution:

$$B_{ij} = \begin{cases} (G_i^T G_i)^{\frac{1}{2}} & \text{if } i = j; \\ 0 & \text{if } i \neq j; \end{cases} \qquad (4)$$

Regarding the regularization parameter $\lambda$, we used the recommended default value ($\frac{1}{SNR}$; $SNR = 3$) included in the Brainstorm toolbox (Tadel et al., 2011).

### 2.3.2. Exact low-resolution brain electromagnetic tomography (eLORETA)

eLORETA belongs to the family of weighted minimum norm inverse solutions. However, in addition to compensating for depth bias, it also has exact-zero error localization in the presence of measurement and structured biological noise (Pascual-Marqui, 2007):

$$B_{ij} = \begin{cases} (G_i^T (G_i B G_i^T + \lambda C)^{-1})^{\frac{1}{2}} & \text{if } i = j; \\ 0 & \text{if } i \neq j; \end{cases} \qquad (5)$$

Regarding the regularization parameter $\lambda$, we used the default value (0.05) of the Fieldtrip toolbox (Oostenveld et al., 2011).

### 2.3.3. Linearly constrained minimum-variance (LCMV) beamformer

The LCMV beamformer (Van Veen et al., 1997) takes an adaptive spatial-filtering approach and estimates the activity for a source at a given location while simultaneously suppressing contributions from all other sources and noise captured in the data covariance.

$$W_{LCMV} = ((G^T \cdot (C + \lambda \cdot I)^{-1}) \cdot G)^{-1} \cdot (G^T \cdot (C + \lambda \cdot I)^{-1}) \qquad (6)$$

The regularization parameter $\lambda$ was set to 0.05.

## 2.4. Functional connectivity

Following the reconstruction of cortical dynamics, functional connectivity was assessed to estimate cortical networks. Two distinct approaches can be adopted to compute functional connectivity: phase-based and amplitude-based techniques. We tested here four different connectivity measures:

### 2.4.1. Phase-locking value (PLV)

For two signals $x(t)$ and $y(t)$, the phase-locking value (Lachaux et al., 2000) is defined as:

$$PLV = \left| E \left\{ e^{i(\phi_x(t) - \phi_y(t))} \right\} \right| \qquad (7)$$

where $E\{\cdot\}$ is the expected value operator, and $\phi(t)$ is the instantaneous phase derived from the Hilbert transform. PLV was computed over consecutive non-overlapping sliding windows, with the length of the window $\delta$ set to $\frac{1}{\text{central frequency}}$ as recommended in (Lachaux et al., 2000), and where 6 is the number of cycles in a given frequency band. Thus, $\delta$ equals 600 ms in the considered alpha band ([8 - 12] Hz) where the central frequency was equal to 10 Hz. PLV values were then averaged over all sliding windows.

### 2.4.2. Phase-lag index (PLI)

The PLI originally proposed in Stam et al. (2007) is a measure of the asymmetry of the distribution of phase differences between two signals. It aims at overcoming source leakage by discarding phase differences centered around 0 and $\pi$, i.e., removing zero-lag connections. For two signals $x(t)$ and $y(t)$, PLI is defined as follows:

$$PLI = \left| E\{\text{sign}(\phi_x(t) - \phi_y(t))\} \right| \qquad (8)$$

where $E\{\cdot\}$ is the expected value operator, and $\phi(t)$ the instantaneous phase derived from the Hilbert transform. Similar to PLV computation, PLI was assessed over consecutive non-overlapping sliding windows (600 ms). PLI values were then averaged over all sliding windows.

### 2.4.3. Amplitude envelope correlation (AEC)

AEC was computed as the Pearson correlation between signals' envelopes derived from the Hilbert transform (Hipp et al., 2012; Brookes et al., 2011). Similar to PLV and PLI computation, a sliding window approach was adopted. Based on (O'Neill et al., 2017), the window length was set to 6 s with an overlap of 0.5 s.

### 2.4.4. Amplitude envelope correlation with source leakage correction (AEC*)

Zero-lag signal overlaps caused by spatial leakage were removed using a multivariate symmetric orthogonalization approach (Colclough et al., 2015). Briefly, the closest orthonormal matrix to the uncorrected regional time courses is first computed; then the magnitudes of the orthogonalized vectors are adjusted iteratively to minimize the least-squares distances between corrected and uncorrected signals. Orthogonalization was applied over the entire data segment. Following the orthogonalization procedure, AEC was computed as previously described in Section 2.4.3.

## 2.5. Results quantification and statistical testing

We used three metrics, detailed below, to assess the performance of different channel densities and the tested inverse solutions and connectivity measures.

### 2.5.1. Pearson correlation

First, the Pearson correlation between the reference and reconstructed connectivity matrices was computed as a global measure of the similarity between networks. Matrices were not thresholded nor binned, and the Pearson correlation was computed between the vector representation of the upper triangular elements.

### 2.5.2. Closeness accuracy

Second, closeness accuracy (CA), defined as follows, was used for a node-wise comparison between the reference and reconstructed networks:

$$CA = \frac{1}{1 + AD} \qquad (9)$$

where $AD$ is the average distance (in cm) between the reference and reconstructed networks given by:

$$AD = \frac{\sum_k d(N_k, N_v)}{M} \quad k \in [1, M]; \quad v \in [1, W] \qquad (10)$$

where $d(N_k, N_v)$ denotes the Euclidean distance between the node $N_k$ in the reconstructed network, and the nearest node $N_v$ in the reference simulated network. $M$ and $W$ represent the total number of nodes detected in reconstructed and reference networks, respectively. Before computing the CA, networks were thresholded by keeping the edges with the highest 0.7% weight values, the choice of this proportion corresponding to the number of edges in the simulated RSNs (30 edges).

### 2.5.3. Edge contribution

Finally, we investigated the contribution of individual edges to the correlation values obtained between the reference and reconstructed networks (Finn et al., 2015; Colclough et al., 2016). The set of edges of reference and reconstructed (non-thresholded) networks were first z-score normalized (mean = 0, std = 1). Second, the edge contribution $\varphi$ was calculated as the element-wise product between the two normalized edge vectors (upper triangular elements). Suppose that $[X_i^{\text{Ref}}]$ and $[X_i^{\text{Rec}}]$ are the set of edges of the reference and reconstructed networks for an epoch $i$, after z-score normalization, then:

$$\varphi_i(e) = X_i^{\text{Ref}}(e) * X_i^{\text{Rec}}(e), e \in [1, E] \qquad (11)$$

where $i$ is the epoch index, $e$ the edge index, and $E$ the total number of edges. The average of $\varphi_i$ across all subjects and all epochs is denoted
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Violin plots of the Pearson correlation values computed between reference and reconstructed DMNs for all electrode montages and inverse methods/connectivity metrics combinations. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC\* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

φ. High positive φ values correspond to the edges that are consistent within- and between- subjects.

All statistical analyses were performed using R (R Core Team, 2022). Linear mixed models implemented using the `{lme4}` package (Bates et al., 2015) were used to assess the effect of (i) number of electrodes, (ii) inverse solution, and (iii) functional connectivity metric as fixed effects on correlation and closeness accuracy. The “subject” was added as a random intercept to account for simulation-related variability (we had 50 subjects and 4 epochs per subject). Assumptions of normality and homoscedasticity of the model’s residuals were graphically checked. We used the same model for correlation and closeness accuracy values with the difference that closeness accuracy was inverse-transformed because of better compliance with the model’s assumptions in that case. Calculation of the significance of fixed effects was performed with *F* tests using the *Anova* function of the `{car}` package (Fox and Weisberg, 2019), and *post-hoc* analyses were performed using *z*-tests with the *glht* function of the `{multcomp}` package (Hothorn et al., 2008) that provides corrected *p*-values. Marginal and conditional *R*² were calculated with the `{MuMin}` package (Barton, 2009) to estimate the variance explained by models. The significance threshold was set at its usual value of *p* = 0.05.

## 3. Results

Fig. 2 and Fig. 3 present, respectively, the distributions of correlation and closeness accuracy values, for all tested conditions (number of channels × inverse solution × connectivity measure) (results obtained for DAN can be found in the Supplementary Materials, Figs. S2, S3). Since the orthogonalization procedure used here is limited by the rank of the data matrix (Colclough et al., 2015), it can only be applied to 128- and 256-electrode configurations. Therefore, the statistical tests were repeated twice:

* Case 1: comparing all electrodes configurations (19, 32, 64, 128, 256), inverse solutions (eLORETA, LCMV, wMNE), and three connectivity measures (PLV, AEC, PLI) where no orthogonalization is applied.

* Case 2: comparing two electrode configurations (128, 256), all inverse solutions (eLORETA, LCMV, wMNE), and all connectivity metrics (PLV, AEC, PLI, AEC\*).

Detailed statistical results are reported in the Supplementary Materials.

### 3.1. Pearson correlation

For all inverse methods and connectivity metrics, the trend in correlation values was similar as the number of electrodes increased, for both DMN and DAN: a greater number of channels was associated with higher correlation values. Pearson correlation results demonstrate a drastic effect of increasing the number of electrodes on enhancing the accuracy of reconstructed networks, regardless of the chosen inverse solution/connectivity measure. Statistical tests comparing all channels configurations, inverse solutions, and connectivity measures (PLV, AEC, PLI) (case 1) revealed that the effect of channel density was significant (*F*(4, 8906)=23,910, *p*<0.0001; full model marginal *R*²=0.97; conditional *R*²=0.97) as well as all pairwise comparisons (*p* < 0.0001). Similarly, comparing all inverse solutions and connectivity measures (PLV, AEC, PLI, AEC\*) at 128 and 256 channels (case 2) revealed a significant effect of electrodes density (*F*(1, 4727) = 4393, *p* < 0.0001; full model marginal *R*² = 0.99; conditional *R*² = 0.99). Significant effects of the inverse method (case 1: *F*(2, 8906) = 3192, *p* < 0.0001; case 2: *F*(2, 4727) = 8604, *p* < 0.0001) and functional connectivity measure (case 1: *F*(2, 8906) = 103, 998, *p* < 0.0001; case 2: *F*(3, 4727) = 175, 314, *p* < 0.0001) were also identified. Moreover, reconstructed networks differed as a function of the inverse solution/connectivity metric combination (case 1: *F*(4, 8906) = 344, *p* < 0.0001; case 2: *F*(6, 4727) = 3341, *p* < 0.0001). Post-hoc analysis showed significant differences between inverse solution algorithms (case 1: *p* < 0.001, except for wMNE vs LCMV: *p* = 0.65; case 2: *p* < 0.001 except for wMNE vs LCMV: *p* = 0.60), as well as connectivity measures (*p* < 0.001). Pairwise comparisons of the different combinations of inverse solution and functional connectivity methods also revealed significant differences between combinations (case 1: *p* < 0.001 except for wMNE/AEC vs LCMV/AEC (*p* = 0.99), wMNE/PLI vs LCMV/PLI (*p* = 0.40); case 2: *p* < 0.001 except for wMNE/AEC vs LCMV/AEC (*p* = 1) and wMNE/PLI vs LCMV/PLI (*p* = 0.41). At high sensor densities, AEC computed between orthogonalized time series had the best reconstruction accuracy. PLI, which also discard zero-lag connections, came second with all three inverse solutions. Regarding measures that do not compensate for spatial leakage, PLV had moderate performance, whereas AEC values were systematically low, regardless of the inverse method and number of electrodes.
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Violin plots of the Pearson correlation values computed between reference and reconstructed DMNs for all electrode montages and inverse methods/connectivity metrics combinations. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC\* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

φ. High positive φ values correspond to the edges that are consistent within- and between- subjects.

All statistical analyses were performed using R (R Core Team, 2022). Linear mixed models implemented using the `{lme4}` package (Bates et al., 2015) were used to assess the effect of (i) number of electrodes, (ii) inverse solution, and (iii) functional connectivity metric as fixed effects on correlation and closeness accuracy. The “subject” was added as a random intercept to account for simulation-related variability (we had 50 subjects and 4 epochs per subject). Assumptions of normality and homoscedasticity of the model’s residuals were graphically checked. We used the same model for correlation and closeness accuracy values with the difference that closeness accuracy was inverse-transformed because of better compliance with the model’s assumptions in that case. Calculation of the significance of fixed effects was performed with *F* tests using the *Anova* function of the `{car}` package (Fox and Weisberg, 2019), and *post-hoc* analyses were performed using *z*-tests with the *glht* function of the `{multcomp}` package (Hothorn et al., 2008) that provides corrected *p*-values. Marginal and conditional *R*² were calculated with the `{MuMin}` package (Barton, 2009) to estimate the variance explained by models. The significance threshold was set at its usual value of *p* = 0.05.

### 3. Results

Fig. 2 and Fig. 3 present, respectively, the distributions of correlation and closeness accuracy values, for all tested conditions (number of channels × inverse solution × connectivity measure) (results obtained for DAN can be found in the Supplementary Materials, Figs. S2, S3). Since the orthogonalization procedure used here is limited by the rank of the data matrix (Colclough et al., 2015), it can only be applied to 128- and 256-electrode configurations. Therefore, the statistical tests were repeated twice:

* Case 1: comparing all electrodes configurations (19, 32, 64, 128, 256), inverse solutions (eLORETA, LCMV, wMNE), and three connectivity measures (PLV, AEC, PLI) where no orthogonalization is applied.

* Case 2: comparing two electrode configurations (128, 256), all inverse solutions (eLORETA, LCMV, wMNE), and all connectivity metrics (PLV, AEC, PLI, AEC\*).

Detailed statistical results are reported in the Supplementary Materials.

#### 3.1. Pearson correlation

For all inverse methods and connectivity metrics, the trend in correlation values was similar as the number of electrodes increased, for both DMN and DAN: a greater number of channels was associated with higher correlation values. Pearson correlation results demonstrate a drastic effect of increasing the number of electrodes on enhancing the accuracy of reconstructed networks, regardless of the chosen inverse solution/connectivity measure. Statistical tests comparing all channels configurations, inverse solutions, and connectivity measures (PLV, AEC, PLI) (case 1) revealed that the effect of channel density was significant (*F*(4, 8906)=23,910, *p*<0.0001; full model marginal *R*²=0.97; conditional *R*²=0.97) as well as all pairwise comparisons (*p* < 0.0001). Similarly, comparing all inverse solutions and connectivity measures (PLV, AEC, PLI, AEC\*) at 128 and 256 channels (case 2) revealed a significant effect of electrodes density (*F*(1, 4727) = 4393, *p* < 0.0001; full model marginal *R*² = 0.99; conditional *R*² = 0.99). Significant effects of the inverse method (case 1: *F*(2, 8906) = 3192, *p* < 0.0001; case 2: *F*(2, 4727) = 8604, *p* < 0.0001) and functional connectivity measure (case 1: *F*(2, 8906) = 103, 998, *p* < 0.0001; case 2: *F*(3, 4727) = 175, 314, *p* < 0.0001) were also identified. Moreover, reconstructed networks differed as a function of the inverse solution/connectivity metric combination (case 1: *F*(4, 8906) = 344, *p* < 0.0001; case 2: *F*(6, 4727) = 3341, *p* < 0.0001). *Post-hoc* analysis showed significant differences between inverse solution algorithms (case 1: *p* < 0.001, except for wMNE vs LCMV: *p* = 0.65; case 2: *p* < 0.001 except for wMNE vs LCMV: *p* = 0.60), as well as connectivity measures (*p* < 0.001). Pairwise comparisons of the different combinations of inverse solution and functional connectivity methods also revealed significant differences between combinations (case 1: *p* < 0.001 except for wMNE/AEC vs LCMV/AEC (*p* = 0.99), wMNE/PLI vs LCMV/PLI (*p* = 0.40); case 2: *p* < 0.001 except for wMNE/AEC vs LCMV/AEC (*p* = 1) and wMNE/PLI vs LCMV/PLI (*p* = 0.41). At high sensor densities, AEC computed between orthogonalized time series had the best reconstruction accuracy. PLI, which also discard zero-lag connections, came second with all three inverse solutions. Regarding measures that do not compensate for spatial leakage, PLV had moderate performance, whereas AEC values were systematically low, regardless of the inverse method and number of electrodes.
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Violin plots of closeness accuracy values computed between the reference and reconstructed DMNs for all electrode montages and inverse methods/connectivity metrics combinations. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC\* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

## 3.2. Closeness accuracy

Differences in closeness accuracy values between evaluated electrode configurations were not as clear as with the Pearson correlation values, except for PLI and AEC\*, which were associated with higher between-epochs variability. We believe that this high variability is a result of the mathematical formulation of the closeness accuracy metric (See Results quantification Section in the Supplementary Materials). Similar to correlation results, closeness accuracy values were significantly affected by electrodes density (case 1: $F(4, 8906) = 775, p < 0.001$; full model marginal $R^2 = 0.77$; conditional $R^2 = 0.78$; case 2: $F(1, 4727) = 51, p < 0.001$; full model marginal $R^2 = 0.88$; conditional $R^2 = 0.88$), inverse method (case 1: $F(2, 8906) = 318, p < 0.001$; case 2: $F(2, 4727) = 334, p < 0.001$), connectivity metric (case 1: $F(2, 8906) = 12,574, p < 0.001$; case 2: $F(3, 4727) = 1156, p < 0.001$), and inverse method/connectivity metric combination (case 1: $F(4, 8906) = 120, p < 0.001$; case 2: $F(6, 4727) = 92, p < 0.001$). Pairwise comparisons of channel density showed significant differences ($p < 0.01$ except for (256 vs 128 ($p = 0.73$) and 64 vs 256 ($p = 0.1$)). When comparing all electrodes configurations, inverse solutions, and (PLV, AEC, PLI), post-hoc analysis showed significant differences between LCMV and eLORETA ($p < 0.001$), as well as between all connectivity metrics ($p < 0.001$). Regarding the inverse solution/connectivity measure combination, differences were not all significant. When comparing results obtained with 128 and 256 electrodes, all inverse solutions and connectivity metrics, significant differences were obtained between LCMV and eLORETA ($p < 0.01$), all connectivity measures ($p < 0.001$), and most inverse solution/connectivity measure combinations.

## 3.3. Edge contribution

We thresholded the averaged edge contribution matrices and kept only the 30 edges (= number of edges in the reference RSNs) having the highest contribution values. The percentile of edges (among those having the highest contribution value) belonging to DMN (DAN) are shown in Fig. 4 (Fig. S4 in Supplementary Materials). With PLV and AEC, the percentage of edges within the DMN increased at higher channel density. With AEC, the percentage did not exceed 50%. In contrast, when using PLI and corrected AEC, the majority of edges having the highest contribution were within the DMN regardless of the number of channels.

Those results were slightly different in the DAN, notably for PLV and AEC, both of which having larger percentages in DAN as compared to DMN. Interestingly, these measures have one thing in common: they are both prone to source leakage. Therefore, one explanation for the observed results could be that source leakage was prominent in DMN, which has a less distributed topology; as a result, spurious connections were detected by PLV and AEC and were consistent across subjects and epochs, resulting in high mean edge contribution values. In contrast, a distributed network topology in DAN led to a lower source leakage effect. Thus, PLV and AEC had higher percentage values, as shown in Fig. S4.

# 4. Discussion

The EEG/MEG source connectivity technique has gained increased interest due to its ability to reconstruct functional networks in the cortical space with a high temporal resolution. However, as aforementioned, there is still no agreement to date over a unified detailed EEG source connectivity pipeline, and several involved parameters/methods require cautious tuning. Such analytical variability is potentially problematic when associated with substantial variability in the outcomes. Here, our objective was to evaluate some parameters related to the EEG source connectivity analysis in the context of simulated RSNs and investigate results variability caused by different choices in the analysis pipeline. We mainly focused on the effect of electrodes density, inverse method, and functional connectivity metric. In a recent study, we tested the effect of electrodes density and tested wMNE and eLORETA along with PLV and wPLI in the context of simulated epileptiform activity (Allouch et al., 2022). In this specific case, epileptiform signals had a sufficiently high signal-to-noise-ratio (SNR) to be distinguished from background noise (Wa 1983; Iwasaki et al., 2005), which might have facilitated the identification of underlying networks. However, the last two decades have witnessed an increased interest in resting-state studies where the significantly lower SNR could be in favor of different methods/parameter tuning in EEG source connectivity analysis. Therefore, in the present study, we used neural mass models to simulate resting-state brain activity (DMN and DAN) in the cortical space and derived the corresponding scalp EEG signals by solving the forward problem. Then, we reconstructed the corresponding cortical networks. This pipeline was repeated using different electrode densities (19, 32, 64, 128, and 256 channels) to solve the forward problem, three algorithms to solve the inverse prob-
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Violin plots of closeness accuracy values computed between the reference and reconstructed DMNs for all electrode montages and inverse methods/connectivity metrics combinations. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC\* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

## 3.2. Closeness accuracy

Differences in closeness accuracy values between evaluated electrode configurations were not as clear as with the Pearson correlation values, except for PLI and AEC\*, which were associated with higher between-epochs variability. We believe that this high variability is a result of the mathematical formulation of the closeness accuracy metric (See Results quantification Section in the Supplementary Materials). Similar to correlation results, closeness accuracy values were significantly affected by electrodes density (case 1: $F(4, 8906) = 775, p < 0.001$; full model marginal $R^2 = 0.77$; conditional $R^2 = 0.78$; case 2: $F(1, 4727) = 51, p < 0.001$; full model marginal $R^2 = 0.88$; conditional $R^2 = 0.88$), inverse method (case 1: $F(2, 8906) = 318, p < 0.001$; case 2: $F(2, 4727) = 334, p < 0.001$), connectivity metric (case 1: $F(2, 8906) = 12,574, p < 0.001$; case 2: $F(3, 4727) = 1156, p < 0.001$), and inverse method/connectivity metric combination (case 1: $F(4, 8906) = 120, p < 0.001$; case 2: $F(6, 4727) = 92, p < 0.001$). Pairwise comparisons of channel density showed significant differences ($p < 0.01$ except for (256 vs 128 ($p = 0.73$) and 64 vs 256 ($p = 0.1$)). When comparing all electrodes configurations, inverse solutions, and (PLV, AEC, PLI), post-hoc analysis showed significant differences between LCMV and eLORETA ($p < 0.001$), as well as between all connectivity metrics ($p < 0.001$). Regarding the inverse solution/connectivity measure combination, differences were not all significant. When comparing results obtained with 128 and 256 electrodes, all inverse solutions and connectivity metrics, significant differences were obtained between LCMV and eLORETA ($p < 0.01$), all connectivity measures ($p < 0.001$), and most inverse solution/connectivity measure combinations.

## 3.3. Edge contribution

We thresholded the averaged edge contribution matrices and kept only the 30 edges (= number of edges in the reference RSNs) having the highest contribution values. The percentile of edges (among those having the highest contribution value) belonging to DMN (DAN) are shown in Fig. 4 (Fig. S4 in Supplementary Materials). With PLV and AEC, the percentage of edges within the DMN increased at higher channel density. With AEC, the percentage did not exceed 50%. In contrast, when using PLI and corrected AEC, the majority of edges having the highest contribution were within the DMN regardless of the number of channels.

Those results were slightly different in the DAN, notably for PLV and AEC, both of which having larger percentages in DAN as compared to DMN. Interestingly, these measures have one thing in common: they are both prone to source leakage. Therefore, one explanation for the observed results could be that source leakage was prominent in DMN, which has a less distributed topology; as a result, spurious connections were detected by PLV and AEC and were consistent across subjects and epochs, resulting in high mean edge contribution values. In contrast, a distributed network topology in DAN led to a lower source leakage effect. Thus, PLV and AEC had higher percentage values, as shown in Fig. S4.

# 4. Discussion

The EEG/MEG source connectivity technique has gained increased interest due to its ability to reconstruct functional networks in the cortical space with a high temporal resolution. However, as aforementioned, there is still no agreement to date over a unified detailed EEG source connectivity pipeline, and several involved parameters/methods require cautious tuning. Such analytical variability is potentially problematic when associated with substantial variability in the outcomes. Here, our objective was to evaluate some parameters related to the EEG source connectivity analysis in the context of simulated RSNs and investigate results variability caused by different choices in the analysis pipeline. We mainly focused on the effect of electrodes density, inverse method, and functional connectivity metric. In a recent study, we tested the effect of electrodes density and tested wMNE and eLORETA along with PLV and wPLI in the context of simulated epileptiform activity (Allouch et al., 2022). In this specific case, epileptiform signals had a sufficiently high signal-to-noise-ratio (SNR) to be distinguished from background noise (Wa 1983; Iwasaki et al., 2005), which might have facilitated the identification of underlying networks. However, the last two decades have witnessed an increased interest in resting-state studies where the significantly lower SNR could be in favor of different methods/parameter tuning in EEG source connectivity analysis. Therefore, in the present study, we used neural mass models to simulate resting-state brain activity (DMN and DAN) in the cortical space and derived the corresponding scalp EEG signals by solving the forward problem. Then, we reconstructed the corresponding cortical networks. This pipeline was repeated using different electrode densities (19, 32, 64, 128, and 256 channels) to solve the forward problem, three algorithms to solve the inverse prob-
<figure><img src="image_6.png" /><figcaption>Fig. 4. Bar plots of the percentage of edges located within the DMN among the 30 edges having the highest contribution values. eLORETA - exact low-resolution electromagnetic tomography. LCMV - linearly constrained minimum norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

lem (wMNE, eLORETA, LCMV), and four metrics to assess functional connectivity (PLV, AEC, PLI, leakage-corrected AEC (AEC*)).

## 4.1. Number of electrodes

To the best of our knowledge, the effect of the number of electrodes on the reconstruction of EEG-based resting-state cortical networks has never been studied before, especially in the presence of a ground truth enabling an objective comparison of the tested electrode configurations. Our results demonstrate a key role of electrodes density on the EEG source connectivity analysis: a more accurate reconstruction of cortical activity was achieved using high-density EEG (hd-EEG). This result is in line with previous simulations and empirical evidence (Song et al., 2015; Sohrabpour et al., 2015; Goran Lantz et al., 2003), as well as with the theoretical foundation (Srinivasan et al., 1998; Song et al., 2015). As established in (Srinivasan et al., 1998), an accurate characterization of spatial electrophysiological information requires a higher number of electrodes. A high inter-electrode distance (i.e., corresponding to a low number of EEG electrodes) can induce aliasing, and therefore high spatial frequency signals are misrepresented as low spatial frequency signals due to violation of the Nyquist criteria ($F_s > 2 \times F_{max}$) (Srinivasan et al., 1998; Song et al., 2015). In Song et al. (2015) results were found to be independent of the inverse method (minimum norm/standardized low-resolution brain electromagnetic tomography). On the other hand (Goran Lantz et al., 2003), showed that 63 electrodes were sufficient for decent source localization with EPIFOCUS (a linear inverse solution that optimally localizes single focal sources (Grave de Peralta Menendez et al. 2001; Lantz et al. (2001)), however, 100 electrodes were required when using wMNE. Consistently, (Goran Lantz et al. 2003) and (Sohrabpour et al., 2015) showed that a dramatic decrease in localization error was achieved when increasing the number of electrodes from 32 to 64 electrodes. Taken together, those studies and the present study support that the number of electrodes substantially affects network reconstruction. Based on our results, we recommend a minimum of 64 electrodes to be used in the course of EEG source reconstruction in the specific context of RSNs.

## 4.2. Inverse solution and connectivity measures

Regarding the source connectivity estimation, a plethora of methods offers the possibility to (1) reconstruct cortical activity dynamics, and (2) assess functional connectivity between reconstructed sources. However, there is no consensus over an optimal approach, nor on whether such a “best” technique exists. Testing different combinations of inverse methods and connectivity measures, we observed significant discrepancies in the results which was also observed in other comparative studies. In general, there is a lack of consistency across studies comparing inverse methods and connectivity measures. For example,

Anzolin et al. (2019) showed that LCMV had a better performance globally as compared to eLORETA. Similarly, in Mahjoory et al. (2017), a relatively strong difference was found between LCMV beamformer on one hand, and eLORETA/wMNE solutions on the other hand. In Hedrich et al. (2017) the coherent maximum entropy on the mean (cMEM) showed similar localization error to MNE, dynamic statistical parametric mapping (DSPM), sLORETA, but lower spatial spread and reduced crosstalk. In Bradley et al. (2016), the use of LORETA for source localization outperformed sLORETA and minimum norm least square (MNLS). Following an extensive comparison between six inverse methods, (Grova et al., 2006) recommended taking into account results from different methods when localizing actual interictal spikes. Results of source localization in Halder et al. (2019) did not identify a clear winner between LCMV, eLORETA, MNE, DISK. Tait et al. (2021) summarized the conditions where each method can be recommended, following a comparison of six inverse methods in resting-state MEG data. It is noteworthy to mention that the discrepancy in results could probably be related to the mathematical and physical constraints imposed by each inverse solution. For instance, wMNE searches for a solution with minimum power (Hämäläinen and Ilmoniemi 1994), while eLORETA tends to achieve exact zero error localization in the presence of measurement and structured biological noise (Pascual-Marqui 2007). On the other hand, LCMV beamformer assumes that only one dipole is active at a time, and tries to minimize the corresponding output energy (Van Veen et al. 1997).

In terms of functional connectivity metrics, (Colclough et al., 2016) assessed the consistency of different measures in experimental MEG resting-state data and recommended using the correlation between orthogonalized, band-limited, power envelopes (AEC). On the other hand, following extensive simulation studies, Wang et al. (2014) and Wendling et al. (2009) both concluded that there is no ideal “one-fit-all” method for all data types: it is rather suggested to evaluate which conditions are appropriate for each method. In Hassan et al. (2014) and Hassan et al. (2017) in the context of epileptic spikes, wMNE combined with PLV had better accuracy as compared to other algorithms. Let us note that, in Allouch et al. (2022), wMNE combined with wPLI performed better in the context of epileptiform activity simulations.

A key limitation of those numerous studies, and ours is no exception, is that different methods were tested in different contexts and using different data types, which further complicates the identification of clear and concise guidelines on the topic. Thus, we are aware that, even with our additional contribution, we are still far from a generalization of the results obtained in the specific context tested here (simulated signals, resting state activity, alpha rhythms, number/location of cortical sources, etc....). In fact, it is entirely possible (and even probable) that the inverse solution/connectivity measure combination is context-specific, and that no ideal method can account for all data types or all research questions. However, this raises the question of whether it is
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Bar plots of the percentage of edges located within the DMN among the 30 edges having the highest contribution values. eLORETA - *exact low-resolution electromagnetic tomography*. LCMV - *linearly constrained minimum* norm beamforming. wMNE - weighted minimum norm estimate. PLV - phase-locking value. AEC - amplitude envelope correlation. PLI - phase-lag index. AEC* - amplitude envelope correlation with source leakage correction.</figcaption></figure>

lem (wMNE, eLORETA, LCMV), and four metrics to assess functional connectivity (PLV, AEC, PLI, leakage-corrected AEC (AEC*)).

## 4.1. Number of electrodes

To the best of our knowledge, the effect of the number of electrodes on the reconstruction of EEG-based resting-state cortical networks has never been studied before, especially in the presence of a ground truth enabling an objective comparison of the tested electrode configurations. Our results demonstrate a key role of electrodes density on the EEG source connectivity analysis: a more accurate reconstruction of cortical activity was achieved using high-density EEG (hd-EEG). This result is in line with previous simulations and empirical evidence (Song et al., 2015; Sohrabpour et al., 2015; Goran Lantz et al., 2003), as well as with the theoretical foundation (Srinivasan et al., 1998; Song et al., 2015). As established in (Srinivasan et al., 1998), an accurate characterization of spatial electrophysiological information requires a higher number of electrodes. A high inter-electrode distance (i.e., corresponding to a low number of EEG electrodes) can induce aliasing, and therefore high spatial frequency signals are misrepresented as low spatial frequency signals due to violation of the Nyquist criteria ($F_s > 2 \times F_{max}$) (Srinivasan et al., 1998; Song et al., 2015). In Song et al. (2015) results were found to be independent of the inverse method (minimum norm/standardized low-resolution brain electromagnetic tomography). On the other hand (Goran Lantz et al., 2003), showed that 63 electrodes were sufficient for decent source localization with EPIFOCUS (a linear inverse solution that optimally localizes single focal sources (Grave de Peralta Menendez et al. 2001; Lantz et al. (2001)), however, 100 electrodes were required when using wMNE. Consistently, (Goran Lantz et al. 2003) and (Sohrabpour et al., 2015) showed that a dramatic decrease in localization error was achieved when increasing the number of electrodes from 32 to 64 electrodes. Taken together, those studies and the present study support that the number of electrodes substantially affects network reconstruction. Based on our results, we recommend a minimum of 64 electrodes to be used in the course of EEG source reconstruction in the specific context of RSNs.

## 4.2. Inverse solution and connectivity measures

Regarding the source connectivity estimation, a plethora of methods offers the possibility to (1) reconstruct cortical activity dynamics, and (2) assess functional connectivity between reconstructed sources. However, there is no consensus over an optimal approach, nor on whether such a “best” technique exists. Testing different combinations of inverse methods and connectivity measures, we observed significant discrepancies in the results which was also observed in other comparative studies. In general, there is a lack of consistency across studies comparing inverse methods and connectivity measures. For example,

Anzolin et al. (2019) showed that LCMV had a better performance globally as compared to eLORETA. Similarly, in Mahjoory et al. (2017), a relatively strong difference was found between LCMV beamformer on one hand, and eLORETA/wMNE solutions on the other hand. In Hedrich et al. (2017) the coherent maximum entropy on the mean (cMEM) showed similar localization error to MNE, dynamic statistical parametric mapping (DSPM), sLORETA, but lower spatial spread and reduced crosstalk. In Bradley et al. (2016), the use of LORETA for source localization outperformed sLORETA and minimum norm least square (MNLS). Following an extensive comparison between six inverse methods, (Grova et al., 2006) recommended taking into account results from different methods when localizing actual interictal spikes. Results of source localization in Halder et al. (2019) did not identify a clear winner between LCMV, eLORETA, MNE, DISK. Tait et al. (2021) summarized the conditions where each method can be recommended, following a comparison of six inverse methods in resting-state MEG data. It is noteworthy to mention that the discrepancy in results could probably be related to the mathematical and physical constraints imposed by each inverse solution. For instance, wMNE searches for a solution with minimum power (Hämäläinen and Ilmoniemi 1994), while eLORETA tends to achieve exact zero error localization in the presence of measurement and structured biological noise (Pascual-Marqui 2007). On the other hand, LCMV beamformer assumes that only one dipole is active at a time, and tries to minimize the corresponding output energy (Van Veen et al. 1997).

In terms of functional connectivity metrics, (Colclough et al., 2016) assessed the consistency of different measures in experimental MEG resting-state data and recommended using the correlation between orthogonalized, band-limited, power envelopes (AEC). On the other hand, following extensive simulation studies, Wang et al. (2014) and Wendling et al. (2009) both concluded that there is no ideal “one-fit-all” method for all data types: it is rather suggested to evaluate which conditions are appropriate for each method. In Hassan et al. (2014) and Hassan et al. (2017) in the context of epileptic spikes, wMNE combined with PLV had better accuracy as compared to other algorithms. Let us note that, in Allouch et al. (2022), wMNE combined with wPLI performed better in the context of epileptiform activity simulations.

A key limitation of those numerous studies, and ours is no exception, is that different methods were tested in different contexts and using different data types, which further complicates the identification of clear and concise guidelines on the topic. Thus, we are aware that, even with our additional contribution, we are still far from a generalization of the results obtained in the specific context tested here (simulated signals, resting state activity, alpha rhythms, number/location of cortical sources, etc....). In fact, it is entirely possible (and even probable) that the inverse solution/connectivity measure combination is context-specific, and that no ideal method can account for all data types or all research questions. However, this raises the question of whether it is
possible to determine the method that is the most adapted in each context, which requests extensive investigation far beyond simple comparative studies. Additionally, because inverse solutions and connectivity measures are interdependent, it is more challenging to determine the "winning" approach. One possibility is that the underlying assumptions of source reconstruction and connectivity algorithms may dictate which inverse solutions are suitable for a given set of connectivity measures, and which ones are not. In an ideal scenario, we would be able to establish that a connectivity measure X should only be applied with the inverse solution Y in a specific context Z. Meanwhile, cross-validation of the results using several methods/measures could be a workable compromise to avoid false positives brought on by such high analytical variability in the absence of a ground-truth and an ideal and clearly defined approach. In other cases, rather than focusing on the performance of a method to reconstruct actual networks, it may be more useful to focus on the method that produces the most useful results in terms of clinical applications. However, in this case, no inference could be drawn regarding the correspondence between the networks computed and the actual networks at the cortical level. In the absence of ground truth, differences between methods and pipelines could also be measured in terms of group-level repeatability, inter-, and intra-subject similarities, which illustrates the degree of variability in outcomes resulting from various analytical decisions. Finally, we believe that the recent approaches aiming at combining both steps (inversion and connectivity) in one step (see Colclough et al. (2018) for instance) is of great interest and can overcome some of the above-mentioned limitations.

A critical issue that affects connectivity estimation is source leakage. If the time lags introduced in our simulations are considered negligible by the algorithms, i.e., we do have zero-lag connections, then, in an ideal scenario, corrected AEC and PLI will omit those (true) connections, and PLV and AEC will perform better, which was not the case. If time lags are not negligible, all connectivity measures (PLV, PLI, AEC, corrected AEC) should perform equally, which was not the case. As a result, we believe that the source leakage introduced when solving the inverse problem has a significant impact on the performance of the connectivity metrics.

### 4.3. Methodological considerations

Taking together the conclusions from the vast majority of analysis and modeling choices, the identified networks might be also sensitive to other factors that were not investigated in this study. For instance, spatial resolution in the cortical space (i.e., number/size of parcellated ROIs) could affect the accuracy of reconstructed networks. Our simulations were all restricted to 66 ROIs (Desikan-Killiany atlas) due to the model design. Many graph-based studies have reported that different parcellation scales resulted in significant differences in network parameters (clustering coefficient, characteristic path length, local and global efficiency, degree distribution, etc....), while inferences about small-world and scale-free properties were maintained across different scales (Hayasaka and Laurienti, 2010; Wang et al., 2009; Zalesky et al., 2010; Fornito et al., 2010). In this study, we did not investigate how analytical variability related to spatial resolution in the source space (i.e., number/size of ROIs) could affect the results, especially with the effect of source leakage that could be more pronounced in very close regions. Therefore, it would be interesting to evaluate if the connectivity estimation was biased by the cortical spatial resolution or whether the optimal number of ROIs depends on the number of recording electrodes and vice versa, i.e., whether there exists an optimal ratio (number of electrodes/number of ROIs) that exhibits higher accuracy regarding source connectivity. Moreover, the effect of parcellation is not limited to the number of regions, but also involves the algorithms by which sources belonging to the same ROI are aggregated (averaging signals across all vertices, averaging their absolute values, power signals, keeping the first mode of PCA decomposition, choosing the maximum value among vertices at each time point...).

Here, we simulated two widely studied RSNs (DMN and DAN). The major difference between those two networks is the position of simulated sources (i.e., the DMN/DAN regions). While the DMN exhibits a less distributed architecture with small to moderate distance between regions (especially the right and left MPFC and PCC), DAN regions are widespread across the cortex with higher inter-region distances. Due to this difference, we were able to test whether the results are affected by the specific spatial location of sources. However, the similarity of results obtained with both networks confirms, to some extent, the absence of a bias caused by the position of sources. For simplicity, we restricted the number of DMN and DAN regions in our simulations to six and focused on the most consistent regions reported in the literature. However, other regions could be involved in the simulated DMN, such as the *precuneus*, *isthmus cingulate*, *rostral anterior cingulate*, and *lateral orbitofrontal cortex*. The DAN could also include the middle temporal gyrus and frontal regions, such as the *pars opercularis*, *pars orbitalis*, and *pars triangularis* of the inferior frontal cortex (Kabbara et al., 2017, 2021). An interesting prospect would be testing the consistency of results when more complex networks are involved.

In the simulations used for this study, alpha activity [8–12] Hz was assigned to RSN, while background activity was attributed to the remaining regions. We acknowledge that such contrast between those two types of activity is not completely realistic. Despite a potential significant interest, we were not able to simulate brain-like broadband signals covering all frequency bands, due to intrinsic model limitations. The first reason for choosing alpha activity is that the resting-state networks are usually identified in the alpha and beta frequency bands. The second reason is related to the model itself: while we were able to generate alpha-like activity, it was not possible to generate physiologically-plausible beta activity. Another limitation of our study is the “static” approach used for reconstructing brain networks (i.e., one network is estimated per epoch), which we could improve upon by investigating the brain networks’ dynamics. Importantly, understanding the dynamics of brain networks provides crucial insights regarding brain functions both in resting-state (Hipp et al., 2012; de Pasquale et al., 2018; De Pasquale et al., 2016; Allen et al., 2018, 2014; Baker et al., 2014; Kabbara et al., 2017, 2021) and task-related paradigms (Shine et al., 2016; O’Neill et al., 2017; Elton and Gao 2015; Braun et al., 2015; Krienen et al., 2014; Fong et al., 2019; Hassan et al., 2015).

As aforementioned, inverse solutions and connectivity algorithms require tuning several parameters, such as the regularization parameter in wMNE, eLORETA, LCMV, and the time window over which connectivity measures are computed. Here, we applied either default and usually used values proposed by the Fieldtrip and Brainstorm toolboxes, or values proposed in previous studies. To test the stability of reconstructed networks across different regularization values, we computed the Pearson correlation between the reference and reconstructed networks for different regularization values (PLV, 256 channels) and presented the results in Fig. S5 in the Supplementary Materials. A decrease in correlation values was observed when increasing the regularization of the inverse solution, emphasizing the importance of tuning the regularization parameter rather than using default values. Several methods are available to address this issue, falling into two general categories: those based on the estimation of the measurement noise, and those which are not (L-curve method, general-cross validation method, composite residual and smoothing operator, minimal product method, zero crossing) (Grech et al., 2008). An additional crucial factor that may introduce variability between results, and that emerges when dealing with experimental EEG data, is the pre-processing procedure that precedes network reconstruction. EEG signals are indeed usually contaminated by several artifacts and noise sources, that may each introduce undesired changes in the measurements and affect the signal of interest (Urigüen and Garcia-Zapirain, 2015), especially in clinical, pediatric and aging populations (Pedroni et al., 2019). To “clean” EEG signals, pre-processing algorithms usually impose different constraints for accepting and rejecting artificial epochs and propose different techniques to eliminate artifacts.
possible to determine the method that is the most adapted in each context, which requests extensive investigation far beyond simple comparative studies. Additionally, because inverse solutions and connectivity measures are interdependent, it is more challenging to determine the "winning" approach. One possibility is that the underlying assumptions of source reconstruction and connectivity algorithms may dictate which inverse solutions are suitable for a given set of connectivity measures, and which ones are not. In an ideal scenario, we would be able to establish that a connectivity measure X should only be applied with the inverse solution Y in a specific context Z. Meanwhile, cross-validation of the results using several methods/measures could be a workable compromise to avoid false positives brought on by such high analytical variability in the absence of a ground-truth and an ideal and clearly defined approach. In other cases, rather than focusing on the performance of a method to reconstruct actual networks, it may be more useful to focus on the method that produces the most useful results in terms of clinical applications. However, in this case, no inference could be drawn regarding the correspondence between the networks computed and the actual networks at the cortical level. In the absence of ground truth, differences between methods and pipelines could also be measured in terms of group-level repeatability, inter-, and intra-subject similarities, which illustrates the degree of variability in outcomes resulting from various analytical decisions. Finally, we believe that the recent approaches aiming at combining both steps (inversion and connectivity) in one step (see Colclough et al. (2018) for instance) is of great interest and can overcome some of the above-mentioned limitations.

A critical issue that affects connectivity estimation is source leakage. If the time lags introduced in our simulations are considered negligible by the algorithms, i.e., we do have zero-lag connections, then, in an ideal scenario, corrected AEC and PLI will omit those (true) connections, and PLV and AEC will perform better, which was not the case. If time lags are not negligible, all connectivity measures (PLV, PLI, AEC, corrected AEC) should perform equally, which was not the case. As a result, we believe that the source leakage introduced when solving the inverse problem has a significant impact on the performance of the connectivity metrics.

### 4.3. Methodological considerations

Taking together the conclusions from the vast majority of analysis and modeling choices, the identified networks might be also sensitive to other factors that were not investigated in this study. For instance, spatial resolution in the cortical space (i.e., number/size of parcellated ROIs) could affect the accuracy of reconstructed networks. Our simulations were all restricted to 66 ROIs (Desikan-Killiany atlas) due to the model design. Many graph-based studies have reported that different parcellation scales resulted in significant differences in network parameters (clustering coefficient, characteristic path length, local and global efficiency, degree distribution, etc...), while inferences about small-world and scale-free properties were maintained across different scales (Hayasaka and Laurienti, 2010; Wang et al., 2009; Zalesky et al., 2010; Fornito et al., 2010). In this study, we did not investigate how analytical variability related to spatial resolution in the source space (i.e., number/size of ROIs) could affect the results, especially with the effect of source leakage that could be more pronounced in very close regions. Therefore, it would be interesting to evaluate if the connectivity estimation was biased by the cortical spatial resolution or whether the optimal number of ROIs depends on the number of recording electrodes and vice versa, i.e., whether there exists an optimal ratio (number of electrodes/number of ROIs) that exhibits higher accuracy regarding source connectivity. Moreover, the effect of parcellation is not limited to the number of regions, but also involves the algorithms by which sources belonging to the same ROI are aggregated (averaging signals across all vertices, averaging their absolute values, power signals, keeping the first mode of PCA decomposition, choosing the maximum value among vertices at each time point...).

Here, we simulated two widely studied RSNs (DMN and DAN). The major difference between those two networks is the position of simulated sources (i.e., the DMN/DAN regions). While the DMN exhibits a less distributed architecture with small to moderate distance between regions (especially the right and left MPFC and PCC), DAN regions are widespread across the cortex with higher inter-region distances. Due to this difference, we were able to test whether the results are affected by the specific spatial location of sources. However, the similarity of results obtained with both networks confirms, to some extent, the absence of a bias caused by the position of sources. For simplicity, we restricted the number of DMN and DAN regions in our simulations to six and focused on the most consistent regions reported in the literature. However, other regions could be involved in the simulated DMN, such as the *precuneus*, *isthmus cingulate*, *rostral anterior cingulate*, and *lateral orbitofrontal cortex*. The DAN could also include the middle temporal gyrus and frontal regions, such as the *pars opercularis*, *pars orbitalis*, and *pars triangularis* of the inferior frontal cortex (Kabbara et al., 2017, 2021). An interesting prospect would be testing the consistency of results when more complex networks are involved.

In the simulations used for this study, alpha activity [8–12] Hz was assigned to RSN, while background activity was attributed to the remaining regions. We acknowledge that such contrast between those two types of activity is not completely realistic. Despite a potential significant interest, we were not able to simulate brain-like broadband signals covering all frequency bands, due to intrinsic model limitations. The first reason for choosing alpha activity is that the resting-state networks are usually identified in the alpha and beta frequency bands. The second reason is related to the model itself: while we were able to generate alpha-like activity, it was not possible to generate physiologically-plausible beta activity. Another limitation of our study is the “static” approach used for reconstructing brain networks (i.e., one network is estimated per epoch), which we could improve upon by investigating the brain networks’ dynamics. Importantly, understanding the dynamics of brain networks provides crucial insights regarding brain functions both in resting-state (Hipp et al., 2012; de Pasquale et al., 2018; De Pasquale et al., 2016; Allen et al., 2018, 2014; Baker et al., 2014; Kabbara et al., 2017, 2021) and task-related paradigms (Shine et al., 2016; O’Neill et al., 2017; Elton and Gao 2015; Braun et al., 2015; Krienen et al., 2014; Fong et al., 2019; Hassan et al., 2015).

As aforementioned, inverse solutions and connectivity algorithms require tuning several parameters, such as the regularization parameter in wMNE, eLORETA, LCMV, and the time window over which connectivity measures are computed. Here, we applied either default and usually used values proposed by the Fieldtrip and Brainstorm toolboxes, or values proposed in previous studies. To test the stability of reconstructed networks across different regularization values, we computed the Pearson correlation between the reference and reconstructed networks for different regularization values (PLV, 256 channels) and presented the results in Fig. S5 in the Supplementary Materials. A decrease in correlation values was observed when increasing the regularization of the inverse solution, emphasizing the importance of tuning the regularization parameter rather than using default values. Several methods are available to address this issue, falling into two general categories: those based on the estimation of the measurement noise, and those which are not (L-curve method, general-cross validation method, composite residual and smoothing operator, minimal product method, zero crossing) (Grech et al., 2008). An additional crucial factor that may introduce variability between results, and that emerges when dealing with experimental EEG data, is the pre-processing procedure that precedes network reconstruction. EEG signals are indeed usually contaminated by several artifacts and noise sources, that may each introduce undesired changes in the measurements and affect the signal of interest (Urigüen and Garcia-Zapirain, 2015), especially in clinical, pediatric and aging populations (Pedroni et al., 2019). To “clean” EEG signals, pre-processing algorithms usually impose different constraints for accepting and rejecting artificial epochs and propose different techniques to eliminate artifacts.
Hence, an important methodological consideration would be to inves-
tigate and systematically quantify the variability induced by different
data preprocessing techniques.

The choice of the most convenient metric for quantifying the results
is also a challenge. Here, to quantify the performance of the tested pa-
rameters, we used the (1) Pearson correlation computed between the
reference and reconstructed networks, (2) closeness accuracy reflect-
ing the distance between the nodes detected in the reconstructed and
reference networks, and (3) the percentage of edges having the highest
contributions to the correlation between networks and falling within the
simulated RSN. It should be noted that each metric measures different
aspects of data, and it is important in the context of comparative stud-
ies to address the relevant aspect of data to prevent biases. In our case,
the general trend in the results was preserved across all used metrics.
While both Pearson correlation and closeness accuracy showed signifi-
cant differences between tested parameters, the latter exhibited higher
variability between epochs (due to its mathematical formulation). Other
metrics such as networks-based metrics have also been developed to
evaluate the distance/similarity between networks, and a promising ap-
proach would be testing whether the differences between networks are
maintained across different aspects investigated by such metrics (global
analysis, edge-wise, node-wise, spectral graph analysis, edit distance,
kernel methods) (Mheich et al., 2020).

Importantly, correlation results presented in this article were com-
puted between non-thresholded weighted connectivity matrices. How-
ever, selecting a set of network edges for subsequent analysis while dis-
carding others remains a debatable subject in the network neuroscience
community. We usually face two main questions: (1) what method to
use to perform thresholding, and (2) the choice of the threshold value
when needed. Since we aimed to reduce further data manipulation and
subjective intervention in the analysis, no threshold was applied before
computing the Pearson correlation and edge contribution matrices in
our study. However, we also tested the effect of a proportional thresh-
old (1%) on correlation values and reported the results in Fig. S6 in
Supplementary materials. Although lower correlation values and higher
between-epochs variability were obtained, the general trend was similar
to that of non-thresholded matrices.

5. Conclusion

To summarize, we simulated EEG data corresponding to RSNs and
tested the effect of several key parameters (number of EEG electrodes
and inverse solution/connectivity measure combination) on network re-
construction accuracy. Different analytical choices led to high variabil-
ity in the resulting networks. In the context of RSN simulations, our re-
sults demonstrate, as expected, that an accurate cortical network recon-
struction requires a high number of EEG electrodes. Therefore, we rec-
ommend using hd-EEG (≥ 64 channels) to infer cortical dynamics from
recorded scalp EEG signals. In addition, we suggest a very careful choice
of the inverse solution/connectivity measure combination, since our re-
sults highlight a significant variability in the networks reconstructed us-
ing different inverse solutions and connectivity measures. Such method-
ological variability and absence of analysis standardization represent a
critical issue for neuroimaging studies that should be prioritized.

Data and codes availability

Data simulated for this study are available at 10.5281/zenodo.
6597385. Codes supporting the results of this study are avail-
able at https://github.com/sahar-allouch/var-RSNs. We used
Matlab (2018), Brainstorm toolbox (Tadel et al., 2011), Fieldtrip
toolbox (Oostenveld et al., 2011; http://fieldtriptoolbox.org), Open-
MEEG (Gramfort et al., 2010) implemented in Fieldtrip, R (R Core
Team, 2022) for statistical analysis, and Seaborn (Waskom, 2021) and
Matplotlib (Hunter, 2007) for visualization.

Declaration of Competing Interest

The authors declare that they have no competing interests.

Credit authorship contribution statement

**Sahar Allouch:** Conceptualization, Methodology, Software, Writing – original draft. **Aya Kabbara:** Conceptualization, Methodology, Supervision, Writing – review & editing. **Joan Duprez:** Conceptualization, Methodology, Supervision, Writing – review & editing. **Mohamad Khalil:** Funding acquisition. **Julien Modolo:** Conceptualization, Methodology, Supervision, Writing – review & editing. **Mahmoud Hassan:** Conceptualization, Methodology, Supervision, Writing – review & editing.

Data Availability

Data will be made available on request.

Acknowledgments

The authors would like to acknowledge the Lebanese National Coun-
cil for Scientific Research (CNRS-L), the Agence Universitaire de la Fran-
cophonie (AUF), and the Lebanese University for granting Ms. Allouch a
doctoral scholarship. This work was also supported by the Labex Comin-
labs project PKSTIM and by the Institute of Clinical Neuroscience of
Rennes (Projects EEGCog and EEGNET3). The authors would also like
to thank Campus France, Programme Hubert Curien CEDRE (PROJET N°
42257YA), and the Lebanese Association for Scientific Research (LASER)
for their support.

Supplementary materials

Supplementary material associated with this article can be found, in
the online version, at doi:10.1016/j.neuroimage.2023.120006.

References

Alavash, M., Doebler, P., Holling, H., Thiel, C.M., Gießing, C., 2015. Is functional integration of resting state brain networks an unspecific biomarker for working memory performance? *Neuroimage* 108 (March), 182–193.

Allen, E.A., Damaraju, E., Eichele, T., Wu, L., Calhoun, V.D., 2018. EEG signatures of dynamic functional network connectivity states. Brain Topogr. 31 (1), 101–116.

Allen, E.A., Damaraju, E., Plis, S.M., Erhardt, E.B., Eichele, T., Calhoun, V.D., 2014. Tracking whole-brain connectivity dynamics in the resting state. Cereb. Cortex 24, 663–676.

Allouch, S., Yochum, M., Kabbara, A., Duprez, J., Khalil, M., Wendling, F., Hassan, M., Modolo, J., 2022. Mean-field modeling of brain-scale dynamics for the evaluation of EEG source-space networks. Brain Topogr. 35 (1), 54–65.

Anzolin, A., Presti, P., Van De Steen, F., Astolfi, L., Haufe, S., Marinazzo, D., 2019. Quantifying the effect of demixing approaches on directed connectivity estimated between reconstructed EEG sources. Brain Topogr. doi:10.1007/s10548-019-00705-z.

Awan, F.G., Saleem, O., Kiran, A., 2019. Recent trends and advances in solving the inverse problem for EEG source localization. Inverse Prob. Sci. Eng. 27 (11), 1521–1536.

Baillet, S., Mosher, J.C., Leahy, R.M., 2001. Electromagnetic brain mapping. IEEE Signal Process. Mag. 14–30 no. November.

Baker, A.P., Brookes, M.J., Rezek, I.A., Smith, S.M., Behrens, T., Smith, P.J.P., Woolrich, M., 2014. Fast transient networks in spontaneous human brain activity. Elife 3, e01867 (March).

Barton, K., 2009. "MuMIN: multi-model inference." http://R-Forge.r-Project.Org/Projects/Mumin/. https://ci.nii.ac.jp/naid/20001420752/.

Bates, D., Mächler, M., Bolker, B.M., Walker, S.C., 2015. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67 (1). doi:10.18637/jss.v067.i01.

Bensaid, S., Modolo, J., Merlet, I., Wendling, F., Benquet, P., 2019. COALIA: a compu-
tational model of human EEG for consciousness research. Front. Syst. Neurosci. 13,
1–18.

Bradley, A., Yao, J., Dewald, J., Richter, C.P., 2016. Evaluation of electroencephalography source localization algorithms with multiple cortical sources. PLoS One 11 (1), 1–14.

Braun, U., Schäfer, A., Walter, H., Erk, S., Romanczuk-Seiferth, N., Haddad, L., Schweiger, J.I., et al., 2015. Dynamic reconfiguration of frontal brain networks during executive cognition in humans. Proc. Natl. Acad. Sci. U.S.A. 112 (37), 11678–11683.

Brookes, M.J., Hale, J.R., Zumer, J.M., Stevenson, C.M., Francis, S.T., Barnes, G.R., Owen, J.P., Morris, P.G., Nagarajan, S.S., 2011. NeuroImage measuring functional connectivity using MEG: methodology and comparison with FcMRI. Neuroimage 56 (3), 1082–1104.
Hence, an important methodological consideration would be to inves-
tigate and systematically quantify the variability induced by different
data preprocessing techniques.

The choice of the most convenient metric for quantifying the results
is also a challenge. Here, to quantify the performance of the tested pa-
rameters, we used the (1) Pearson correlation computed between the
reference and reconstructed networks, (2) closeness accuracy reflect-
ing the distance between the nodes detected in the reconstructed and
reference networks, and (3) the percentage of edges having the highest
contributions to the correlation between networks and falling within the
simulated RSN. It should be noted that each metric measures different
aspects of data, and it is important in the context of comparative stud-
ies to address the relevant aspect of data to prevent biases. In our case,
the general trend in the results was preserved across all used metrics.
While both Pearson correlation and closeness accuracy showed signifi-
cant differences between tested parameters, the latter exhibited higher
variability between epochs (due to its mathematical formulation). Other
metrics such as networks-based metrics have also been developed to
evaluate the distance/similarity between networks, and a promising ap-
proach would be testing whether the differences between networks are
maintained across different aspects investigated by such metrics (global
analysis, edge-wise, node-wise, spectral graph analysis, edit distance,
kernel methods) (Mheich et al., 2020).

Importantly, correlation results presented in this article were com-
puted between non-thresholded weighted connectivity matrices. How-
ever, selecting a set of network edges for subsequent analysis while dis-
carding others remains a debatable subject in the network neuroscience
community. We usually face two main questions: (1) what method to
use to perform thresholding, and (2) the choice of the threshold value
when needed. Since we aimed to reduce further data manipulation and
subjective intervention in the analysis, no threshold was applied before
computing the Pearson correlation and edge contribution matrices in
our study. However, we also tested the effect of a proportional thresh-
old (1%) on correlation values and reported the results in Fig. S6 in
Supplementary materials. Although lower correlation values and higher
between-epochs variability were obtained, the general trend was similar
to that of non-thresholded matrices.

5. Conclusion

To summarize, we simulated EEG data corresponding to RSNs and
tested the effect of several key parameters (number of EEG electrodes
and inverse solution/connectivity measure combination) on network re-
construction accuracy. Different analytical choices led to high variabil-
ity in the resulting networks. In the context of RSN simulations, our re-
sults demonstrate, as expected, that an accurate cortical network recon-
struction requires a high number of EEG electrodes. Therefore, we rec-
ommend using hd-EEG (≥ 64 channels) to infer cortical dynamics from
recorded scalp EEG signals. In addition, we suggest a very careful choice
of the inverse solution/connectivity measure combination, since our re-
sults highlight a significant variability in the networks reconstructed us-
ing different inverse solutions and connectivity measures. Such method-
ological variability and absence of analysis standardization represent a
critical issue for neuroimaging studies that should be prioritized.

Data and codes availability

Data simulated for this study are available at 10.5281/zenodo.
6597385. Codes supporting the results of this study are avail-
able at https://github.com/sahar-allouch/var-RSNs. We used
Matlab (2018), Brainstorm toolbox (Tadel et al., 2011), Fieldtrip
toolbox (Oostenveld et al., 2011; http://fieldtriptoolbox.org), Open-
MEEG (Gramfort et al., 2010) implemented in Fieldtrip, R (R Core
Team, 2022) for statistical analysis, and Seaborn (Waskom, 2021) and
Matplotlib (Hunter, 2007) for visualization.

Declaration of Competing Interest

The authors declare that they have no competing interests.

Credit authorship contribution statement

**Sahar Allouch:** Conceptualization, Methodology, Software, Writing – original draft. **Aya Kabbara:** Conceptualization, Methodology, Supervision, Writing – review & editing. **Joan Duprez:** Conceptualization, Methodology, Supervision, Writing – review & editing. **Mohamad Khalil:** Funding acquisition. **Julien Modolo:** Conceptualization, Methodology, Supervision, Writing – review & editing. **Mahmoud Hassan:** Conceptualization, Methodology, Supervision, Writing – review & editing.

Data Availability

Data will be made available on request.

Acknowledgments

The authors would like to acknowledge the Lebanese National Coun-
cil for Scientific Research (CNRS-L), the Agence Universitaire de la Fran-
cophonie (AUF), and the Lebanese University for granting Ms. Allouch a
doctoral scholarship. This work was also supported by the Labex Comin-
labs project PKSTIM and by the Institute of Clinical Neuroscience of
Rennes (Projects EEGCog and EEGNET3). The authors would also like
to thank Campus France, Programme Hubert Curien CEDRE (PROJET N°
42257YA), and the Lebanese Association for Scientific Research (LASER)
for their support.

Supplementary materials

Supplementary material associated with this article can be found, in
the online version, at doi:10.1016/j.neuroimage.2023.120006.

References

Alavash, M., Doebler, P., Holling, H., Thiel, C.M., Gießing, C., 2015. Is functional integration of resting state brain networks an unspecific biomarker for working memory performance? *Neuroimage* 108 (March), 182–193.

Allen, E.A., Damaraju, E., Eichele, T., Wu, L., Calhoun, V.D., 2018. EEG signatures of dynamic functional network connectivity states. Brain Topogr. 31 (1), 101–116.

Allen, E.A., Damaraju, E., Plis, S.M., Erhardt, E.B., Eichele, T., Calhoun, V.D., 2014. Tracking whole-brain connectivity dynamics in the resting state. Cereb. Cortex 24, 663–676.

Allouch, S., Yochum, M., Kabbara, A., Duprez, J., Khalil, M., Wendling, F., Hassan, M., Modolo, J., 2022. Mean-field modeling of brain-scale dynamics for the evaluation of EEG source-space networks. Brain Topogr. 35 (1), 54–65.

Anzolin, A., Presti, P., Van De Steen, F., Astolfi, L., Haufe, S., Marinazzo, D., 2019. Quantifying the effect of demixing approaches on directed connectivity estimated between reconstructed EEG sources. Brain Topogr. doi:10.1007/s10548-019-00705-z.

Awan, F.G., Saleem, O., Kiran, A., 2019. Recent trends and advances in solving the inverse problem for EEG source localization. Inverse Prob. Sci. Eng. 27 (11), 1521–1536.

Baillet, S., Mosher, J.C., Leahy, R.M., 2001. Electromagnetic brain mapping. IEEE Signal Process. Mag. 14–30 no. November.

Baker, A.P., Brookes, M.J., Rezek, I.A., Smith, S.M., Behrens, T., Smith, P.J.P., Woolrich, M., 2014. Fast transient networks in spontaneous human brain activity. Elife 3, e01867 (March).

Barton, K., 2009. "MuMIN: multi-model inference." http://R-Forge.r-Project.Org/Projects/Mumin/. https://ci.nii.ac.jp/naid/20001420752/.

Bates, D., Mächler, M., Bolker, B.M., Walker, S.C., 2015. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67 (1). doi:10.18637/jss.v067.i01.

Bensaid, S., Modolo, J., Merlet, I., Wendling, F., Benquet, P., 2019. COALIA: a compu-
tational model of human EEG for consciousness research. Front. Syst. Neurosci. 13,
1–18.

Bradley, A., Yao, J., Dewald, J., Richter, C.P., 2016. Evaluation of electroencephalography source localization algorithms with multiple cortical sources. PLoS One 11 (1), 1–14.

Braun, U., Schäfer, A., Walter, H., Erk, S., Romanczuk-Seiferth, N., Haddad, L., Schweiger, J.I., et al., 2015. Dynamic reconfiguration of frontal brain networks during executive cognition in humans. Proc. Natl. Acad. Sci. U.S.A. 112 (37), 11678–11683.

Brookes, M.J., Hale, J.R., Zumer, J.M., Stevenson, C.M., Francis, S.T., Barnes, G.R., Owen, J.P., Morris, P.G., Nagarajan, S.S., 2011. NeuroImage measuring functional connectivity using MEG: methodology and comparison with FcMRI. Neuroimage 56 (3), 1082–1104.
Cao, J., Zhao, Y., Shan, X., Wei, H.-L., Guo, Y., Chen, L., Erkoyuncu, J.A., Sarrigian-
nis, P.G., 2022. Brain functional and effective connectivity based on electroen-
cephalography recordings: a review. Hum. Brain Mapp. 43 (2), 860-879.

Colclough, G.L., Brookes, M.J., Smith, S.M., Woolrich, M.W., 2015. A symmetric multivari-
ate leakage correction for MEG connectomes. Neuroimage 117 (August), 439-448.

Colclough, G.L., Woolrich, M.W., Tewarie, P.K., Brookes, M.J., Quinn, A.J., Smith, S.M., 2016. How reliable are MEG resting-state connectivity metrics? Neuroimage 138, 284-293.

Colclough, G.L., Woolrich, M.W., Harrison, S.J., López, P.A.R., Valdes-Sosa, P.A., Smith, S.M., 2018. Multi-subject hierarchical inverse covariance modelling improves estimation of functional brain networks. Neuroimage 178 (September), 370-384.

Damoiseaux, J.S., Rombouts, S.A.R.B., Barkhof, F., Scheltens, P., Stam, C.J., Smith, S.M., Beckmann, C.F., 2006. Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U.S.A. 13 (37), 13848-13853.

De Pasquale, F., Penna, S.D., Sporns, O., Romani, G.L., Corbetta, M., 2016. A dynamic core network and global efficiency in the resting human brain. Cereb. Cortex 26 (10), 4015-4033.

Desikan, R.S., Ségonne, F., Fischl, B., Quinn, B.T., Dickerson, B.C., Blacker, D., Buck-
ner, R.L., et al., 2006. An automated labeling system for subdividing the human
cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage 31,
968-980.

Elton, A., Gao, W., 2015. Task-related modulation of functional connectivity variability and its behavioral correlations. Hum. Brain Mapp. 36 (8), 3260-3272.

Finn, E.S., Shen, X., Scheinost, D., Rosenberg, M.D., Huang, J., Chun, M.M., Pa-
pademetris, X., Constable, R.T., 2015. Functional connectome fingerprinting: iden-
tifying individuals using patterns of brain connectivity. Nat. Neurosci. 18 (11),
1664-1671.

Fong, A.H.C., Yoo, K., Rosenberg, M.D., Zhang, S., Li, C.S.R., Scheinost, D., Todd Consta-
ble, R., Chun, M.M., 2019. Dynamic functional connectivity during task performance
and rest predicts individual differences in attention across studies. Neuroimage 188
(March), 14-25.

Fornito, A., Zalesky, A., Bullmore, E.T., 2010. Network scaling effects in graph analytic studies of human resting-state fMRI data. Front. Syst. Neurosci. 4, 1-16.

Fox, J., Weisberg, S., 2019. An R Companion to Applied Regression (Third). Sage
doi:10.1177/0049124105277200.

Fox, M.D., Raichle, M.E., 2007. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8 (9), 700-711.

Friston, K.J., 2011. Functional and effective connectivity: a review. Brain Connect 1 (1),
13-36.

Fuchs, M., Wagner, M., Köhler, T., Wischmann, H.A., 1999. Linear and nonlinear current density reconstructions. J. Clin. Neurophysiol. Off. Publ. Am. Electroencephalogr. Soc. 16 (3), 267-295.

Gramfort, A., Papadopoulo, T., Olivi, E., Clerc, M., 2010. OpenMEEG: open-source software for quasistatic bioelectromagnetics. Biomed. Eng. Online 9 (45), doi:10.1186/1475-925X-8-1.

Gratton, C., Koller, J.M., Shannon, W., Greene, D.J., Maiti, B., Snyder, A.Z., Petersen, S.E., Perlmutter, J.S., Campbell, M.C., 2018. Emergent functional network effects in Parkinson disease. Cereb. Cortex 29 (6), 2509-2523.

Grave de Peralta Menendez, R., Gonzalez Andino, S., Lantz, G., Michel, C.M., Landis, T., 2001. Noninvasive localization of electromagnetic epileptic activity. I. Method descriptions and simulations. Brain Topogr. 14 (2), 131-137.

Grech, R., Cassar, T., Muscat, J., Camilleri, K.P., Fabri, S.G., Zervakis, M., Xanthopoulos, S., Sakkalis, V., Vanrumste, B., 2008. Review on solving the inverse problem in EEG source analysis. J. Neuroeng. Rehabil. 5 (25). doi:10.1186/1743-0003-5-25.

Greicius, M.D., Krasnow, B., Reiss, A.L., Menon, V., 2003. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl. Acad. Sci. U.S.A. 100 (1), 253-258.

Grova, C., Daunizeau, J., Lina, J.M., Bénar, C.G., Benali, H., Gotman, J., 2006. Evaluation of EEG localization methods using realistic simulations of interictal spikes. Neuroimage 29, 734-753.

Halder, T., Talwar, S., Kumar Jaiswal, A., Banerjee, A., 2019. Quantitative evaluation in estimating sources underlying brain oscillations using current source density methods and beamformer approaches. eNeuro 6 (4), 1-14.

Hämäläinen, M.S., Ilmoimiemi, R.J., 1994. Interpreting magnetic fields of the brain: minimum-norm estimates. Med. Biol. Eng. Comput. 32, 35-42.

Hassan, M., Benquet, P., Biraben, A., Berrou, C., Dufor, O., Wendling, F., 2015. Dynamic reorganization of functional brain networks during picture naming. Cortex 73, 276-288.

Hassan, M., Chaton, L., Benquet, P., Delval, A., Leroy, C., Plomhause, L., Moonen, A.J.H., et al., 2017a. Functional connectivity disruptions correlate with cognitive phenotypes in Parkinson's disease. Neuroimage Clin. 14 (March), 591-601.

Hassan, M., Dufor, O., Merlet, I., Berrou, C., Wendling, F., 2014. EEG source connectivity analysis: from dense array recordings to brain networks. PLoS One 9 (8), e105041.

Hassan, M., Merlet, I., Mheich, A., Kabbara, A., Biraben, A., Nica, A., Wendling, F., 2017b. Identification of interictal epileptic networks from dense-EEG. Brain Topogr. 30 (1), 60-76.

Hassan, M., Wendling, F., 2018. Electroencephalography source connectivity: aiming for high resolution of brain networks in time and space. IEEE Signal Process. Mag. 35 (3), 81-96.

Hayasaka, S., Laurienti, P.J., 2010. Comparison of characteristics between region-and vo-
el-based network analyses in resting-state fMRI data. Neuroimage 50, 499-508.

Hedrich, T., Pellegrino, G., Kobayashi, E., Lina, J.M., Grova, C., 2017. Comparison of the spatial resolution of source imaging techniques in high-density EEG and MEG. Neuroimage 157 (August), 531-544.

Hipp, J.F., Hawellek, D.J., Corbetta, M., Siegel, M., Engel, AK., 2012. Large-scale corti-
cal correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15 (6),
884-890.

Hothorn, T., Bretz, F., Westfall, P., 2008. Simultaneous inference in general parametric models. *Biom. J.* 50 (3), 346-363.

Hunter, 2007. Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9 (May), 90-95.

Iwasaki, M., Pestana, E., Burgess, RC., Lüders, HO., Shamoto, H., Nakasato, N., 2005. Detection of epileptiform activity by human interpreters: blinded comparison between electroencephalography and magnetoencephalography. Epilepsia 46 (1), 59-68.

Jimenez, AM., Riedel, F., Lee, J., Reavis, EA., Green, MF., 2019. Linking resting-state networks and social cognition in schizophrenia and bipolar disorder. Hum. Brain Mapp. 40 (16), 4703-4715.

Jockowitz, C., Caspers, S., Lux, S., Eickhoff, SB., Jütten, K., Lenzen, S., Moebus, S., et al., 2017. Influence of age and cognitive performance on resting-state brain networks of older adults in a population-based cohort. Cortex 89 (April), 28-44.

Kabbara, A., Eid, H., El Falou, W., Khalil, M., Wendling, F., Hassan, M., 2018. Reduced integration and improved segregation of functional brain networks in Alzheimer's disease. J. Neural Eng. 15 (2), 026023.

Kabbara, A., El Falou, W., Khalil, M., Wendling, F., Hassan, M., 2017. The dynamic func-
tional core network of the human brain at rest. Sci. Rep. 7 (1), 2936.

Kabbara, A., Paban, V., Hassan, M., 2021. The dynamic modular fingerprints of the human brain at rest. Neuroimage 227 (February), 117674.

Keller, J.B., Hedden, T., Thompson, T.W., Anteraper, S.A., Gabrieli, J.D.E., Whitfield-Gabrieli, S., 2015. Resting-state anticorrelations between medial and lateral prefrontal cortex: association with working memory, aging, and individual differences. Cortex 64 (March), 271-280.

Krienen, F.M., Thomas Yeo, B.T., Buckner, R.L., 2014. Reconfigurable task-dependent functional coupling modes cluster around a core functional architecture. Philos. Trans. R. Soc. Lond. B Biol. Sci. 369 (1653). doi:10.1098/rstb.2013.0526.

Lachaux, J.-P., Rodriguez, E., Le Van Quyen, M., Lutz, A., Martinerie, J., Varela, F.J., 2000. Studying single-trials of phase synchronous activity in the brain. Int. J. Bifurcat. Chaos 10 (10), 2429-2439.

Lantz, G., Grave de Peralta Menendez, R., Gonzalez Andino, S., M.ichel, C.M., 2001. Nonin-
vasive localization of electromagnetic epileptic activity. II. Demonstration of sublobar
accuracy in patients with simultaneous surface and depth recordings. Brain Topogr.
14 (2), 139-147.

Lantz, G., Grave de Peralta, R., Spinelli, L., Seeck, M., M.ichel, C.M., 2003. Epileptic source localization with high density EEG: how many electrodes are needed? Clin. Neurophysiol. Off. J. Int. Fed. Clin. Neurophysiol. 114 (1), 63-69.

Lin, F.H., Witzel, T., Ahlfors, S.P., Stufflebeam, S.M., Belliveau, J.W., Hämäläinen, M.S., 2006. Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates. Neuroimage 31, 160-171.

Mahjooy, K., Nikulin, V.V., Botrel, L., Linkenkaer-Hansen, K., Fato, M.M., Haufe, S., 2017. Consistency of EEG source localization and connectivity estimates. Neuroimage 152 (May), 590-601.

Matlab, 2018. "9.5.0.944444 (R2018b)". Natick, Massachusetts: The MathWorks Inc.

Mheich, A., Wendling, F., Hassan, M., 2020. Brain network similarity: methods and applications. Netw. Neurosci. 4 (3), 507-527.

O'Neill, G.C., Tewarie, P.K., Colclough, G.L., Gascoyne, L.E., Hunt, B.A.E., Morris, P.G., Woolrich, M.W., Brookes, M.J., 2017. Measurement of dynamic task related functional networks using MEG. Neuroimage 146, 667-678.

Oostenveld, R., Fries, P., Maris, E., Schoffelen, J.M., 2011. FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 2011. doi:10.1155/2011/156869.

Pascual-Marqui, R.D., 2007. "Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization." http://arxiv.org/abs/0710.3341.

Pasquale, F., Corbetta, M., Betti, V., Della Penna, S., 2018. Cortical cores in network dynamics. Neuroimage 180 (Pt B), 370-382.

Pasquale, F., Penna, S.D., Snyder, A.Z., Marzetti, L., Pizzella, V., Romani, G.L., Corbetta, M., 2012. A cortical core for dynamic integration of functional networks in the resting human brain. Neuron 74 (4), 753-764.

Pedroni, A., Bahreini, A., Langer, N., 2019. Automagic: standardized preprocessing of big EEG data. Neuroimage 200 (October), 460-473.

Pereda, E., Quiroga, R.Q., Bhattacharya, J., 2005. Nonlinear multivariate analysis of neu-
rophysiological signals. Prog. Neurobiol. 77, 1-37.

R Core Team. 2022. *R: A Language and Environment for Statistical Computing*. Vienna, Austria. https://www.R-project.org/.

Raichle, M.E., MacLeod, A.M., Snyder, A.Z., Powers, W.J., Gusnard, D.A., Shulman, G.L., 2001. A default mode of brain function. Proc. Natl. Acad. Sci. U.S.A. 98 (2), 676-682.

Raichle, M.E., Snyder, A.Z., 2007. A default mode of brain function: a brief history of an evolving idea. Neuroimage 37 (4), 1083-1090.

Sheline, Y.I., Raichle, M.E., 2013. Resting state functional connectivity in preclinical Alzheimer's disease. Biol. Psychiatry 74 (5), 340-347.

Shen, X., Cox, S.R., Adams, M.J., Howard, D.M., Lawrie, S.M., Ritchie, S.J., Bastin, M.E., Deary, I.J., McIntosh, A.M., Whalley, H.C., 2018. Resting-state connectivity and its association with cognitive performance, educational attainment, and household income in the UK Biobank. Biol. Psychiatry 3 (10), 878-886 *Cognitive Neuroscience and Neuroimaging*.

Shine, J.M., Bissett, P.G., Bell, P.T., Koyejo, O., Balsters, J.H., Gorgolewski, K.J., Moodie, C.A., Poldrack, R.A., 2016. The dynamics of functional brain networks: integrated network states during cognitive task performance. Neuron 92 (2), 544-554.

Shirer, W.R., Ryali, S., Rykhlevskiaia, E., Menon, V., G.reicius, M.D., 2012. Decoding sub-
ject-driven cognitive states with whole-brain connectivity patterns. Cereb. Cortex 22
(1), 158-165.
Cao, J., Zhao, Y., Shan, X., Wei, H.-L., Guo, Y., Chen, L., Erkoyuncu, J.A., Sarrigian-
nis, P.G., 2022. Brain functional and effective connectivity based on electroen-
cephalography recordings: a review. Hum. Brain Mapp. 43 (2), 860-879.

Colclough, G.L., Brookes, M.J., Smith, S.M., Woolrich, M.W., 2015. A symmetric multivari-
ate leakage correction for MEG connectomes. Neuroimage 117 (August), 439-448.

Colclough, G.L., Woolrich, M.W., Tewarie, P.K., Brookes, M.J., Quinn, A.J., Smith, S.M., 2016. How reliable are MEG resting-state connectivity metrics? Neuroimage 138, 284-293.

Colclough, G.L., Woolrich, M.W., Harrison, S.J., López, P.A.R., Valdes-Sosa, P.A., Smith, S.M., 2018. Multi-subject hierarchical inverse covariance modelling improves estimation of functional brain networks. Neuroimage 178 (September), 370-384.

Damoiseaux, J.S., Rombouts, S.A.R.B., Barkhof, F., Scheltens, P., Stam, C.J., Smith, S.M., Beckmann, C.F., 2006. Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U.S.A. 13 (37), 13848-13853.

De Pasquale, F., Penna, S.D., Sporns, O., Romani, G.L., Corbetta, M., 2016. A dynamic core network and global efficiency in the resting human brain. Cereb. Cortex 26 (10), 4015-4033.

Desikan, R.S., Ségonne, F., Fischl, B., Quinn, B.T., Dickerson, B.C., Blacker, D., Buck-
ner, R.L., et al., 2006. An automated labeling system for subdividing the human
cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage 31,
968-980.

Elton, A., Gao, W., 2015. Task-related modulation of functional connectivity variability and its behavioral correlations. Hum. Brain Mapp. 36 (8), 3260-3272.

Finn, E.S., Shen, X., Scheinost, D., Rosenberg, M.D., Huang, J., Chun, M.M., Pa-
pademetris, X., Constable, R.T., 2015. Functional connectome fingerprinting: iden-
tifying individuals using patterns of brain connectivity. Nat. Neurosci. 18 (11),
1664-1671.

Fong, A.H.C., Yoo, K., Rosenberg, M.D., Zhang, S., Li, C.S.R., Scheinost, D., Todd Consta-
ble, R., Chun, M.M., 2019. Dynamic functional connectivity during task performance
and rest predicts individual differences in attention across studies. Neuroimage 188
(March), 14-25.

Fornito, A., Zalesky, A., Bullmore, E.T., 2010. Network scaling effects in graph analytic studies of human resting-state fMRI data. Front. Syst. Neurosci. 4, 1-16.

Fox, J., Weisberg, S., 2019. An R Companion to Applied Regression (Third). Sage
doi:10.1177/0049124105277200.

Fox, M.D., Raichle, M.E., 2007. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8 (9), 700-711.

Friston, K.J., 2011. Functional and effective connectivity: a review. Brain Connect 1 (1),
13-36.

Fuchs, M., Wagner, M., Köhler, T., Wischmann, H.A., 1999. Linear and nonlinear current density reconstructions. J. Clin. Neurophysiol. Off. Publ. Am. Electroencephalogr. Soc. 16 (3), 267-295.

Gramfort, A., Papadopoulo, T., Olivi, E., Clerc, M., 2010. OpenMEEG: open-source software for quasistatic bioelectromagnetics. Biomed. Eng. Online 9 (45), doi:10.1186/1475-925X-8-1.

Gratton, C., Koller, J.M., Shannon, W., Greene, D.J., Maiti, B., Snyder, A.Z., Petersen, S.E., Perlmutter, J.S., Campbell, M.C., 2018. Emergent functional network effects in Parkinson disease. Cereb. Cortex 29 (6), 2509-2523.

Grave de Peralta Menendez, R., Gonzalez Andino, S., Lantz, G., Michel, C.M., Landis, T., 2001. Noninvasive localization of electromagnetic epileptic activity. I. Method descriptions and simulations. Brain Topogr. 14 (2), 131-137.

Grech, R., Cassar, T., Muscat, J., Camilleri, K.P., Fabri, S.G., Zervakis, M., Xanthopoulos, S., Sakkalis, V., Vanrumste, B., 2008. Review on solving the inverse problem in EEG source analysis. J. Neuroeng. Rehabil. 5 (25). doi:10.1186/1743-0003-5-25.

Greicius, M.D., Krasnow, B., Reiss, A.L., Menon, V., 2003. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl. Acad. Sci. U.S.A. 100 (1), 253-258.

Grova, C., Daunizeau, J., Lina, J.M., Bénar, C.G., Benali, H., Gotman, J., 2006. Evaluation of EEG localization methods using realistic simulations of interictal spikes. Neuroimage 29, 734-753.

Halder, T., Talwar, S., Kumar Jaiswal, A., Banerjee, A., 2019. Quantitative evaluation in estimating sources underlying brain oscillations using current source density methods and beamformer approaches. eNeuro 6 (4), 1-14.

Hämäläinen, M.S., Ilmoimiemi, R.J., 1994. Interpreting magnetic fields of the brain: minimum-norm estimates. Med. Biol. Eng. Comput. 32, 35-42.

Hassan, M., Benquet, P., Biraben, A., Berrou, C., Dufor, O., Wendling, F., 2015. Dynamic reorganization of functional brain networks during picture naming. Cortex 73, 276-288.

Hassan, M., Chaton, L., Benquet, P., Delval, A., Leroy, C., Plomhause, L., Moonen, A.J.H., et al., 2017a. Functional connectivity disruptions correlate with cognitive phenotypes in Parkinson's disease. Neuroimage Clin. 14 (March), 591-601.

Hassan, M., Dufor, O., Merlet, I., Berrou, C., Wendling, F., 2014. EEG source connectivity analysis: from dense array recordings to brain networks. PLoS One 9 (8), e105041.

Hassan, M., Merlet, I., Mheich, A., Kabbara, A., Biraben, A., Nica, A., Wendling, F., 2017b. Identification of interictal epileptic networks from dense-EEG. Brain Topogr. 30 (1), 60-76.

Hassan, M., Wendling, F., 2018. Electroencephalography source connectivity: aiming for high resolution of brain networks in time and space. IEEE Signal Process. Mag. 35 (3), 81-96.

Hayasaka, S., Laurienti, P.J., 2010. Comparison of characteristics between region-and vo-
el-based network analyses in resting-state fMRI data. Neuroimage 50, 499-508.

Hedrich, T., Pellegrino, G., Kobayashi, E., Lina, J.M., Grova, C., 2017. Comparison of the spatial resolution of source imaging techniques in high-density EEG and MEG. Neuroimage 157 (August), 531-544.

Hipp, J.F., Hawellek, D.J., Corbetta, M., Siegel, M., Engel, AK., 2012. Large-scale corti-
cal correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15 (6),
884-890.

Hothorn, T., Bretz, F., Westfall, P., 2008. Simultaneous inference in general parametric models. *Biom. J.* 50 (3), 346-363.

Hunter, 2007. Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9 (May), 90-95.

Iwasaki, M., Pestana, E., Burgess, RC., Lüders, HO., Shamoto, H., Nakasato, N., 2005. Detection of epileptiform activity by human interpreters: blinded comparison between electroencephalography and magnetoencephalography. Epilepsia 46 (1), 59-68.

Jimenez, AM., Riedel, F., Lee, J., Reavis, EA., Green, MF., 2019. Linking resting-state networks and social cognition in schizophrenia and bipolar disorder. Hum. Brain Mapp. 40 (16), 4703-4715.

Jockowitz, C., Caspers, S., Lux, S., Eickhoff, SB., Jütten, K., Lenzen, S., Moebus, S., et al., 2017. Influence of age and cognitive performance on resting-state brain networks of older adults in a population-based cohort. Cortex 89 (April), 28-44.

Kabbara, A., Eid, H., El Falou, W., Khalil, M., Wendling, F., Hassan, M., 2018. Reduced integration and improved segregation of functional brain networks in Alzheimer's disease. J. Neural Eng. 15 (2), 026023.

Kabbara, A., El Falou, W., Khalil, M., Wendling, F., Hassan, M., 2017. The dynamic func-
tional core network of the human brain at rest. Sci. Rep. 7 (1), 2936.

Kabbara, A., Paban, V., Hassan, M., 2021. The dynamic modular fingerprints of the human brain at rest. Neuroimage 227 (February), 117674.

Keller, J.B., Hedden, T., Thompson, T.W., Anteraper, S.A., Gabrieli, J.D.E., Whitfield-Gabrieli, S., 2015. Resting-state anticorrelations between medial and lateral prefrontal cortex: association with working memory, aging, and individual differences. Cortex 64 (March), 271-280.

Krienen, F.M., Thomas Yeo, B.T., Buckner, R.L., 2014. Reconfigurable task-dependent functional coupling modes cluster around a core functional architecture. Philos. Trans. R. Soc. Lond. B Biol. Sci. 369 (1653). doi:10.1098/rstb.2013.0526.

Lachaux, J.-P., Rodriguez, E., Le Van Quyen, M., Lutz, A., Martinerie, J., Varela, F.J., 2000. Studying single-trials of phase synchronous activity in the brain. Int. J. Bifurcat. Chaos 10 (10), 2429-2439.

Lantz, G., Grave de Peralta Menendez, R., Gonzalez Andino, S., M.ichel, C.M., 2001. Nonin-
vasive localization of electromagnetic epileptic activity. II. Demonstration of sublobar
accuracy in patients with simultaneous surface and depth recordings. Brain Topogr.
14 (2), 139-147.

Lantz, G., Grave de Peralta, R., Spinelli, L., Seeck, M., Michel, C.M., 2003. Epileptic source localization with high density EEG: how many electrodes are needed? Clin. Neurophysiol. Off. J. Int. Fed. Clin. Neurophysiol. 114 (1), 63-69.

Lin, F.H., Witzel, T., Ahlfors, S.P., Stufflebeam, S.M., Belliveau, J.W., Hämäläinen, M.S., 2006. Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates. Neuroimage 31, 160-171.

Mahjooy, K., Nikulin, V.V., Botrel, L., Linkenkaer-Hansen, K., Fato, M.M., Haufe, S., 2017. Consistency of EEG source localization and connectivity estimates. Neuroimage 152 (May), 590-601.

Matlab, 2018. "9.5.0.944444 (R2018b)". Natick, Massachusetts: The MathWorks Inc.

Mheich, A., Wendling, F., Hassan, M., 2020. Brain network similarity: methods and applications. Netw. Neurosci. 4 (3), 507-527.

O'Neill, G.C., Tewarie, P.K., Colclough, G.L., Gascoyne, L.E., Hunt, B.A.E., Morris, P.G., Woolrich, M.W., Brookes, M.J., 2017. Measurement of dynamic task related functional networks using MEG. Neuroimage 146, 667-678.

Oostenveld, R., Fries, P., Maris, E., Schoffelen, J.M., 2011. FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput. Intell. Neurosci. 2011. doi:10.1155/2011/156869.

Pascual-Marqui, R.D., 2007. "Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization." http://arxiv.org/abs/0710.3341.

Pasquale, F., Corbetta, M., Betti, V., Della Penna, S., 2018. Cortical cores in network dynamics. Neuroimage 180 (Pt B), 370-382.

Pasquale, F., Penna, S.D., Snyder, A.Z., Marzetti, L., Pizzella, V., Romani, G.L., Corbetta, M., 2012. A cortical core for dynamic integration of functional networks in the resting human brain. Neuron 74 (4), 753-764.

Pedroni, A., Bahreini, A., Langer, N., 2019. Automagic: standardized preprocessing of big
EEG data. Neuroimage 200 (October), 460-473.

Pereda, E., Quiroga, R.Q., Bhattacharya, J., 2005. Nonlinear multivariate analysis of neu-
rophysiological signals. Prog. Neurobiol. 77, 1-37.

R Core Team. 2022. *R: A Language and Environment for Statistical Computing*. Vienna, Austria. https://www.R-project.org/.

Raichle, M.E., MacLeod, A.M., Snyder, A.Z., Powers, W.J., Gusnard, D.A., Shulman, G.L., 2001. A default mode of brain function. Proc. Natl. Acad. Sci. U.S.A. 98 (2), 676-682.

Raichle, M.E., Snyder, A.Z., 2007. A default mode of brain function: a brief history of an evolving idea. Neuroimage 37 (4), 1083-1090.

Sheline, Y.I., Raichle, M.E., 2013. Resting state functional connectivity in preclinical Alzheimer's disease. Biol. Psychiatry 74 (5), 340-347.

Shen, X., Cox, S.R., Adams, M.J., Howard, D.M., Lawrie, S.M., Ritchie, S.J., Bastin, M.E., Deary, I.J., McIntosh, A.M., Whalley, H.C., 2018. Resting-state connectivity and its association with cognitive performance, educational attainment, and household income in the UK Biobank. Biol. Psychiatry 3 (10), 878-886 *Cognitive Neuroscience and Neuroimaging*.

Shine, J.M., Bissett, P.G., Bell, P.T., Koyejo, O., Balsters, J.H., Gorgolewski, K.J., Moodie, C.A., Poldrack, R.A., 2016. The dynamics of functional brain networks: integrated network states during cognitive task performance. Neuron 92 (2), 544-554.

Shirer, W.R., Ryali, S., Rykhlevskiaia, E., Menon, V., G.reicius, M.D., 2012. Decoding subject-driven cognitive states with whole-brain connectivity patterns. Cereb. Cortex 22 (1), 158-165.
Sohrabpour, A., Lu, Y., Kankirawatana, P., Blount, J., Kim, H., He, B., 2015. Effect of EEG electrode number on epileptic source localization in pediatric patients. Clin. Neurophysiol. Off. J. Int. Fed. Clin. Neurophysiol. 126, 472–480.

Song, J., Davey, C., Poulsen, C., Luu, P., Turovets, S., Anderson, E., Li, K., Tucker, D., 2015. EEG source localization: sensor density and head surface coverage. J. Neurosci. Methods 256, 9–21.

Srinivasan, R., Tucker, DM., Murias, M., 1998. Estimating the spatial nyquist of the human EEG. Behav. Res. Methods Instrum. Comput. A J. Psychon. Soc. Inc 30 (1), 8–19.

Stam, C.J., Nolte, G., Daffertshofer, A., 2007. Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Hum. Brain Mapp. 28, 1178–1193.

Tadel, F., Baillet, S., Mosher, JC., Pantazis, D., Leahy, RM., 2011. Brainstorm: a user-friendly application for MEG/EEG analysis. Comput. Intell. Neurosci. 2011. doi:10.1155/2011/879716.

Tait, L., Özkan, A., Szul, MJ., Zhang, J., 2021. A systematic evaluation of source reconstruction of resting MEG of the human brain with a new high-resolution Atlas: performance, precision, and parcellation. Hum. Brain Mapp. 42 (14), 4685–4707.

Urigüen, J.A., Garcia-Zapirain, B., 2015. EEG artifact removal - state-of-the-art and guidelines. J. Neural Eng. 12, 1–23.

van den Heuvel, M., Mandl, R., Pol, H.H., 2008. Normalized cut group clustering of resting-state FMRI data. PLoS One 3 (4). doi:10.1371/journal.pone.0002001.

Van Essen, DC., Smith, SM., Barch, DM., Behrens, TE.J., Yacoub, E., Ugurbil, K.WU-Minn HCP Consortium, 2013. The WU-Minn human connectome project: an overview. Neuroimage 80 (October), 62–79.

Van Veen, B.D., van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9), 867–880.

Wa, C., 1983. IFCN Recommendations for the Practice of Clinical Neurophysiology Amsterdam. Elsevier, New York, NY.

Wang, H.E., Bénar, C.G., Quilichini, P.P., Friston, K.J., Jirsa, V.K., Bernard, C., 2014. A systematic framework for functional connectivity measures. Front. Neurosci. 8. doi:10.3389/fnins.2014.00405.

Wang, J., Wang, L., Zang, Y., Yang, H., Tang, H., Gong, Q., Chen, Z., Zhu, C., He, Y., 2009. Parcellation-dependent small-world brain functional networks: a resting-state FMRI study. Hum. Brain Mapp. 30, 1511–1523.

Waskom, M., 2021. Seaborn: statistical data visualization. J. Open Source Softw. 6 (60), 3021.

Wendling, F., Ansari-Asl, K., Bartolomei, F., Senhadji, L., 2009. From EEG signals to brain connectivity: a model-based evaluation of interdependence measures. J. Neurosci. Methods 183, 9–18.

Zalesky, A., Fornito, A., Harding, I.H., Cocchi, L., Yücel, M., Pantelis, C., Bullmore, E.T., 2010. Whole-brain anatomical networks: does the choice of nodes matter? Neuroimage 50, 970–983.
Sohrabpour, A., Lu, Y., Kankirawatana, P., Blount, J., Kim, H., He, B., 2015. Effect of EEG electrode number on epileptic source localization in pediatric patients. Clin. Neurophysiol. Off. J. Int. Fed. Clin. Neurophysiol. 126, 472–480.

Song, J., Davey, C., Poulsen, C., Luu, P., Turovets, S., Anderson, E., Li, K., Tucker, D., 2015. EEG source localization: sensor density and head surface coverage. J. Neurosci. Methods 256, 9–21.

Srinivasan, R., Tucker, DM., Murias, M., 1998. Estimating the spatial nyquist of the human EEG. Behav. Res. Methods Instrum. Comput. A J. Psychon. Soc. Inc 30 (1), 8–19.

Stam, C.J., Nolte, G., Daffertshofer, A., 2007. Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Hum. Brain Mapp. 28, 1178–1193.

Tadel, F., Baillet, S., Mosher, JC., Pantazis, D., Leahy, RM., 2011. Brainstorm: a user-friendly application for MEG/EEG analysis. Comput. Intell. Neurosci. 2011. doi:10.1155/2011/879716.

Tait, L., Özkan, A., Szul, MJ., Zhang, J., 2021. A systematic evaluation of source reconstruction of resting MEG of the human brain with a new high-resolution Atlas: performance, precision, and parcellation. Hum. Brain Mapp. 42 (14), 4685–4707.

Urigüen, J.A., Garcia-Zapirain, B., 2015. EEG artifact removal - state-of-the-art and guidelines. J. Neural Eng. 12, 1–23.

van den Heuvel, M., Mandl, R., Pol, H.H., 2008. Normalized cut group clustering of resting-state FMRI data. PLoS One 3 (4). doi:10.1371/journal.pone.0002001.

Van Essen, DC., Smith, SM., Barch, DM., Behrens, TE.J., Yacoub, E., Ugurbil, K.WU-Minn HCP Consortium, 2013. The WU-Minn human connectome project: an overview. Neuroimage 80 (October), 62–79.

Van Veen, B.D., van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9), 867–880.

Wa, C., 1983. IFCN Recommendations for the Practice of Clinical Neurophysiology Amsterdam. Elsevier, New York, NY.

Wang, H.E., Bénar, C.G., Quilichini, P.P., Friston, K.J., Jirsa, V.K., Bernard, C., 2014. A systematic framework for functional connectivity measures. Front. Neurosci. 8. doi:10.3389/fnins.2014.00405.

Wang, J., Wang, L., Zang, Y., Yang, H., Tang, H., Gong, Q., Chen, Z., Zhu, C., He, Y., 2009. Parcellation-dependent small-world brain functional networks: a resting-state FMRI study. Hum. Brain Mapp. 30, 1511–1523.

Waskom, M., 2021. Seaborn: statistical data visualization. J. Open Source Softw. 6 (60), 3021.

Wendling, F., Ansari-Asl, K., Bartolomei, F., Senhadji, L., 2009. From EEG signals to brain connectivity: a model-based evaluation of interdependence measures. J. Neurosci. Methods 183, 9–18.

Zalesky, A., Fornito, A., Harding, I.H., Cocchi, L., Yücel, M., Pantelis, C., Bullmore, E.T., 2010. Whole-brain anatomical networks: does the choice of nodes matter? Neuroimage 50, 970–983.