 Methodology article
 Open Access
 Published:
Distortions in EEG interregional phase synchrony by spherical spline interpolation: causes and remedies
Neuropsychiatric Electrophysiology volume 1, Article number: 9 (2015)
Abstract
Background
Measures of the coherence of electroencephalography (EEG) timeseries recorded at spatially distant points on the scalp are often used by researchers to characterize the dynamic interactions of brain regions. In densearray 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 distortionfree 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 timefrequency crossspectra between two timeseries, 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 bandpass 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 rereferencing. In a simulation study that included coherence values between 111 electrodes across various reference electrodes, coherence data using an average reference closely aligned with referenceindependent 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 denseelectrode 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, densearray 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), thinplate spline interpolation (Harder & Desmarais, 1972), spherical spline interpolation (SSI; Perrin et al., 1989), and 3D 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 timeseries 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 crossvalidation 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. Twentytwo 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 45s blocks alternating conditions with eyes open or closed. Similarly, the 128 channel resting session alternated conditions, but included six 60s 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 equidistantlayout. EEG data were recorded with a sampling rate of 1024 Hz, and downsampled to 256 Hz in offline processing. In offline processing, EEG data were referenced to a linked earlobe reference signal and highpass 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 timeseries of ICs (McMenamin et al., 2010). Denoised data were reconstructed by removing noise IC variances from the EEG data, then rereferenced 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
where c _{ 0 } 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
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:
with P _{ 0 } = 1 and P _{ 1 } = 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 _{ 0 } 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 _{ 0 } 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^{−6} 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^{−5} (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^{−5}, 10^{−6}, 10^{−7}, and 10^{−8}) 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 timeseries 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^{−5}and 10^{−8}.
Interregional phase coherence
The preprocessed EEG data were further downsampled to a 128 Hz sampling rate and resegmented into two second epochs, then timevarying complex energy timefrequency distribution (TFD) of only the eyesclosed 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 timefrequency resolution without the tradeoff 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)^{Footnote 1} covering most scalp area were chosen for the interpolation. They were used as referenceelectrodes whose interregional phase locking values (PLVs; Lachaux et al., 1999) were computed as the consistency of TFD phase differences with all other pairelectrodes across sweeps.
In the resting EEG data where no specific event time existed, time information was meaningless. Thus, PLVTFD 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.
Statistical analysis
We conducted a mixed repeatedmeas. multivariate analysis of variance (MANOVA) with a dependent variable of the PLV differences (dPLV), one betweengroup factor of number of electrodes (N_{e}), and six withingroup factors including frequency (alpha and gamma), m (3 and 4), n (7 and 50), λ (10^{−5}and 10^{−8}), 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 repeatedmeasure 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 repeatedmeasure MANOVA using only dPLV of distant electrode group as a dependent variable. In this analysis, any other significant factors of the first repeatedmeasure 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 η ^{2} was reported as effect size for significant main and interaction effects. Significant interactions were further examined with followup 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, η ^{2} = .62; m: F[1,20] = 224.79, p < .001, η ^{2} = .92; n: F[1,20] = 310.20, p < .001, η ^{2} = .94, λ: F[1,20] = 259.90, p < .001, η ^{2} = .93, electrodegroup: F[1,20] = 174.82, p < .001, η ^{2} = .90, location of reference electrodes: F[7.9,200] = 3.37, p = .001, η ^{2} = .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^{−8} had lower dPLV than m = 4, n = 50, and λ = 10^{−8}, 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.
Interaction effects of frequency with distance of pairelectrodes and SSI parameters on EEG interpolation error
Although the main effect of frequency was not significant, its interaction with distance of pairelectrodes and SSI parameters m, n, and λ were significant (Frequency × DIST_{pair}: F[1,20] = 34.75, p < .001, η ^{2} = .64; Frequency × m × DIST_{pair}: F[1,20] = 33.71, p < .001, η ^{2} = .66; Frequency × n × DIST_{pair}: F[1,20] = 27.08, p < .001, η ^{2} = .58; Frequency × λ × DIST_{pair}: F[1,20] = 38.23, p < .001, η ^{2} = .66). As depicted in Fig. 4a), dPLV of gamma band was significantly larger than that of alpha band only in the near pairelectrodes. In contrast, alpha band dPLVs were slightly larger than that of gamma band in distant pairelectrodes, although, the difference was not statistically significant. The significant threeway 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 pairelectrodes were observed, especially when SSI parameters that caused larger interpolation errors (i.e., m = 4, n = 7, and λ = 1e5) were used. In other words, the SSIinduced large PLV distortions in high frequency band compared to low frequency band tend to occur in near pairelectrodes depending on the selected SSI parameters.
Interaction effects of location of referenceelectrodes with distance of pairelectrodes and SSI parameters on interpolation error
In addition to the significant main effect, location of referenceelectrodes (LOC_{ref}) had significant interactions with DIST_{pair} (F[8.9,177.6] = 3.35, p = .001, η ^{2} = .14), indicating that the effect of LOC_{ref} was significant only in the distant pairelectrodes. 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, η ^{2} = .15) and SSI parameters (LOC_{ref} × DIST_{pair} × m: F[7.9,158.3] = 3.82, p < .001, η ^{2} = .16; LOC_{ref} × DIST_{pair} × λ: F[8.0,159.6] = 3.49, p = .001, η ^{2} = .15). As shown in Fig. 5, these threeway interactions indicated that the effect of LOC_{ref} was even larger in distant pairelectrodes with 128 channel data and some SSI parameters (i.e., m = 4 and λ = 10^{−5}) than with 64 channel data and alternative SSI parameters (i.e., m = 3 and λ = 10^{−8}).
Specific analyses of longrange 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 repeatedmeasure MANOVA with only dPLVs of the distant pairelectrode 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 pairelectrodes, this analysis included only dPLVs of the alpha band which is the dominant frequency in the eyesclosed resting state. However, the location of the interpolated referenceelectrodes was included in this analysis because this variable was found to be more significant in the distant than in near pairelectrodes. Therefore, in the second repeated measure MANOVA for dPLV of the distant pairelectrodes had five independent variables, including number of electrodes (N_{e}), the three SSI parameters, and location of referenceelectrodes (LOC_{ref}).
The effects of SSI parameters
The threeway interaction of m × n × λ was significant (F[1,20] = 118.59, p < .001, η ^{2} = .86), in addition to the twoway interactions between them (m × n: F[1,20] = 215.96, p < .001, η ^{2} = .92; m × λ: F[1,20] = 168.04, p < .001, η ^{2} = .89; n × λ: F[1,20] = 237.51, p < .001, η ^{2} = .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, η ^{2} = .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).
The effects of location of interpolated electrodes
The fiveway interaction of N_{e} × m × n × λ × LOC_{ref} was also significant (F[7.8,156.2] = 3.85, p < .001, η ^{2} = .16). By taking the threeway interaction of m × n × λ as a single factor of SSI parameter (P_{SSI}) with eight levels, this complex interaction can be understood as a threeway 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 referenceelectrodes) 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 referenceelectrodes were somewhat different between 64 and 128 channels. However, as shown in Fig. 7, the variation of dPLV across the referenceelectrode 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 λ = 1e8) in both 64 and 128 channel data.
Discussion
In this systematic examination of the effect of interpolated data on EEGbased indices of interregional 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 pairelectrodes than for distant pairelectrodes. 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., 1e8) had better SSI results than m of 4, smaller n (i.e., 7), and larger λ (i.e., 1e5), respectively. As indicated by the significant threeway 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. SSIinduced PLV inflation also varied across the location of interpolated referenceelectrodes, 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 pairelectrodes. Overestimates of PLVs for near pairelectrodes 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 λ = 1e8) 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 timeseries 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 timeseries (Soong et al., 1993). Actually, m = 2 was the best in terms of the temporal similarity, but the interpolated timeseries with m = 2 increased biased amplitudes of the timeseries, 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.unileipzig.de/~biocog/content/widmann/eeglabplugins/, retrieved on March 6, 2015), while Brain Vision Analyzer defaults n to 10 (http://erpinfo.org/bootcamp2007materials/2011bootcampmaterials/bootcamptutorials/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, η ^{2} = .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 λ = 1e5 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 λ = 1e8 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 (λ = 1e8) but nonzero λ is recommended to ensure temporal accuracy of interpolated timeseries.
Location of interpolated electrodes
We found the amount of PLV inflation varied across the location of interpolated referenceelectrodes. 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 pairelectrodes 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 threedimensional 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 densearray 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 electrodespecific artifacts (i.e., artifactual electrode noise or electrolyte bridges). A missing electrode leads to data loss in multiple electrodepair 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 densearray 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 timeseries, 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 timeseries. 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 timeseries 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 densearray 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, densearray 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 1e8 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 timeseries data.
Notes
 1.
Since the Biosemi 128 channel electrode cap does not have exact electrode sites for F5, F6, P5, and P6 of the international 10–10 electrode system, the closest frontolateral and parietolateral electrodes were chosen to have matched electrodes for interpolation between the 64 and 128 channel EEG data.
Abbreviations
 DIST_{pair} :

Distance of pairelectrodes
 dPLV:

Difference phase locking value
 ECI:

Electrode cap international
 EEG:

Electroencephalography
 IC:

Independent component
 ICA:

Independent component analysis
 iPLV:

Phase locking value computed using interpolated electrode signals
 LOC_{ref} :

Location of referenceelectrodes
 MANOVA:

Multivariate analysis of variance
 N_{e} :

Number of electrodes
 oPLV:

Phase locking value computed using original recorded electrode signals
 PLV:

Phase locking value
 P_{SSI} :

Single factor of spherical spline interpolation parameters (i.e., P_{SSI} = m × n × λ)
 RID:

Reduced interference distribution
 RMSE:

Root of mean square error
 SSI:

Spherical spline interpolation
 TFD:

Timefrequency distribution
References
Alschuler DM, Tenke CE, Bruder GE, Kayser J (2014) Identifying electrode bridging from electgrode distance distributions: a survey of publiclyavailable EEG data using a new method. Clin Neurophysiol 112(3):545–550. doi:10.1016/j.clinph.2013.08.024
Aviyente S, Bernat EM, Evans WS, Sponheim SR (2011) A phase synchrony measure for quantifying dynamic functional integration in the brain. Hum Brain Mapp 32:80–93. doi:10.1002/hbm.21000
Biggins CA, Fein G, Raz J, Amir A (1991) Artifactually high coherences result from using spherical spline computation of scalp current density. Electroencephalogr Clin Neurophysiol 79:413–419. doi:10.1016/00134694(91)90206J
Chorlian DB, Rangaswamy M, Porjesz B (2009) EEG coherence: topography and frequency structure. Exp Brain Res 198:59–83. doi:10.1007/s0022100919369
Cohen MX (2014) Analyzing neural time series data: theory and practice. MIT Press, Massachusetts
de Boor C (1978) A practical guide to splines. Springer, Heidelberg
Ferree TC (2006) Spherical splines and average referencing in scalp electroencephalography. Brain Topogr 19:43–52. doi:10.1007/s1054800600110
Fletcher EM, Kussmaul CL, Mangun GR (1996) Estimation of interpolation errors in scalp topographic mapping. Electroencephalogr Clin Neurophysiol 98:422–434. doi:10.1016/00134694(96)951354
Greischar LL, Burghy CA, Van Reekum CM, Jackson DC, Pizzagalli DA, Mueller C, Davidson RJ (2004) Effects of electrode density and electrolyte spreading in dense array electroencephalographic recording. Clin Neurophysiol 115(3):710–720. doi:10.1016/j.clinph.2003.10.028
Griffiths DJ (2007) Introduction to electrodynamics, 3rd edn. Hall, Prentice
Harder RL, Desmarais RN (1972) Interpolation using surface splines. J Aircraft 9:189–191. doi:10.2514/3.44330
Hyvärinen A, Oja E (2000) Independent component analysis: algorithms and applications. Neural Netw 13:411–30. doi:10.1016/S08936080(00)000265
Junghöfer M, Elbert T, Tucker DM, Rockstroh B (2000) Statistical control of artifacts in dense array EEG/MEG studies. Psychophysiology 37(4):523–532. doi:10.1111/14698986.3740523
Kayser J, Tenke CE (2006) Principal components analysis of Laplacian waveforms as a generic method for identifying ERP generator patterns: II. Adequacy of lowdensity estimates. Clin Neurophysiol 117:369–380. doi:10.1016/j.clinph.2005.08.033
Lachaux JP, Rodriguez E, Martinerie J, Varela FJ (1999) Measuring phase synchrony in brain signals. Hum Brain Mapp 8:194–208. doi:10.1002/(SICI)10970193(1999)8:4<194::AIDHBM4>3.0.CO;2C
Law SK, Nunez PL, Wijesinghe RS (1993) Highresolution EEG using spline generated surface Laplacians on spherical and ellipsoidal surfaces. IEEE Trans Biomed Eng 40:145–153. doi:10.1109/10.212068
McMenamin BW, Shackman AJ, Maxwell JS, Bachhuber DRW, Koppenhaver AM, Greischar LL et al (2010) Validation of ICAbased myogenic artifact correction for scalp and sourcelocalized EEG. Neuroimage 49:2416–2432. doi:10.1016/j.neuroimage.2009.10.010
Nunez PL, Srinivasan R (2006) Electrical fields of the brain: the neurophysics of EEG, 2nd edn. Oxford Press, New York
Nunez PL, Srinivasan R, Westdorp AF, Wijesinghe RS, Tucker DM, Silberstein RB, Cadusch PJ (1997) EEG coherency I: Statistics, reference electrode, volume conduction, Laplacians, cortical imaging, and interpretation at multiple scales. Electroencephalogr Clin Neurophysiol 103:499–515. doi:10.1016/S00134694(97)000667
Perrin F, Pernier J, Bertrand O, Giard MH, Echallier JF (1987) Mapping of scalp potentials by surface spline interpolation. Electroencephalogr Clin Neurophysiol 66:75–81. doi:10.1016/00134694(87)901416
Perrin F, Pernier J, Bertrand O, Echallier JF (1989) Spherical splines for scalp potential and current density mapping. Electroencephalogr Clin Neurophysiol 72:184–187. doi:10.1016/00134694(89)901806
Roach BJ, Mathalon DH (2008) Eventrelated EEG timefrequency analysis: an overview of measures and an analysis of early gamma band phase locking in schizophrenia. Schizophr Bull 34:907–926. doi:10.1093/schbul/sbn093
Schwarz G (1978) Estimating the dimension of a model. Ann Stat 6:461–464. doi:10.1214/aos/1176344136
Shepard, D. (1968). A twodimensional interpolation function for irregularlyspaced data. In Proceedings of 23rd ACM National Conference (pp. 517–524). doi:10.1145/800186.810616
Soong AC, Lind JC, Shaw GR, Koles ZJ (1993) Systematic comparisons of interpolation techniques in topographic brain mapping. Electroencephalogr Clin Neurophysiol 87:185–195. doi:10.1016/00134694(93)90018Q
Soufflet L, Toussaint M, Luthringer R, Gresser J, Minot R, Macher JP (1991) A statistical evaluation of the main interpolation methods applied to 3dimensional EEG mapping. Electroencephalogr Clin Neurophysiol 79:393–402. doi:10.1016/00134694(91)90204H
Srinivasan R, Winter WR, Ding J, Nunez PL (2007) EEG and MEG coherence: measures of functional connectivity at distinct spatial scales of neocortical dynamics. J Neurosci Methods 166:41–52. doi:10.1016/j.jneumeth.2007.06.026
Tenke CE, Kayser J (2001) A convenient method for detecting electrolyte bridges in multichannel electroencephalogram and eventrelated potential recordings. Clin Neurophysiol 112(3):545–550. doi:10.1016/S13882457(00)005538
Varela F, Lachaux JP, Rodriguez E, Martinerie J (2001) The brainweb: phase synchronization and largescale integration. Nat Rev Neurosci 2:229–239. doi:10.1038/35067550
Acknowledgement
This work was supported by Merit Review grants from the Veterans Health Administration Clinical Science Research and Development Program (1 I01 CX00068301) and Rehabilitation Research and Development Program (1 I01 RX00062201A1) to SRS.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
SK conceived of the study, developed its design, carried out data analyses, performed statistical analysis, and drafted the manuscript. TJL processed EEG data, carried out part of data analyses, and wrote part of the manuscript. SRS participated in the study design and coordination, and helped to draft the manuscript. All authors read and approved the final manuscript.
Rights and permissions
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Functional Connectivity
 Independent Component Analysis
 Channel Data
 Phase Coherence
 Scalp Electrode