Distortions in EEG interregional phase synchrony by spherical spline interpolation: causes and remedies

Abstract

Background

Measures of the coherence of electroencephalography (EEG) time-series recorded at spatially distant points on the scalp are often used by researchers to characterize the dynamic interactions of brain regions. In dense-array EEG recordings, one or more electrode signals often contain prominent artifact necessitating replacement of the recorded data with an estimated signal using interpolation from data in valid recordings from the surrounding electrodes. Typically the signal estimation is carried out using spherical spline interpolation (SSI; Perrin et al., 1989); however, it is shown that this can introduce an erroneous increase in coherence between signals because the estimated signal is derived from elements in common recordings from other electrode sites. Although SSI performance depends on three SSI parameters, including interpolation order m, the Legendre polynomial order n, and regularization parameter λ, clear guidelines on how to optimally choose parameters have yet to be established for ensuring the temporal features of interpolated signals are accurate.

Methods

Using a cross validation approach we compared interpolated and original signals in 64 and 128 channel resting EEG data systematically examining the effect of data interpolation across various sets of SSI parameters, frequency bands, and locations of interpolated electrodes and electrode pairs on a standard measure of phase coherence (phase locking value: PLV).

Results

We found substantial PLV inflation with SSI, which was greater in 128 channel than in 64 channel EEG data and greater for neighboring electrodes than distant electrodes. Furthermore, the amount of PLV distortions varied across the locations of interpolated signals. Performance of SSI was highly dependent on the three SSI parameters, and the effect of each parameter was dependent on the levels of other parameters.

Conclusions

Although SSI parameters commonly recommended in the literature led to notable and concerning PLV inflation, application of SSI with alternative and optimal parameter values led to accurate and nearly distortion-free phase coherence values (i.e., PLV) from estimated signals.

Background

In order to better understand how the coordination of neural activity between brain regions supports cognitive and motor functions in health and disease, scientists have examined the synchronization of electroencephalography (EEG) signals measured from different points on the scalp. EEG synchronization is thought to reflect the functional interaction (i.e., connectivity) of spatially separated groups of neurons (Varela et al., 2001) reflecting coordinated neural computation. An analysis of functional connectivity specifically requires quantitative measurement of the temporal relationships between two neural signals. The millisecond temporal resolution of EEG gives researchers the ability to richly characterize functional coordination across regions of the brain. Coherence indices are the most widely used measures of EEG interregional functional connectivity. There are several ways to compute coherence (Chorlian et al. 2009; Roach & Mathalon, 2008) but it is generally computed based on the time-frequency cross-spectra between two time-series, reflecting not only synchronization of phase but also amplitude of the EEG signal. Unlike this conventional computation of coherence, phase coherence is isolated to the phase relationships of two signals and excludes amplitude information, thus it provides a measure of functional connectivity uncontaminated by signal strength (Lachaux et al., 1999; Aviyente et al., 2011).

In functional connectivity studies using EEG signals from spatially separate scalp electrodes, the conduction of electrical currents through the brain is a serious problem that limits the utility of coherence measures in characterizing functional interactions between groups of neurons (Srinivasan et al., 2007). Because EEG coherence measures are contaminated by common signals from volume conduction in the brain (Nunez & Srinivasan, 2006), the surface Laplacian transformation has been used to better estimate the neural sources of electrical signals (Law et al., 1993; Perrin et al., 1987; 1989; Nunez et al., 1997). The surface Laplacian is a spatial band-pass filter that is applied to effectively remove spatially broad features of the data so that coherence values derived from the transform are less likely to be contaminated by common signals. However, the Laplacian as a spatial filter removes not only common signals related to volume conduction, but also genuine signals from distributed coherent cortical sources. Therefore, the surface Laplacian is a less than ideal method for studies of the functional connectivity between brain areas that are further apart than the spatial scale of the surface Laplacian. As an alternative interregional coherence of scalp EEG can still be quite useful when measured between distant (>12 cm) electrodes with average re-referencing. In a simulation study that included coherence values between 111 electrodes across various reference electrodes, coherence data using an average reference closely aligned with reference-independent potentials (Nunez & Srinivasan, 2006).

Dense placement of EEG electrodes across the scalp commonly results in one or more electrodes with entirely artifactual data due to very high impedance, poor contact with the scalp, or poor electrical conductivity due to electrode overuse or insufficient cleaning. With dense-electrode EEG recordings it is usually recommended to replace the errant data with an estimated signal by computing an interpolated EEG recording from predicted correlations with other electrode sites on the scalp. On the other hand, dense-array electrodes could result in electrical bridges between EEG electrodes due to electrolyte spreading across adjacent electrode sites, causing distortions in EEG topography as well as serious power loss in electrodes that electrically bridged to a reference electrode. It has been suggested that researchers can detect such artifacts based on similarity measures of EEG signals (Tenke & Kayser, 2001; Alschuler et al., 2014) and replace bridged signals using spherical spline interpolations from unaffected electrodes (Greischar et al., 2004). Given that interpolated data are derived from other electrode data, it seems likely that when phase relationships are computed with the interpolated data there will be artifactual connectivity involving the estimated EEG data. To date, no published papers have systematically examined the influence of interpolated EEG data on phase coherence values of functional connectivity. The goal of the present study was to examine how distortion of interregional phase coherence is related to the selection of interpolation parameters and provide guidelines for the best means of interpolating data in studies of interregional EEG phase coherence.

There are several interpolation methods employed by published papers, such as nearest neighbor methods (Shepard, 1968), thin-plate spline interpolation (Harder & Desmarais, 1972), spherical spline interpolation (SSI; Perrin et al., 1989), and 3-D spline Interpolation (Law et al., 1993). These interpolation methods were developed for the topographical mapping of scalp voltages derived from EEG recordings, but have been extended to interpolation of estimated time-series for electrodes with artifactual data. Since Perrin and his colleagues firstly introduced the SSI method for scalp EEG data, SSI has been the most widely used. In applying SSI in investigations of functional connectivity we observed artificially inflated interregional coherence for the estimated EEG signals. This is because the spline weighting matrix with common signals from other electrodes can artificially increase similarity of interpolated and original EEG signals (Biggins et al., 1991). Given that the locations of bad electrodes vary across subjects unsystematic interregional coherence distortions are introduced which may yield false condition or group effects.

Although the accuracy of SSI and the surface Laplacian transformation based on spherical spline can significantly vary depending on the choice of several parameters (Perrin et al., 1989; Soong et al., 1993; Ferree, 2006) the relationships of the parameters and the amount of interregional coherence distortion has not been systematically examined. Several studies investigated the performance of spherical spline interpolation in their comparisons with other interpolation methods (Soufflet et al., 1991; Soong et al., 1993; Fletcher et al., 1996), but none, or at most one, of the parameters were manipulated in the tests. Furthermore, investigations were largely focused on the accuracy of interpolated values in the spatial domain or at limited time points, providing little information about possible problems and the appropriate choice of SSI parameters for EEG interregional functional connectivity studies. In the present study, we investigated the effect of SSI on interregional phase coherence measures using a systematic cross-validation approach where differences of interregional phase coherence between original and interpolated electrode signals were examined. We performed tests across different numbers of channels, frequency bands, locations of interpolated electrodes, locations of paired electrodes, and sets of SSI parameters.

Methods

EEG recordings and processing

Participants from this sample were taken from two larger projects studying United States military veterans. Twenty-two male participants completed a resting session using an EEG system with 64 channels (n = 11, mean age = 56.9, SD = 11.9) or 128 channels (n = 11, mean age = 31.8, SD = 7.7). The 64 channel resting session consisted of four 45-s blocks alternating conditions with eyes open or closed. Similarly, the 128 channel resting session alternated conditions, but included six 60-s blocks. They were selected based on the following criteria: No history of traumatic brain injury and having no bad electrodes (there was one exception, but the bad electrode did not occur at a site of interest). EEG data were collected using a 64 or 128 channel Ag/AgCl Biosemi ActiveTwo system (www.biosemi.com) and an elastic cap (Electrode Cap International [ECI]). The 64 channel cap used a 10–10 configuration while the 128 channel cap was configured in an equidistant-layout. EEG data were recorded with a sampling rate of 1024 Hz, and down-sampled to 256 Hz in offline processing. In offline processing, EEG data were referenced to a linked earlobe reference signal and high-pass filtered with a cutoff frequency of .5 Hz. EEG data were segmented to 1 s epochs, and bad epochs with strong low or high frequency noise were visually inspected and removed. EEG data were denoised using independent component analysis (ICA) procedure. To obtain optimal decomposition using ICA, data dimension was reduced to the dimensionality estimated by Bayesian information criterion (Schwarz, 1978) using principle component analysis (PCA). Fast ICA algorithm (Hyvärinen & Oja, 2000) was applied to the EEG with reduced dimension to decompose EEG into independent components from various sources. Ocular, heart, muscular, and electrical noise ICs were identified based on the topography, power spectrum, and time-series of ICs (McMenamin et al., 2010). Denoised data were reconstructed by removing noise IC variances from the EEG data, then re-referenced to an average signal of all EEG channels.

Spherical spline interpolation

Spline interpolation is a type of interpolation where a special type of piecewise polynomial, called a spline, is used for data fitting and the estimation of missing values using existing values (de Boor, 1978). The spherical spline method was first introduced for scalp EEG by Perrin et al. (1989), where interpolation and the surface Laplacian can be estimated by approximating the electrode positions on a sphere. In the definition of spherical splines, let ${\overrightarrow{r}}_j$ denote the location of a measurement electrode, the spherical scalp surface (j = 1,2, …, N_e), $\overrightarrow{r}$ denote the location of an arbitrary surface point, and $V\left(\overrightarrow{r}\right)$ denote the potential at the point. Spherical splines assume that the potential at any point r on the surface of the sphere can be represented by

$$ V\left(\overrightarrow{r}\right)={c}_0+{\displaystyle \sum_{i=1}^{Ne}{c}_j{g}_m\left( \cos \left(\widehat{r},{\widehat{r}}_i\right)\right)} $$

(1)

where c ₀ and c _j are constants fit to the data and the g _m (x) is the function of the cosine of the angle between the interpolation point $\overrightarrow{r}$ and an electrode location $\overrightarrow{r}$ (variable x), which is given by

$$ {g}_m(x)=\frac{1}{4\pi }{\displaystyle \sum_{n=1}^{\infty}\frac{2n+1}{{\left(n\left(n+1\right)\right)}^m}{P}_n(x)} $$

(2)

where P _n (x) is the Legendre polynomials of order n. The Legendre polynomials are ordinary differential equations that have singular points between −1 and 1. P _n (x) is given by following recurrence relation:

$$ \left(n+1\right){P}_{n+1}(x)=\left(2n+1\right)x{P}_n(x)-n{P}_{n-1(x)} $$

(3)

with P ₀  = 1 and P ₁   = x. Legendre polynomials are useful to approximate the electrical potential at a scalp location due to a point charge within a head because the potential expands in the Legendre polynomial depending on angles between the locations (Griffiths, 2007). P _n (x) converges where −1 < x < 1. Interpolated voltage $V\left(\overrightarrow{r}\right)$ at a point r on the surface of the sphere can be computed by fitting N + 1 coefficients c ₀ and c _j with two conditions; 1) the function $V\left(\overrightarrow{r}\right)$ should be equal to the measured potential V _i , and 2) the sum of all c _j coefficients should be zero. In the matrix solution of the N + 1, equations of c ₀ and c _j can be obtained by computing the inverse of the matrix of the concatenated equations.

SSI parameters

In the spherical spline interpolation method, the two parameters m and n have critical effects in the accuracy of the interpolation. The constant m is the “order” of the interpolation. By increasing m, the sum in g _m (x) of (2) converges more rapidly, leading to a smoother interpolation curve (Ferree, 2006). The summation in (2) should be truncated with a certain number of n in practice, and Perrin et al. (1989) stated that only the first seven terms in the summation were adequate to obtain an accuracy of 10⁻⁶ on g(x) with m = 4. Since m = 4 and n = 7 were suggested the original paper of spherical spline interpolation for EEG, these values have been used frequently in EEG studies although some papers suggested to use other values (Ferree, 2006; Soong et al., 1993; Kayser & Tenke, 2006). In addition, Perrin et al. (1989) also introduced another parameter λ, which is small positive value added to the diagonal in the matrix solutions of (1) to regularize. The larger λ, the more spatial smoothing in the interpolation.

Importantly, most EEG studies do not specify the parameter values in spite of their importance. Appropriate values for the SSI parameters can vary across different scalp electrode densities. Additionally, the three aforementioned parameters can interact with each other across varying number of scalp electrodes. However, the appropriate parameter range has not been clearly addressed in the literature. In terms of the m parameter, there is general agreement that m of 3 or 4 works best based on systematic examination (Perrin et al., 1989; Ferree, 2006; Soong et al., 1993), but there is a lack of clear guidelines for parameters n and λ for spherical spline methods. For example, some studies used a very large λ value of 1.0 (Soong et al., 1993), and others used 10⁻⁵ (Kayser & Tenke, 2006). In terms of n, studies evaluating SSI performance (Fletcher et al., 1996, Soong et al., 1993) used 7, while more recent study employing EEG Laplacian transformation using spherical spline method used 50 (Kayser & Tenke, 2006).

In the present study, we used 2 levels of m (3 and 4), 4 levels of n (7, 10, 20, and 50), and 4 levels of λ (10⁻⁵, 10⁻⁶, 10⁻⁷, and 10⁻⁸) for spherical spline interpolation. We firstly examined the accuracy of interpolation by computing Pearson correlation as a similarity measure and root of mean square error (RMSE) as a dissimilarity measure between the recorded and interpolated time-series obtained by the 32 combinations of the m, n, and λ parameters. Every electrode signal was interpolated one at a time using all the other electrode signals. In the preliminary analyses, both the similarity index of Pearson correlation and dissimilarity index RMSE indicated substantial variations of the accuracy of interpolated signals across the number of scalp electrodes (N_e) and spherical spline interpolation parameter m, n, and λ (see Fig. 1). In general, the interpolated signal was more accurate in 64 channel data than in 128 channel data. In addition, m = 3, larger n, and smaller λ had higher accuracy than m = 4, smaller n, and larger λ, respectively (see Fig. 1). Based on these preliminary results, the best and the worst parameter values were chosen for the subsequent analyses of the effect of SSI on interregional PLV; m: 3 and 4, n: 7 and 50, λ: 10⁻⁵and 10⁻⁸.

Fig. 1

Indices of the accuracy of spherical spline interpolation across various conditions. Top panel: similarity index of temporal Pearson correlations between original and interpolated signals, which were averaged across all electrodes. Bottom panel: dissimilarity index of root of mean square error (RMSE) computed using the differences between original and interpolated signals, which were averaged across all electrodes. Both the similarity and dissimilarity indices strongly varied across the levels of number of electrodes, m, n, and λ parameters. All of the between-level differences were statistically significant in paired t-tests (p < .001), even between the n = 20 and n = 50 that look nearly identical due to very little variability in the difference values

Interregional phase coherence

The preprocessed EEG data were further down-sampled to a 128 Hz sampling rate and re-segmented into two second epochs, then time-varying complex energy time-frequency distribution (TFD) of only the eyes-closed resting EEG data were obtained using reduced interference distribution (RID) Rihaczek method (Aviyente et al. 2011). The RID Rihaczek method computes a complex TFD with uniform time-frequency resolution without the trade-off between time and frequency resolution inherent to wavelet analysis. Eleven scalp electrodes (FPz, Fz, F5, F6, Cz, C5, C6, Pz, P5, P6, and Oz in the international 10–10 electrode system)¹ covering most scalp area were chosen for the interpolation. They were used as reference-electrodes whose interregional phase locking values (PLVs; Lachaux et al., 1999) were computed as the consistency of TFD phase differences with all other pair-electrodes across sweeps.

In the resting EEG data where no specific event time existed, time information was meaningless. Thus, PLV-TFD data were averaged across time points, then alpha (8 ~ 13 Hz) and gamma (31 ~ 64 Hz) frequency PLV was computed as the mean within the frequency bands. The difference in interregional PLVs (dPLV) between interpolated signals and the original signals were quantified by simple subtraction (dPLV = PLV_interpolated - PLV_original). As shown in Fig. 2, the dPLV of the reference electrodes were observed not only in their near electrodes, but also in distant electrodes. Therefore, the dPLV for the alpha and gamma frequency bands were further grouped into near and distant pair electrodes.

Fig. 2

Topographies of grand average interregional phase locking values (PLV) for 11 reference electrodes. The grand average PLV of the reference electrodes with all the other electrodes were obtained from 11 subjects in 64 and 128 channel resting EEG data. The PLVs were computed using original recorded reference electrode signals (oPLV) and the interpolated signal computed by SSI (iPLV). For illustration purposes, only alpha (8 ~ 13 Hz) band oPLV and iPLV were computed using interpolated signals with SSI parameter m = 4, n = 7, and λ = 10⁻⁵ in this figure. The differences between iPLV and oPLV were computed by simple subtraction (dPLV = oPLV - iPLV). As shown in dPLV topographies, there were clear differences not only in neighboring pair-electrodes close to the reference-electrodes but also in distant pair-electrodes across all of the 11 reference-electrodes. dPLVs were averaged across pair-electrodes in near and distant groups for the locations marked with red dots in the topographies

Statistical analysis

We conducted a mixed repeated-meas. multivariate analysis of variance (MANOVA) with a dependent variable of the PLV differences (dPLV), one between-group factor of number of electrodes (N_e), and six within-group factors including frequency (alpha and gamma), m (3 and 4), n (7 and 50), λ (10⁻⁵and 10⁻⁸), distance of pair electrode groups (DIST_pair: near and distant), and the locations of the reference electrodes (LOC_ref: FPz, Fz, F5, F6, Cz, C5, C6, Pz, P5, P6, and Oz). With this complex repeated-measure MANOVA, the main effects of the 7 independent variables were examined. Since our main interest was to examine interregional PLVs between distant scalp electrodes that are less vulnerable to volume conduction effects, we conducted a subsequent repeated-measure MANOVA using only dPLV of distant electrode group as a dependent variable. In this analysis, any other significant factors of the first repeated-measure MANOVA were included as independent variables in addition to the 3 SSI parameters to examine interactions among the factors. To estimate the magnitude of effects, partial η ² was reported as effect size for significant main and interaction effects. Significant interactions were further examined with follow-up tests.

Results

Main effects of number of electrodes, spherical spline parameters, and electrode locations on EEG interpolation error

Except frequency, all the other independent variables had significant main effects (N_e: F[1,20] = 32.85, p < .001, η ² = .62; m: F[1,20] = 224.79, p < .001, η ² = .92; n: F[1,20] = 310.20, p < .001, η ² = .94, λ: F[1,20] = 259.90, p < .001, η ² = .93, electrode-group: F[1,20] = 174.82, p < .001, η ² = .90, location of reference electrodes: F[7.9,200] = 3.37, p = .001, η ² = .14). As shown in Fig. 3, 128 channel data had larger dPLV than 64 channel. Consistent with the preliminary results of temporal correlation and RMSE between the original and the interpolated signals with various SSI parameters, m = 3, n = 50, and λ = 10⁻⁸ had lower dPLV than m = 4, n = 50, and λ = 10⁻⁸, respectively. Electrodes distant from the reference electrodes had lower dPLV than those close to the reference electrodes. Lastly, dPLV significantly varied across the locations of the interpolated reference electrodes.

Fig. 3

Main effects of the seven factors on SSI-induced distortions in interregional phase coherence measured by PLV difference between original and interpolated data (dPLV). As shown in the bar graphs, SSI with a larger number of electrodes had larger dPLV than those with a smaller number of electrodes, but dPLV was not different between low (alpha) and high (gamma) frequencies in general. All of the three SSI parameters had a strong effect in dPLV, where SSIs with m = 4, smaller n, and larger λ resulting in larger dPLVs than those with the other parameter values. Near pair-electrodes had larger dPLV than distant pair-electrodes. SSI-induced dPLV also significantly varied across the locations of the interpolated electrodes

Interaction effects of frequency with distance of pair-electrodes and SSI parameters on EEG interpolation error

Although the main effect of frequency was not significant, its interaction with distance of pair-electrodes and SSI parameters m, n, and λ were significant (Frequency × DIST_pair: F[1,20] = 34.75, p < .001, η ² = .64; Frequency × m × DIST_pair: F[1,20] = 33.71, p < .001, η ² = .66; Frequency × n × DIST_pair: F[1,20] = 27.08, p < .001, η ² = .58; Frequency × λ × DIST_pair: F[1,20] = 38.23, p < .001, η ² = .66). As depicted in Fig. 4a), dPLV of gamma band was significantly larger than that of alpha band only in the near pair-electrodes. In contrast, alpha band dPLVs were slightly larger than that of gamma band in distant pair-electrodes, although, the difference was not statistically significant. The significant three-way interaction effects involving each of the three SSI parameters in addition to the frequency by DIST_pair interaction further demonstrated that the larger gamma band dPLV in near pair-electrodes were observed, especially when SSI parameters that caused larger interpolation errors (i.e., m = 4, n = 7, and λ = 1e-5) were used. In other words, the SSI-induced large PLV distortions in high frequency band compared to low frequency band tend to occur in near pair-electrodes depending on the selected SSI parameters.

Fig. 4

Interaction of frequency with distance between electrode pairs (DIST_pair) and SSI parameters. a Significant frequency by DIST_pair interaction indicated that the frequency effect (larger dPLV in gamma than in alpha band) was significant only in near electrode pairs. Significant 3-way interactions among frequency, DIST_pair, and the 3 SSI-parameters are depicted with bar graphs of frequency by SSI-parameter interactions in near (left panels: b, d, f) and distant electrode pairs (right panels: c, e, g). The frequency effects (i.e., larger dPLVs in gamma than in alpha band) in near electrode pairs were observed especially with SSI parameters that led to larger interpolation errors (i.e., m = 4, n = 4, and λ = 1e-5)

Interaction effects of location of reference-electrodes with distance of pair-electrodes and SSI parameters on interpolation error

In addition to the significant main effect, location of reference-electrodes (LOC_ref) had significant interactions with DIST_pair (F[8.9,177.6] = 3.35, p = .001, η ² = .14), indicating that the effect of LOC_ref was significant only in the distant pair-electrodes. Furthermore, the LOC_ref by DIST_pair interaction varied across levels of number of electrodes (LOC_ref × DIST_pair × Ne: F[8.9,177.6] = 3.60, p < .001, η ² = .15) and SSI parameters (LOC_ref × DIST_pair × m: F[7.9,158.3] = 3.82, p < .001, η ² = .16; LOC_ref × DIST_pair × λ: F[8.0,159.6] = 3.49, p = .001, η ² = .15). As shown in Fig. 5, these three-way interactions indicated that the effect of LOC_ref was even larger in distant pair-electrodes with 128 channel data and some SSI parameters (i.e., m = 4 and λ = 10⁻⁵) than with 64 channel data and alternative SSI parameters (i.e., m = 3 and λ = 10⁻⁸).

Fig. 5

Effect of the scalp location of the interpolated signal (i.e., reference-electrodes; LOC_ref) and its interaction with distance between electrode pairs (DIST_pair), number of electrodes (N_e), and SSI parameters m and λ. The p-values in the titles or legends of each graph stands for those of repeated measure ANOVAs testing significance of LOC_ref effect in the corresponding condition. a Error-bars of PLV difference (dPLV) across 11 locations of reference-electrodes in near and distant pair-electrodes, depicting a significant LOC_ref by DIST_pair interaction; dPLV was significantly different across the 11 reference electrodes in distant pair-electrodes, but not in near pair-electrodes. Two-way interactions of LOC_ref by N_e (b and c), LOC_ref by m (d and e), and LOC_ref by λ (f and g) depicted for near (left panels) and distant pair-electrodes (right panels) illustrate the three-way interactions among them. Significant LOC_ref effect in distant pair-electrodes were even larger with 128 channel electrodes, m = 4, and λ = 1e-8, which caused large interpolation errors. LOC_ref effect was even significant in near pair-electrodes when SSI was conducted with these levels of N_e, m and λ

Specific analyses of long-range phase synchrony error from interpolation of EEG data

Interregional PLV of scalp electrodes can be useful when they are measured between distant electrodes, thus less contaminated by volume conduction effects (Srinivasan et al., 2007). Therefore, we conducted a subsequent repeated-measure MANOVA with only dPLVs of the distant pair-electrode groups in order to examine effects of the SSI parameters, number of electrodes, and their interactions in PLV measures of long distance interregional functional connectivity. Given that frequency effect was not significant in dPLVs of the distant pair-electrodes, this analysis included only dPLVs of the alpha band which is the dominant frequency in the eyes-closed resting state. However, the location of the interpolated reference-electrodes was included in this analysis because this variable was found to be more significant in the distant than in near pair-electrodes. Therefore, in the second repeated measure MANOVA for dPLV of the distant pair-electrodes had five independent variables, including number of electrodes (N_e), the three SSI parameters, and location of reference-electrodes (LOC_ref).

The effects of SSI parameters

The three-way interaction of m × n × λ was significant (F[1,20] = 118.59, p < .001, η ²  = .86), in addition to the two-way interactions between them (m × n: F[1,20] = 215.96, p < .001, η ² = .92; m × λ: F[1,20] = 168.04, p < .001, η ² = .89; n × λ: F[1,20] = 237.51, p < .001, η ² = .92), indicating that an effect of one of the SSI parameters is dependent on the level of the other two SSI parameters. In other words, dPLV significantly varied across the eight combinations of the SSI parameters. Moreover, the three way interaction had significant interaction with number of electrodes (m × n × λ × N_e: F[1,20] = 10.38, p = .004, η ² = .34), suggesting that the dPLV variation across the eight conditions of SSI parameters was larger in 128 channel data than in 64 channel data (see Figs. 6 and 7).

Fig. 6

Interactions among the SSI parameters in 64 (top row) and 128 channel data (bottom row). Three-way interactions of the three SSI parameters were depicted for 64 (top row) and 128 channel (bottom row) data. Significant three-way interactions were observed in both 64 and 128 channel data, in which effect size of a parameter depended on the levels of other two parameters. For example, m = 3 had generally lower dPLV than m = 4. But when n = 7 and λ = 1e-5 were used for SSI, SSI-induced dPLV with m = 3 had larger dPLV than those of m = 4 in some sets of n and λ (i.e., n = 7 and λ = 1e-8 or with n = 50 and λ = 1e-8). The pattern of three-way interaction was somewhat different between 64 and 128 channel data as suggested by the significant four-way interaction among m, n, λ, and N_e

Fig. 7

Effect of location of reference-electrodes on dPLV across eight sets of SSI parameters in 64 (top) and 128 channel data (bottom). The p-values of repeated-measure ANOVAs testing significance of LOC_ref effect as well as the effect size (η ²) of LOC_ref are shown in each error-bar plot of the corresponding conditions. The eight SSI parameter conditions were sorted so that the amount of dPLV decreases from right to left. Note that the order of SSI sets were the same between 64 and 128 channel data although the patterns of dPLV variation across the 11 reference-electrodes are different between them. Importantly, PLV inflation by SSI was close to zero in both 64 and 128 channel data when the best set of SSI parameters was used (i.e., m = 3, n = 50, and λ = 1e-8)

The effects of location of interpolated electrodes

The five-way interaction of N_e × m × n × λ × LOC_ref was also significant (F[7.8,156.2] = 3.85, p < .001, η ² = .16). By taking the three-way interaction of m × n × λ as a single factor of SSI parameter (P_SSI) with eight levels, this complex interaction can be understood as a three-way interaction of N_e × P_SSI × LOC_ref . That is, the effect size of LOC_ref (i.e., amount of dPLV variation across the 11 reference-electrodes) varied across the eight P_SSI conditions, and the effect size of LOC_ref was larger in 128 channel than in 64 channel. Moreover, the pattern of dPLV variation across the 11 reference-electrodes were somewhat different between 64 and 128 channels. However, as shown in Fig. 7, the variation of dPLV across the reference-electrode locations as well as the amount of dPLV were minimal with an optimally chosen set of SSI parameters (i.e., m = 3, n = 50, and λ = 1e-8) in both 64 and 128 channel data.

Discussion

In this systematic examination of the effect of interpolated data on EEG-based indices of inter-regional functional connectivity, we found a substantial artifactual increase of phase coherence values with the more “standard” choice of EEG interpolation parameters. In general, overestimates of phase coherence were more evident in interpolations involving a larger numbers of electrodes (i.e., greater for 128 channel than 64 channel EEG data) and for near pair-electrodes than for distant pair-electrodes. Consistent with the preliminary results in the evaluations of SSI performances using temporal correlation and RMSE measures, different levels of the three SSI parameters resulted in large differences in the interchannel PLV inflation. Particularly, m of 3, larger n (i.e., 50), and smaller λ (i.e., 1e-8) had better SSI results than m of 4, smaller n (i.e., 7), and larger λ (i.e., 1e-5), respectively. As indicated by the significant three-way interaction between the three SSI parameters, the effects of a single SSI parameter depended on the levels of the other SSI parameters. That is, in order to optimize SSI accuracy, it is required to consider all three parameters and their interactions rather than one or two parameters at once. SSI-induced PLV inflation also varied across the location of interpolated reference-electrodes, especially when inappropriate SSI parameters were used. Additionally, PLV inflation varied based on the frequency of the EEG signal of interest for inappropriate SSI parameters in near pair-electrodes. Overestimates of PLVs for near pair-electrodes were larger in gamma than in alpha frequencies. Most importantly, it was found that SSI was able to produce interpolated signals nearly identical to original signals with the best set of parameters (i.e., m = 3, n = 50, and λ = 1e-8) for both 64 and 128 channel EEG data.

Effects of the number of channels on interpolation accuracy

In order to have accurate topographic mapping through interpolation techniques, a large number of electrodes are required. In comparisons of several interpolation methods on global and local accuracy for reproducing original values across many spatial reference set points, a simulation study indicated that a greater number of channels resulted in more accurate spatial interpolation (Fletcher et al., 1996). Our finding of the inverse relationship between SSI accuracy and number of channels should be considered with respect to the idea that interpolation basically estimates a missing signal by weighted summation of other electrode signals. More channels should result in more common signals in the weighted sum, and thus interpolation of the EEG time-series based on more cosine angles between many electrode locations could artifactually increase the interregional coherence measures estimating phase consistency between electrodes.

Effect of SSI parameters

Although we found that performance of SSI depended on the selection of the three SSI parameters, EEG studies using SSI usually do not report the parameter values. Given that studies typically referred to early studies of SSI literature for EEG (Perrin et al., 1989), most of the published work to date likely uses m = 4 and n = 7 which Perrin and colleagues recommended. Early SSI studies provide no specific recommended value for the regularization parameter λ, thus it is hard to know what values have been used. Few studies have systematically examined SSI performance with various m parameter values. In an examination of temporal accuracy of interpolation methods, it was found that smaller m values produced more accurate interpolated time-series (Soong et al., 1993). Actually, m = 2 was the best in terms of the temporal similarity, but the interpolated time-series with m = 2 increased biased amplitudes of the time-series, thus m = 3 was recommended. More recent work confirmed that m = 3 yielded acceptable SSI performance (Ferree, 2006).

In terms of the Legendre polynomial order parameter n, we failed to find any recommendations based on objective and systematic evaluations in the literature. Only open source or commercial EEG analysis programs have recommendations for the choice of n, which are not 7. For example, the spherical spline interpolation toolbox plugin for EEGlab has a default n of 50 (http://www.uni-leipzig.de/~biocog/content/widmann/eeglab-plugins/, retrieved on March 6, 2015), while Brain Vision Analyzer defaults n to 10 (http://erpinfo.org/boot-camp-2007-materials/2011-boot-camp-materials/boot-camp-tutorials/documentation/Analyzer.pdf. retrieved on March 6, 2015). The order of Legendre polynomial defines spatial harmonic frequencies with respect to each electrode and determines spatial frequency precision of the SSI results (Cohen, 2014); a higher order parameter will increase the accuracy of SSI. Our finding of a significant n × N_e interaction (F[1,20] = 40.78, p < .001, η ²  = .67) where smaller n caused larger PLV inflation in 128 than 64 channel EEG data is consistent with the higher order polynomial being required for a larger number of channels. Although it is hard to find a clear guideline in published studies, results from the present study indicate that n = 50 yields sufficiently accurate results for both 64 and 128 channel data.

Lastly, the regularization or spatial smoothing parameter λ also had critical effects on PLV inflation through SSI. As shown in Fig. 3, different levels of λ had the largest dPLV differences compared those observed between levels of other variables. Although λ = 1e-5 is a very small value, this level of regularization parameter led to substantial increases in PLV. In our preliminary analyses, we even tried to avoid regularization of the SSI solution (i.e., λ = 0) which led to slightly higher temporal correlations and smaller RMSE between original and interpolated signals than λ = 1e-8 in several cases. However, no regularization sometimes generated a singular matrix of c _j equations that were not invertible in the SSI solution, causing very inaccurate results. Therefore, very small (λ = 1e-8) but non-zero λ is recommended to ensure temporal accuracy of interpolated time-series.

Location of interpolated electrodes

We found the amount of PLV inflation varied across the location of interpolated reference-electrodes. The patterns of PLV inflation across the electrode sites were different between 64 and 128 channels (see Figs. 5 and 7). This might be related to the spatial features of neural sources such as depth and eccentricity (Perrin et al., 1987; Fletcher et al. 1996) as well as head shapes, which should vary across individuals. Both the 64 and 128 channel EEG data were collected during resting states, thus their interregional PLV patterns were quite similar but not identical, especially in the locations of distant pair-electrodes that showed the greatest PLV inflation by SSI (see Fig. 2). Many factors may contribute to their dissimilarities including slight differences in the location of the reference electrodes (i.e., electrodes chosen to approximate the locations of F5, F6, P5, and P6 in 128 channel data), a larger number of electrode pairs covering extended scalp areas in 128 channels, older ages of participants whose EEG data were collected with 64 channels, and so on. In general, interpolated signals at F5, F6, and Pz electrodes had the largest PLV distortion in distant pair electrodes (see Fig. 5). Interestingly, the actual locations of these electrodes may deviate most from what is assumed in a spherical head model. As Law et al. (1993) pointed out, head shape is not exactly spherical, and thus a three-dimensional spline interpolation that accepts realistic head models might perform better with less variation in interpolation accuracy. Notably, we found that with the optimal set of SSI parameters, variations of PLV inflation across electrode locations were effectively removed.

Conclusions

EEG collected from dense-array electrode systems have greatly improved the spatial resolution of measures used to characterize brain dynamics, allowing more precise localization of sources of scalp EEG features with adequate spatial sampling. Along with this advance comes increased probability of electrode-specific artifacts (i.e., artifactual electrode noise or electrolyte bridges). A missing electrode leads to data loss in multiple electrode-pair variables in interelectrode coherence studies. Furthermore, advanced multivariate analysis techniques such as PCA and ICA require no missing channels across subjects, thus interpolation to replace artifactual electrode signals is needed to prevent significant loss of data for individual trials or subjects. Although it has been believed that SSI is acceptably accurate in dense-array electrodes (>64 channels; Junghöfer et al., 2000; Greischar et al., 2004), the findings of the present study demonstrate that interpolation can only be considered accurate when performed after careful selection of interpolation parameters (m, n, and λ).

Specifically, we found that conventional parameter recommendations for interpolating missing data for an EEG channel failed to produce temporally accurate EEG time-series, causing substantial errors in interregional functional connectivity measures. Such interpolation errors appear to be due to two main reasons. First, interpolation methods were originally developed for topographic mapping of scalp EEG activities, thus most studies evaluated interpolation performances in spatial accuracy of similarity between true (simulated) and interpolated values across many spatial points (Perrin et al., 1987, 1989; Soufflet et al., 1991; Fletcher et al., 1996) rather than temporal accuracy of similarity between true and interpolated time-series. Given that SSI using a larger spatial smoothing λ parameter caused more distortions in the temporal features of interpolated signals, employing SSI to replace an artifactual electrode time-series requires other parameter settings. Second, the recommendations for SSI parameters in an original article (Perrin et al., 1989) were made based on results for a much smaller number of electrodes (<32 channels) than dense-array EEG recordings that are used today. As indicated by the significant interactions between SSI parameters and number of channels, EEG recordings with a larger number of channels requires different parameter settings from recordings with small numbers of channels (i.e., larger n and smaller λ). Since our investigation included only 64 and 128 channel data, dense-array EEG data with even larger numbers of channels might require further adjustments to SSI parameter values beyond our recommendation. However, given the relationships between SSI parameters and temporal accuracy of SSI revealed in the current investigation, the parameter settings of m = 3, a large n around 50, and a small λ around 1e-8 or lower may apply for EEG recordings with an even higher density of electrodes. The present study evaluated SSI performance specifically for the accuracy of interregional PLV measures, but these findings are likely to apply to other EEG measures of time-series data.