Reliable evaluation of functional connectivity and graph theory measures in source-level EEG: How many electrodes are enough?

Objective: Using EEG to characterise functional brain networks through graph theory has gained significant interest in clinical and basic research. However, the minimal requirements for reliable measures remain largely unaddressed. Here, we examined functional connectivity estimates and graph theory metrics obtained from EEG with varying electrode densities. Methods: EEG was recorded with 128 electrodes in 33 participants. The high-density EEG data were subsequently subsampled into three sparser montages (64, 32, and 19 electrodes). Four inverse solutions, four measures of functional connectivity, and five graph theory metrics were tested. Results: The correlation between the results obtained with 128-electrode and the subsampled montages decreased as a function of the number of electrodes. As a result of decreased electrode density, the network metrics became skewed: mean network strength and clustering coefficient were overestimated, while characteristic path length was underestimated. Conclusions: Several graph theory metrics were altered when electrode density was reduced. Our results


Introduction
For nearly a century, the electroencephalogram (EEG) has held an undisputed position in clinical neurological diagnostics. Facilitated by the advent of accessible high-performance computers over the past two decades, the field faces a paradigm shift from the expert visual analysis and classification of well-characterised signal compounds into automated interpretation of complex, large-scale temporal and spatial signal dynamics. Among the most promising approaches are methods for capturing functional brain connectivity, i.e., the interactions over time between spatially separated modules in the brain (Rossini et al., 2019;van Diessen et al., 2015). From the EEG signals, we can extract statistical dependencies between brain regions independent of the anatomical connections between them. Through functional connectivity, brain networks at various scales may be modelled with the application of graph theoretical concepts (Rubinov and Sporns, 2010). Such approaches have already gained clinical interest, including localisation of the seizure onset zone in epilepsy (Baroumand et al., 2021) and as a biomarker for cognitive changes in various neurological diseases (Hassan et al., 2017;Hatlestad-Hall et al., 2021a;Rodríguez-Cruces et al., 2020;Vlooswijk et al., 2011). However, the methodological complexity and required computational power of these approaches make their current use in clinical neurology scarce.
The main advantage of EEG lies in its ability to represent electrical activity from the brain on a millisecond scale. However, EEG is limited in its spatial resolution: Firstly, the electrical potentials recorded from scalp-located electrodes reflect neuronal signal transmission mainly in the pyramidal neurons of the cortex (Cohen, 2017); and secondly, the current flow between the signal source and the electrode is filtered through the different tissues of the head, i.e., the signals are altered and mixed by the volume conduction effect (Brunner et al., 2016). While reports suggest that subcortical brain structures may also contribute to the EEG signal (Seeber et al., 2019), its sensitivity to the latter remains a persistent caveat of EEG. To mitigate this effect, source reconstruction, i.e., the estimation of source-space signals given the individual's head anatomy (Hassan and Wendling, 2018;Michel et al., 2004;Schoffelen and Gross, 2009), is warranted. In this regard, it is important to note that the anatomies of both the human brain (e.g., convolution, grey/white matter distribution) and the skull (e.g., thickness, curvature) have profound individual variation. As a consequence, scalp EEG signals are essentially passed through an individual mixing and attenuating filter. Thus, while anatomical templates may sometimes be used, individual anatomical images obtained with magnetic resonance imaging (MRI) or computerised tomography (CT), can greatly enhance the source reconstruction by introducing individualised information about tissue distribution (Céspedes-Villar et al., 2020;Douw et al., 2018).
Electrophysiological source reconstruction relies on solving the inverse problem of localising the generators of the measured signals. This problem is ill-posed, as it requires applying a set of additional constraints to yield a unique solution, including restricting the solution to physiologically plausible sources (e.g., cortical grey matter), presuming a set of independent, uncorrelated sources (used e.g., in beamformer solutions), and discarding solutions that require an excessive amount of energy (e.g., minimum norm estimates). The inverse solution algorithms can broadly be categorised into those based on localised (current dipoles) or distributed sources (Michel and Brunet, 2019). The former model is based on the a priori assumption that a limited number of rather focal brain sites generate the electrical potential measured at the scalp, and it is typically used, e.g., for modelling responses in sensory brain regions where the number of activated, point-like sources can be assumed to be limited. In contrast, the distributed source models assume superficial, distributed current estimates, and are used, e.g., in the present work with little a priori information of the activated sites. For a comprehensive review on source reconstruction methods, the reader is referred to, e.g., He et al. (2018). As constraints and assumptions regarding the underlying sources vary amongst inverse solution approaches, it is crucial to understand the effect of the chosen source reconstruction method on estimates of functional connectivity.
Antecedent to the importance of selecting an appropriate mathematical approach for the source reconstruction process are parameters inherent to the raw scalp EEG data. Among these, and perhaps the most challenging to accommodate in clinical settings, is the spatial coverage of the EEG recording (i.e., the number and spatial distribution of electrodes). Importantly, reports from both empirical (Lantz et al., 2003) and simulated (Song et al., 2015) EEG data suggest that source localisation is considerably enhanced by employing dense electrode arrays. However, we know relatively little about how the electrode density affects estimates of source-space functional connectivity and subsequently derived graph theory metrics, and whether the potential effect depends on the employed source reconstruction technique or functional connectivity measure.
These are important issues, considering the increasing use of graph theory concepts, such as the small world index (Bassett and Bullmore, 2017), in describing functional brain networks in various applications. Furthermore, for modern analysis methods to be implemented in clinical neurology, we must find the optimal balance between high-quality research and realistic clinical practice. Most clinics are incapable of conducting high-density EEG due to the time and resources required.
Here, we investigated the effect of electrode density on different functional connectivity and graph theory-based metrics while employing various source reconstruction methods. Specifically, we aimed to establish if sparse and clinically feasible electrode montages (64, 32 and 19 electrodes) can replicate the results obtained with a high-density 128-electrode montage without a significant loss of information. We calculated functional connectivity in source-space, based both on phase (phaselocking value, PLV; Bruña et al., 2018;Lachaux et al., 1999) and amplitude (amplitude envelope correlation, AEC; Brookes et al., 2011a;Hipp et al., 2012), for the different electrode densities. From the resulting functional connectivity matrices, we calculated several frequently reported graph theory metrics and examined whether they differed between the outlined methodological scenarios.
Our principal hypothesis was that the investigated metrics would demonstrate instability dependent on the sparseness of the electrodes. Yet, the effects were expected to be qualitatively similar across individuals, resulting in relatively strong correlations between the graph metrics derived from the different electrode montages.

EEG data and spatial subsampling
An EEG dataset comprising 33 participants (mean age: 55.7 ± 5.9 years; 57.6% female; 87.9% righthanded) was analysed in the present study. The participants comprised both chronic phase focal epilepsy patients and healthy controls; the sample is detailed in previous reports (Hatlestad-Hall et al., 2021a, 2021b. The data were originally recorded with 128 electrodes organised in a radial positional scheme during five minutes of wakeful resting using a BioSemi ActiveTwo system (sampling frequency: 2048 Hz). The data preprocessing included re-referencing to an average signal, resampling to 512 Hz, removal of channels with high levels of noise and independent components related to ocular or muscular activity, band-pass filtering (1-45 Hz), and visual inspection of the cleaned four-seconds non-overlapping data segments. The procedures are detailed in Hatlestad-Hall et al. (2021b). To enhance the accuracy of the subsequent source reconstruction (Homölle and Oostenveld, 2019), individual electrode positions were obtained with a Structure Sensor 3D scanner (Occipital Inc., Boulder, CA). The preprocessed 128-channel data were then subsampled into lowerdensity montages by extracting the electrodes from the original montage approximately corresponding to the standard 10-10 and 10-20 montages (Oostenveld and Praamstra, 2001) for 64, 32 and 19 electrodes. Figure 1 shows an overview of the applied montages, and Table 1 contains the across-montage electrode mappings. In this scenario, the data associated with each electrode is identical across the montages; only the absolute number of electrodes is reduced in order to

Source reconstruction
Four source reconstruction techniques were investigated in the present work. The methods are extensively detailed elsewhere, and thus introduced only briefly here.
The linearly constrained, minimum variance beamformer (LCMV, van Veen et al., 1997) is a spatial filtering-based, distributed dipole solution. Each estimated source is independently reconstructed based on the premise that the sensor-space measurements are composed of activity from a specific dipolar source and interference from other sources. Consequently, signal compounds not coming from the defined sources are effectively removed, producing a noninvertible result (i.e., not all sensor-space activity is represented in the source estimation).
Minimum norm estimate (MNE, Hämäläinen and Ilmoniemi, 1994) returns the source-level activity that most accurately represents the sensor-space data (i.e., with the smallest error), while minimising the overall energy. Thus, the solution represents the most efficient source configuration capable of generating the measured sensor-space activity. In contrast to the LCMV method, the MNE solution is invertible. The MNE solution has an inherent bias towards superficial sources to which the sensors are most sensitive (Dale et al., 2000), and depth-weighting is usually applied (Lin et al., 2006).
Low-resolution electrical tomography (LORETA) is a modified MNE solution specifically designed for low-resolution EEG (Pascual-Marqui et al., 1994). It adds an extra constraint to the solution, forcing it to be spatially smooth. Both the sLORETA (standardised, Pascual-Marqui, 2002) and eLORETA (exact, Pascual-Marqui, 2007) algorithms return standardised current density estimates, effectively resolving the bias towards superficial generators in the original MNE solution. Importantly, under the right circumstances, these approaches provide zero localization error (i.e., the maximum of the reconstructed activity overlaps with the physiological source). The difference between the two LORETA algorithms lies in how the source-space activity is reconstructed: in sLORETA, the activity is reconstructed independently for each source (i.e., similarly to beamforming), yielding a noninvertible solution; whereas in eLORETA, a single inverse operator is used to reconstruct the activity in all the sources simultaneously (i.e., similarly to the original MNE), resulting in an invertible solution.

Forward modelling
The source and forward models used for source reconstruction were identical across the four applied inverse solutions and are described in detail in Hatlestad-Hall et al. (2021b). In brief, the source model consisted of a uniform grid with 10 mm spacing defined in MNI space and subsequently transformed to the individual anatomy using an individual T1-weighted MRI, which was available from each participant. Only source positions in cortical regions according to the automated anatomical labelling atlas (AAL, Tzourio-Mazoyer et al., 2002) were evaluated for this purpose, resulting in 1210 sources. The conduction model was based on a three-layer (inner skull, outer skull, and scalp) model, and the boundaries were defined individually using the unified segmentation algorithm (Ashburner and Friston, 2005). The subject-specific electrode positions were extracted from the 3D scan of the subject's head while wearing the EEG cap. Finally, the forward model was solved numerically using the symmetrical boundary element method as implemented in OpenMEEG (Gramfort et al., 2010).

Phase-locking value
The phase-locking value (PLV) is defined as the instantaneous phase difference between two signals (Bruña et al., 2018;Lachaux et al., 1999), and measures functional connectivity under the assumption that brain areas generating signals whose phases evolve dependently are functionally connected (Rosenblum et al., 1996). The phase in the oscillatory activity of the brain is related to cycles of activation and deactivation of neuronal populations (Buzsáki et al., 2012). Consequently, a dependency between the phases of signals implies a relationship between their activation and deactivation cycles and can be considered an indication of functional coupling. In the case of PLV, the level of coupling is estimated from the stability of the phase difference between the signals: two signals with a constant phase difference are very likely to be coupled (PLV is 1), while two signals with a random phase difference are likely independent (PLV is 0); any other case ranges between 0 and 1.

Amplitude envelope correlation
The amplitude envelope correlation (AEC) is defined as the correlation between the power variations in two signals (Brookes et al., 2011a;Hipp et al., 2012). If two regions exhibit an increase or decrease in their amplitude concurrently, their activity is likely connected, and their AEC is high. In contrast, the AEC is low if the power variations are unrelated. The measure is specified between 0 (uncorrelated) and 1 (perfect correlation).

Leakage corrections
Functional connectivity estimates are known to be affected by source leakage (Stam et al., 2007), which spuriously increases the estimated values. This effect relates to the spatial resolution of the source reconstruction, and it is thus amplified when the number of EEG electrodes is limited. Since source leakage can be modelled as a linear mixing of true source activities, removing the instantaneous projection of one signal over the other can eliminate this effect. In this work, the effect of source leakage was controlled separately for PLV and AEC. For the PLV, the instantaneous (i.e., zero-lag) coupling was eliminated by considering only the (corrected) imaginary component of the PLV, rather than the original value (Bruña et al., 2018). For the AEC, orthogonalisation was used (Brookes et al., 2011a). Although correcting for spatial leakage between sources is generally recommended (Colclough et al., 2016), removing instantaneous interactions has been shown to also eliminate actual physiological information (Kovach, 2017) and to reduce the test-retest reliability of the functional connectivity metrics (Garcés et al., 2016). Therefore, both the original and leakagecorrected versions of the estimates were investigated here. In the following, the leak-corrected measures are termed ciPLV (corrected imaginary part of the PLV) and lcAEC (leak-corrected AEC).

Measures calculations
As both PLV and AEC are known to exhibit frequency-specific traits (Brookes et al., 2011b;Bruña et al., 2018), the functional connectivity and the derived graph metrics were calculated separately for two frequency bands; alpha (α; 8-12 Hz) and beta (β; 12-30 Hz). The PLV and AEC measures were calculated between each pair of source positions, resulting initially in a 1220-by-1220 matrix for each measure. By computing the root-mean-squared value of all the connections between each pair of cortical areas defined in the AAL atlas, these matrices were condensed to represent the functional connectivity between each pair of 80 (40 in each hemisphere) cortical areas (Bruña and Pereda, 2021). This resulted in an 80-by-80 matrix per participant, electrode density, frequency band, and functional connectivity measure. Please see the Supplementary Material on functional connectivity for a detailed description of the PLV and AEC (and their leak-corrected versions) calculations.

Functional network construction
An important supposition in network analysis of functional connectivity is that all the areas of the brain are constantly active and thus, to some degree, are involved in neurophysiological processes.
To reflect this, the estimated adjacency matrix should be fully connected, i.e., there should be no disconnected brain areas, or nodes. However, as most graph theory metrics were designed to characterise sparse networks (i.e., networks in which not all nodes are interconnected), a predefined proportion of the network connections, or edges, is commonly discarded. In most cases, the weakest edges are removed, as these are putatively more likely to be spurious (Fornito et al., 2012).
Considering that there exists no consensus regarding the size of this proportion (Garrison et al., 2015;Hatlestad-Hall et al., 2021b;van Wijk et al., 2010), the calculation of graph metrics in this work was initially performed on a range of pre-defined network density (i.e., proportion of considered connections) thresholds (0.5 to 0.9, with 0.1 increments). As the network density had a negligible effect on reproducibility across electrode montages, the remaining analyses were conducted using a network density threshold of 0.7. To preserve the full connectedness characteristic in the thresholdimposed networks, the sparse matrices were obtained by incrementally adding edges by descending weight order to the minimum spanning tree backbone of the inverse dense functional connectivity matrix until the specified sparsity level was reached (Hatlestad-Hall et al., 2021b). The minimum spanning tree is defined as the subnetwork that connects all nodes while minimising edge weights and avoiding loops Tewarie et al., 2015).
For the calculation of the small world index, the clustering coefficient and the characteristic path length were normalising against the mean of the corresponding measure calculated over 25 pseudorandom null model networks re-wired from the dense functional connectivity matrix, while preserving weight distributions (Rubinov and Sporns, 2011).

Graph theory metrics
In the present work, five commonly reported global graph metrics were considered, including the average node strength, the characteristic path length, the clustering coefficient, and the assortativity coefficient. Also, the small world index, here defined as the ratio between the normalised measures of the clustering coefficient and the characteristic path length (Humphries and Gurney, 2008;Watts and Strogatz, 1998), was examined. In brief, the characteristic path length is the average shortest path between all pairs of nodes in the network, whereas the clustering coefficient is the fraction of triangles around an individual node (i.e., the fraction of the node's neighbours that are also neighbours of each other). The small world index is a measure of the balance between local segregation (clustering) and global integration (long paths). Assortativity, also known as the assortative mixing coefficient, is a measure of a network's nodes' preference to connect to other nodes with a similar degree, i.e., number of edges (Newman, 2002). It is defined as the Pearson correlation between all nodes on opposite ends of an edge (Leung and Chau, 2007;Rubinov and Sporns, 2010). The mean network strength is equal to the mean node strength, defined as the sum of the weights of the edges connected to the node. For comprehensive accounts of graph metrics and their mathematical definitions, the reader is referred to Rubinov and Sporns (2010) and Newman (2008).

Reproducibility of graph theory metrics across montages
To test the reproducibility of estimated EEG-derived graph theory metrics across montages of varying electrode density, we used the intraclass correlation coefficient (ICC, McGraw and Wong, 1996). The ICC reflects the ratio of between-subject variance to within-subject variance, and is normalised between 0 and 1, where a value above 0.5 indicates greater between-subject variance (Bassett et al., 2011). ICCs above 0.5 are commonly suggested to indicate moderate reproducibility, and values above 0.75 good reproducibility (Koo and Li, 2016). In this case, where the measurements are independent and the montage factor has fixed levels (type C-1), the ICC is defined as (McGraw and Wong, 1996): The ICC was calculated separately for each unique combination of source reconstruction method, functional connectivity measure, network density, frequency band, and graph theory metric. To assess if the reproducibility of a graph metric across electrode densities depends on any of the first three factors, we carried out two separate repeated measures ANOVAs for the two frequency bands (α and β). The measurements were represented by the ICCs associated with the five different graph metrics.
Subsequently, one-sample t tests (one-tailed) were conducted to test whether each estimated ICC was greater than 0.5, i.e., indicating moderate consistency (Koo and Li, 2016). As the previous analysis revealed the ICCs independence of the network density level, the ICCs were collapsed across the five network density levels for this analysis. The tests' p values were corrected for multiple comparisons using the false discovery rate (FDR) procedure (Benjamini and Yekutieli, 2001).

Correlation of graph theory metrics between the baseline and subsampled EEG montages
In a second series of analyses, we narrowed the scope to investigate the interrelationship between graph theory metrics derived from the baseline, high-density montage (128 electrodes) and the subsampled montages (64, 32, and 19 electrodes). For these analyses, an a priori assumption was made about the high-density montage having a superior signal-to-noise ratio compared to those composed of sparser electrode constellations. Pearson correlation coefficients were calculated between arrays of values obtained from the baseline and subsampled montages separately for every combination of source reconstruction method, functional connectivity measure, and graph theory metric (strength, path length, clustering coefficient, assortativity, and small world index).
Considering previous findings that graph metrics derived from adjacency matrices with varying sparsity levels exhibit different absolute values while retaining their relative difference (Hatlestad-Hall et al., 2021b;van Wijk et al., 2010), and that the ICC in this study was not significantly dependent on the network density level, the analyses were carried out using adjacency matrices with a set network density of 0.7, as explained above. To test for main effects of source reconstruction method, functional connectivity measure, and montage, on the magnitude of the correlation between graph theory metrics derived from the baseline and subsampled electrode densities, a three-way repeated measures ANOVA was conducted for each of the frequency bands (α and β). The five graph theory metrics were defined as observations. Finally, we analysed the correlations between the baseline and subsampled montages broken down by the individual graph metric. As the correlation coefficients were independent of the functional connectivity statistic in the repeated measures ANOVA, the graph theory metrics were averaged across functional connectivity measures for this analysis. Pearson correlation coefficients were calculated separately for each combination of graph theory metric and source reconstruction method for each of the frequency bands (α and β). The correlation p values were computed nonparametrically with 1 million permutations and then corrected for multiple comparisons with the FDR procedure (60 comparisons per frequency band). Lastly, the correlation coefficients were compared pairwise between the source reconstruction methods using t tests with Bonferronicorrected p values.

Relative effect of the inverse solution and electrode density on graph metrics
In a second analysis, we tested the effect of source reconstruction method and electrode density on derived graph theory metric values, using repeated measures ANOVAs, implementing these as within-subject factors. The repeated measures models were defined separately for each frequency band, graph metric, and functional connectivity measure.

Reproducibility of nodal functional connectivity across montages
Mathematically, equally sized functional connectivity matrices with different distributions may give rise to similar estimates of several global graph theory metrics. Thus, we conducted a series of analyses to delineate the consistency and reproducibility of functional connectivity estimates across electrode density levels and source reconstruction methods. As the different functional connectivity measures give estimates with variable means, ranges, and variance, and are thus not directly comparable, the measures were analysed separately. The analysis was two-fold. First, we calculated the ICC across all four electrode density levels for each source reconstruction method, and secondly, we calculated the Pearson correlation between functional connectivity estimates obtained from the baseline electrode density level and the subsampled levels. The coefficients were calculated separately for each network edge, then averaged for each node (AAL region). Considering the symmetricity of the electrode montages, and that we did not expect systematic differences in estimated connectivity between hemispheres, the nodal coefficients were averaged across the hemispheres. For the ICCs, we conducted a series of nonparametric significance tests based on 1 million permutations to evaluate if each coefficient was greater than 0.5. The resulting p values were corrected for multiple comparisons using the FDR procedure (16 comparisons per frequency band).

Statistics and software
The network analysis was conducted with functions available in the Brain Connectivity Toolbox (BCT; ver. 2019-03-03, Rubinov and Sporns, 2010) and with in-house MATLAB code (see Data availability statement). Null model networks were generated with the BCT function null_model_und_sign. To quantify the proportion of explained variance by the repeated measures model factors, we used the effect size estimate eta squared, η 2 (Lakens, 2013). The corresponding effect size estimate used for Pearson correlations was r 2 . For post hoc comparisons based on the repeated measures models' estimated marginal means, effect size estimates are provided in the form of Hedges' g with the mean standard deviation of the samples as the standardiser (g AV ). In the violin plots (see Figures), the marker represents the median value, the vertical black line represents the interquartile range, the height of the violin represents the range, and the violin shape represents a kernel density estimate of the data.

Ethical approval
Ethical approval for the study was granted by the Regional Committees for Medical Research Ethics of Norway (reference no. 2015/2470).

Participant informed consent
Following the Declaration of Helsinki, all participants provided informed written consent.

Data availability
The functional connectivity and graph metric data that support the findings of this study, and the MATLAB code used for calculating graph metrics and statistics, are available upon request.

Graph theory metrics
Figures 2 to 5 display the graph theory metric values obtained from the various constellations of functional connectivity measure, source reconstruction method, and montage density, collapsed across the α and β bands.

Reproducibility of graph theory metrics across EEG electrode montages and source reconstruction methods
The intraclass correlation coefficient (ICC) was computed across electrode montages separately for every combination of source reconstruction method, functional connectivity measure and network density level.
There was no significant main effect of the functional connectivity measure in either frequency band.
Subsequently, we analysed the correlations between the baseline 128-electrode montage and the

Relative effect of the source reconstruction method and electrode density on graph metrics
The repeated measures ANOVA revealed significant main effects of both within-subject factors (source reconstruction method and electrode density; subscripted SR and ED, respectively, below), as well as interaction effects, across most of the functional connectivity measures and graph theory metrics in both frequency bands (α, β). All η 2 estimates broken down by graph theory metric, and frequency band are presented in Figure 8.

Reproducibility of nodal functional connectivity across different electrode densities
In general, all combinations of source reconstruction method and functional connectivity measure

Discussion
In the last decade, there has been a significant growth in applying graph theory to characterise functional brain networks derived from EEG. To date, however, there is no consensus regarding the necessary methodological implementations. In this study, we demonstrate that methodological choices have a nontrivial effect on the estimation of functional connectivity and graph theory metrics. Particular attention was here paid to the number of EEG electrodes used in the recordings, i.e., the spatial sampling density, as this element may be the most pertinent to the potential clinical value of these approaches. Our results show that several widely employed graph theoretical concepts, such as the clustering coefficient and the small world index, are inconsistent between high-density data and data sampled with fewer electrodes: the fewer the electrodes used, the larger the difference to the baseline montage of 128 electrodes. Furthermore, we found that the reproducibility of graph theory metrics across different electrode density arrays depends on both the applied source reconstruction method and the statistic used to estimate functional connectivity.
Although it is essential to emphasise that reproducibility indices do not imply validity, it might be argued that reliability across spatial electrode density is crucial for guiding the methodological parameters that frame a network study.

Consistency across EEG electrode densities
The same fundamental principles govern spatial frequency as frequency in the time domain: the Nyquist theorem of discretization asserts that the sampling rate must be at least twice (2.5 to account for the possibility of phase-locking) as frequent as the maximum frequency to be described in order to avoid aliased signals (i.e., an artificial increase in energy at lower frequencies caused by contamination from higher frequencies). Unlike in the time domain, the spatial signal is intrinsically a discrete sampling of the continuous scalp surface potential field, constrained by the density of the electrode array. Importantly, prior to digitization, the temporal signal can be low-pass filtered to remove aliasing information, whereas the spatial signal cannot, eliminating the opportunity to rectify undersampling (Srinivasan et al., 1998). The poorly conductive skull is the only spatial anti-aliasing filter available. Previous studies suggest that an electrode separation of approximately 2-3 cm is the minimum required to obtain spatial EEG patterns without undersampling (Freeman et al., 2003;Srinivasan et al., 1998). To meet this need for an average adult human head, a minimum of 100 electrodes are required (Michel et al., 2004). However, electrode arrays with inter-electrode spacing between 0.5 and 1 cm have been shown to yield even more functional information (Petrov et al., 2014). This suggests that the threshold for preventing undersampling in scalp EEG may be more stringent than previously assumed.
Several studies have demonstrated advantages of using high-density electrode arrays. Epileptic source localisation accuracy rose with the number of electrodes up to 128, after which the effect of further increases became marginal (Sohrabpour et al., 2015;Song et al., 2015). Montages with fewer than 64 electrodes yielded poor outcomes. Furthermore, research utilising source-level functional connectivity (Staljanssens et al., 2017a(Staljanssens et al., , 2017b and high frequency oscillations (Kuhnke et al., 2018) to localise seizure onset zones in epilepsy patients revealed a considerable advantage of applying high density electrode montages. Here, we demonstrate that compared to the 128-electrode montage, the functional connectivity estimates and graph theory metrics show fairly strong correlations for the 64-electrode montage, lesser correlations for the 32-electrode montage, and the weakest correlations for the 19-electrode montage. In keeping with the source localisation accuracy from other studies, the correlation magnitude increase per electrode is substantially greater when increasing from 19 to 32 than when increasing from 32 to 64. These results may indicate that certain functional brain information is shared across spatial frequency levels and that the extent of spatial aliasing distortion follows a trajectory.
In our results, the intraclass correlation coefficients (ICC) measured across electrode densities were slightly higher for node-wise functional connectivity estimates compared to graph theory metrics.
For the former, ICCs (across nodes) ranged from around 0.5 to 0.9, whereas for the latter, the lowest values were close to 0.2. This implies that graph theory metrics, irrespective of the employed source reconstruction method and functional connectivity measure, are more susceptible to electrode sparsity than their precursor in the analysis pipeline, i.e., the applied functional connectivity measures. However, it remains uncertain what, in this context, constitutes a reasonable ICC size. For example, it is known that the ICC is reliant on the range and variance in the data (Lee et al., 2012), i.e., parameters that were not thoroughly studied in this paper. Yet, the variance associated with the majority of graph theory metrics was larger for sparser montages.

Graph theory metrics across electrode densities
Basic and complex topological characteristics of a functional brain network can be efficiently summarised using graph theory metrics. However, understanding the impact of electrode density on EEG-based graph theory metrics is crucial for their informative use. In the current study, we found that several of the widely used graph theory metrics, calculated in source-space, varied as a function of the electrode density. The average node strength displayed a predictable pattern across all functional connectivity measurements, with sparser electrode montages overestimating its magnitude. This also held true for the clustering coefficient, whereas the characteristic path length was underestimated for montages with low spatial density. All three observations are probably related: the lower the spatial sampling, the less activity is unique to each of the reconstructed sources. Regarding node strength, this results in an overestimation of connectivity between neighbouring nodes, which in turn increases the number of local edges retained when a threshold is applied to the network. As a consequence, more nodes form local triplets in the network, which inflates the clustering coefficient. In direct relation to the overestimation of local connectivity, the characteristic path length of the network is decreased; more short paths are retained than longer ones. This does not necessarily imply that functional networks constructed from sparse electrode montages are without long paths. However, due to the overrepresentation of within-module pathways, the analysis will be less sensitive to the presence of longer paths. This issue should be studied further using network analysis at the node level.
Our data is less decisive regarding the small world index and network assortativity. The small world index is intriguing because it is a normalised metric that is fairly extensively employed. In principle, the small world index is more directly comparable among montages due to the normalisation that accounts for the confounding variances in clustering coefficient and characteristic path length.
Despite this, our findings suggest that the small world index increases as a function of electrode sparsity. This was the case for all inverse solutions and functional connectivity measures, with the exception of the corrected imaginary part of the PLV (ciPLV). These results contrast with a prior study that addresses this same subject and suggested that the small world index was similar in lower and higher electrode densities (Miraglia et al., 2021). It should be emphasised, however, that the measure of functional connectivity applied there, lagged linear coherence, was not examined here, and in their study, the most dense montage comprised 59 electrodes. The small world index quantifies the extent to which a network consists of local, highly connected modules while maintaining longer, direct connections between modules. Our results demonstrate that functional networks created from sparse electrode montages skew this ratio between local and global connectedness, even for normalised metrics. In contrast, the electrode density had no effect on the network's assortativity. Assortativity, which is a measure of a network's nodes' preference to connect to other nodes with a similar degree, was the only graph theory measure investigated here that was not notably more dependent on the montage than the applied inverse solution.

Source reconstruction methods and functional connectivity measures
Interestingly, when targeting a particular neural source (i.e., well-defined focal activity), reasonable localisation accuracy can be achieved with as few as 19 electrodes, but only by employing beamforming techniques to solve the inverse problem (Nguyen-Danse et al., 2021). This is consistent with the current findings, in which the beamformer approach (LCMV) displayed greater reproducibility across the majority of functional connectivity measures and graph theory metrics than the other methods. Beamformers reconstruct source-level activity source-by-source, isolating focal activity as a result, whereas other techniques that reconstruct all sensor-level activity (i.e., produces invertible solutions) may be more sensitive to spatial aliasing. The remaining inverse solutions tested here exhibited generally significant correlations between graph theory metrics obtained from the 128-electrode montage and the 64-and 32-electrode tiers. In assessing the 19electrode montage, however, all inverse solutions except for LCMV failed to demonstrate significant correlations with the baseline density for several key metrics, including the small world index. This was true for the two frequency bands studied here.
The various functional connectivity measures differ in their underlying assumptions and the magnitude of their estimates (van Diessen et al., 2015). Phase and amplitude coupling, both physiologically plausible in the brain, were considered in this study. Their values seem to be independent, at least up to a certain level; however, phase and amplitude coupling might reflect different processes, as the energy required to alter the phase of a system is negligible in comparison to the energy required to increase its amplitude (Siems and Siegel, 2020). Graph theory metrics will inherently vary based on the functional connectivity measure from which they are derived, due to these fundamental underlying differences. Yet, in our results, the choice of functional connectivity measure had no significant effect on the ICC between the four decreasing density montages. This suggests a similar degree of reproducibility among the statistics under consideration. However, regarding the dependency of graph theory metrics on electrode density, the PLV, AEC, and leakcorrected AEC all demonstrated comparable patterns, whereas the ciPLV did not. The ciPLV statistic measures the portion of the PLV that is not due to zero-lag connectivity (Bruña et al., 2018), and phase synchronisation metrics removing zero-lag have demonstrated low test-retest reliability of resting-state activity in prior research (Colclough et al., 2016;Garcés et al., 2016). This could imply that these measures may neglect some physiologically plausible connectivity (Kovach, 2017), which could have an effect on the estimation of network-based properties such as those considered in the current work. In addition, ciPLV typically returns very low values of synchronisation (here implied by the comparatively low mean node strength), thus the network thresholding procedure performed prior to the estimation of graph theory metrics may be more susceptible to noise than in other measures. Overall, this may indicate that this type of network analysis is less reliable for ciPLV, although the metric can still provide useful information in other contexts (see, for example, Bruña, 2019).

Clinical relevance
Ultimately, a context-specific cost-benefit analysis must identify the ideal electrode array densities in practice. Increasing the number of recording sites results in higher-dimensional data at the cost of both material and human resources. Recently, Stoyell et al. (2021) reported a series of epileptic patient cases in which high-density EEG significantly aided clinical diagnostics, including localisation of epileptogenic zones in close proximity to eloquent motor/language cortex and identification of high-frequency ripples. In generalising from these cases, the authors highlighted future clinical applications for high-density EEG. According to current recommendations, a minimum of 64 electrodes should be utilised when source localisation is important (Seeck et al., 2017). By demonstrating a relatively strong correspondence between functional connectivity and graph theory metrics derived from 64 and 128 electrodes, our current findings support this recommendation of expanding the use of high-density (at least 64 electrodes) EEG into clinical practice. Importantly, recent research has identified a range of potential biomarker applications for EEG-based brain network analysis, such as the identification of disease-specific patterns (for a review, see Stam, 2014) and neurodegeneration risk (for a review, see Rossini et al., 2020). Our findings suggest that low-density montages may lack the spatial sensitivity required to capture the key signal features needed in these applications. Also, significant and distinct associations between the functional brain network topology and cognitive functions, both in the presence and absence of brain pathology, as demonstrated by previous research from our group (Hatlestad-Hall et al., 2021a) and others (Hassan et al., 2017;Langer et al., 2012), might not be detectable with sparse electrode arrays. To determine the sensitivity of such biomarkers under sparse electrode conditions, further investigations are required.
A frequent objection against the utility of high-density EEG in clinical practice, in addition to the resource-intensive application of more dense electrode montages, is the sophisticated data processing and analyses required for source reconstruction. Indeed, reports frequently recommend the use of detailed forward models based on individual anatomy (from MRI and/or CT) and precise representations of electrode positions relative to the head surface (e.g., Liu et al., 2018;Taberna et al., 2019;Vorwerk et al., 2014). However, other studies suggest that, at least for epileptic source localisation, the importance of complicated forward modelling in practical application may be reduced (Birot et al., 2014). The same study demonstrated the feasibility of adapting template electrode positions to the head shape derived from an MRI, as opposed to individually obtaining each position. It remains to be determined whether these findings also hold true for functional connectivity estimation and network topology characterisation.

Limitations
We would like to highlight some important limitations of this study. First, we did not address the issue of electrode coverage. Previous research suggests that accurate source localisation is strongly influenced by the head surface coverage, with a particular emphasis on the inferior temporal lobes, which are rarely covered by normal EEG montages (Song et al., 2015). Second, our data is based solely on one recording session, from which all the subsampled montages were derived from. Actual repeated measurements may have offered additional and important insight into the internal consistency of the graph theory metrics; for instance, does a 128-electrode montage produce more consistent observations over time than a 19-electrode montage? Third, analysing the advantage of increasing to 256 electrodes (and beyond) would be useful in light of recent findings that suggest even denser montages (more than 128 electrodes) may capture meaningful functional information.
Therefore, in a strict sense, the densest electrode array investigated here does not represent "ground truth," limiting the conclusions that may be drawn from the current data. In addition, the conventional 10-10 and 10-20 electrode placements are not represented in the montages examined here. Fourth, the current study was limited to global network properties. Considering that local graph measurements may theoretically have clinical promise for identifying specific abnormal brain regions (nodes) participating in functional networks, we recommend additional research into the effect of electrode density on such measures. Finally, to limit the number of analyses conducted in the current investigation, only the alpha and beta frequency ranges were selected for examination.
Consequently, caution must be exercised when extrapolating the current findings to slower (delta, theta) or faster (gamma) frequency ranges.

Conclusion
Functional connectivity estimates and graph theory metrics derived from EEG data are inherently dependent on the number of electrodes used to record the original data due to variations in spatial frequency. The purpose of this study was to contribute to the investigation of this subject, which, in our opinion, has not received sufficient attention.
First, our results demonstrate that the correlation between graph theory metrics and functional connectivity estimates from the high-density baseline montage (128 electrodes) and the subsampled montages decreases with electrode sparsity. Second, when low-density electrode montages are used, short-distance network edges are overestimated in comparison to long-distance edges, resulting in an inflated small world index. Third, the choice of source reconstruction method influences the reproducibility of graph theory metrics and, to a lesser extent, estimates of functional connectivity. Further research is required to determine the precision of functional brain networkbased biomarkers at various electrode densities. Nevertheless, in order to strike a balance between resource needs and accuracy, we recommend using at least 64 electrodes when characterising functional brain networks with graph theory metrics derived from EEG.

Competing interests
CHH and IHH are shareholders in the medical technology start-up BrainSymph AS. The remaining authors report no competing interests.      those marked with an asterisk are significantly larger than 0.5 (p FDR < 0.05).

Figure Legends
Abbreviations: PLV (phase-locking value); ciPLV (corrected imaginary part of the PLV); AEC (amplitude envelope correlation); lcAEC (leak-corrected AEC).    The plots represent the distribution of correlations (r 2 ) between node-level functional connectivity obtained from the high-density electrode montage and the subsampled montages. The scatters on the right-hand axis represent the individual node correlation p FDR values. Note that p FDR values above 0.1 are not displayed.

Highlights
• Functional connectivity estimates and graph metrics based on source EEG depend on electrode density.
• Compared to high-density, low-density EEG gives skewed graph metric estimation.
• The reproducibility of graph metrics across electrode densities depends on the inverse solution.