Unified Bayesian Estimator of EEG Reference at Infinity: rREST (Regularized Reference Electrode Standardization Technique)

The choice of reference for the electroencephalogram (EEG) is a long-lasting unsolved issue resulting in inconsistent usages and endless debates. Currently, both the average reference (AR) and the reference electrode standardization technique (REST) are two primary, apparently irreconcilable contenders. We propose a theoretical framework to resolve this reference issue by formulating both (a) estimation of potentials at infinity, and (b) determination of the reference, as a unified Bayesian linear inverse problem, which can be solved by maximum a posterior estimation. We find that AR and REST are very particular cases of this unified framework: AR results from biophysically non-informative prior; while REST utilizes the prior based on the EEG generative model. To allow for simultaneous denoising and reference estimation, we develop the regularized versions of AR and REST, named rAR and rREST, respectively. Both depend on a regularization parameter that is the noise to signal variance ratio. Traditional and new estimators are evaluated with this framework, by both simulations and analysis of real resting EEGs. Toward this end, we leverage the MRI and EEG data from 89 subjects which participated in the Cuban Human Brain Mapping Project. Generated artificial EEGs—with a known ground truth, show that relative error in estimating the EEG potentials at infinity is lowest for rREST. It also reveals that realistic volume conductor models improve the performances of REST and rREST. Importantly, for practical applications, it is shown that an average lead field gives the results comparable to the individual lead field. Finally, it is shown that the selection of the regularization parameter with Generalized Cross-Validation (GCV) is close to the “oracle” choice based on the ground truth. When evaluated with the real 89 resting state EEGs, rREST consistently yields the lowest GCV. This study provides a novel perspective to the EEG reference problem by means of a unified inverse solution framework. It may allow additional principled theoretical formulations and numerical evaluation of performance.

Introduction

The human electroencephalogram (EEG) has been an indispensable technology for both cognitive and clinical neuroscience for almost 90 years. Ultrahigh temporal resolution, low cost, and non-invasiveness single it out as a translational tool of choice to study the brain. Nevertheless, two main drawbacks of EEG detract from its ability to localize the brain activity: (i) spatial blurring due to volume conduction; (ii) the inherent indeterminacy of potentials measurements which are always carried out with respect to a given reference (Teplan, 2002). Spatial blurring is being addressed by advanced source imaging techniques that however will not be the focus of our attention. We will rather concentrate on the vexing and still incompletely resolved “EEG reference problem.” To precisely define this issue, we note that it is due to the intrinsic nature of EEG recordings that are the measurement of potential differences between two sites shown in Figure 1. Ideally, one would like to record the potentials of an “active electrode” that is only picking up the activities due to a few brain structures in comparison to a neutral “reference electrode” with zero activity. One might think that such a reference electrode could be placed at infinity, yielding the ideal potentials φ_∞. However, this would not work in practice, since this configuration would serve as an antenna, picking unwanted activity from the environment. Some researchers therefore experimented with reference electrode placed on the body so that EEG differential amplimers could eliminate environmental noise with high common mode rejection ratio. Unfortunately, because there is no neutral or inactive point upon the body, these proposals are also inadequate. A physical neutral reference seems therefore to be out of our reach.

FIGURE 1

Figure 1. EEG reference problem. EEG recordings measure the potential differences between active electrodes marked in yellow (only 3 for illustration) and reference marked in red which may be (A) a point at infinity, (B) an electrode placed on the body, e.g., the base of neck, (C) the cephalic reference (here, Cz), and (D) average reference: the mean of potentials over all active electrodes. All reference techniques (B–D) attempt to approximate infinity reference in (A). The issue is that all the proposals results in different EEG waveforms originating the “EEG reference problem.” It should be noted that with the identical source activities j, the different waveforms (φ_∞, v_BR, v_Cz,v_AR) could be taken as well as the outcome of different lead fields (K_∞, K_BR, K_Cz,K_AR).

However, the non-neutrality of the reference has consequences cascading through the following processing stages, including the final statistical result. In view of the failures of physical references, attention was turned to the construction of “virtual” estimators of the neutral references, namely, virtual estimators of φ_∞.

One popular virtual estimator is the “average reference” (AR, Figure 1D), which based on the following logic: (i) the integral of the electrical potential over a sphere, due to a current source inside it, is zero (Goldman, 1950; Offner, 1950); (ii) the head can be approximated as a sphere; (iii) therefore, a neutral reference may be obtained by summing or averaging the activities of all electrodes. Re-referencing proceeds by subtracting this average from all channels. Unfortunately, recent work (Yao, 2017) has shaken the theoretical foundation of AR: the potential integral for a realistic head surface is not zero.

A more biophysics-based virtual estimator of φ_∞ can be obtained by the reference electrode standardization technique (REST, Figure 2) which directly estimates the ideal potentials referenced to a point at infinity (Yao, 2001). REST uses a head model and equivalent sources to localize source activities, then project the source activities to electrodes—now with reference to infinity. Early work on REST was based upon a simple spherical head model. It was soon shown that EEG power maps (Yao et al., 2005), ERP peak values and latencies(Li and Yao, 2007) did, in fact, critically depend on the choice of reference. In a further study, it was shown (Tian and Yao, 2013) that scalp statistical parametric mapping with REST for audiovisual stimulus evoked potentials provided closer correspondence to the source localization by low resolution electromagnetic tomography than that with AR. These encouraging results about REST have been bolstered by several simulation experiments. Using a spherical head model for simulation, Marzetti et al. (2007) and Qin et al. (2010) indicated that REST led to better estimates of EEG spectra and coherence than AR. Several papers unsurprisingly showed that realistic head model for REST gives superior results for the reconstruction of simulated EEG scalp topographies (Liu et al., 2015), functional connectivity (Chella et al., 2016), and bispectral analysis (Chella et al., 2017).

FIGURE 2

Figure 2. Diagram of REST. (A) measured EEG potentials v_Cz with reference to Cz; j and K_Cz are the actual source activities and the actual lead field with reference to Cz, respectively; (B) for REST, one first builds the lead field ${\widehat{\text{K}}}_{\infty }$ with estimated head model and equivalent sources, then transforms the lead field with reference to Cz as ${\widehat{\text{K}}}_{\text{Cz}}$; with this lead field, one enables to estimate the equivalent source activities $\widehat{\text{j}}$; after, the equivalent source activities are taken the forward calculation through ${\widehat{\text{K}}}_{\infty }$ to ${\widehat{\phi }}_{\infty }$ which is the approximation of EEG potentials at infinity.

In spite of these suggestive results in favor of REST, there is still an intense and to a certain extent unresolved debate on which reference is preferable (Nunez, 2010; Kulaichev, 2016). The lack of resolution is due to that simulation studies, while useful, are not enough to demonstrate the superiority of one reference technique over another. In addition to simulations, the evidence is needed on which reference achieves the “best fit” to actual data. The choice of the best model is a well-studied problem in modern statistics (Robert, 2007) and can be resolved by model selection criteria that approximates the Bayesian model evidence (Konishi and Kitagawa, 2008). However, to apply these techniques, an explicit Bayesian model of the “EEG reference problem” must be stated. Thus, one of the primal goals for this study is to uncover a unified estimator of EEG reference at infinity.

In the current study, we formulate, to our best knowledge for the first time, the “EEG reference problem” as a generalized Bayesian inverse problem. One surprising consequence of this approach is the insight that AR and REST share the same model and just differ in the prior distribution for the covariance of EEG potentials at infinity. On the one hand, assuming uncorrelated activities over electrodes leads to the AR estimator. On the other hand, if the correlations between electrodes are assumed to be caused by sources filtered through a volume conductor model, the resulting estimator is REST.

Our theoretical formulation will allow us to examine different reference estimators within a common statistical framework. We note that the REST estimator (Yao, 2001) was originally defined for the case of exact noise-free. In some situations, this might be unrealistic, since the scalp EEG may have quite low signal to noise ratio (Ferree et al., 2001; Lemm et al., 2006; Guruvareddy, 2013; Bigdely-Shamlo et al., 2015). Our framework allows using regularization technique as a way to accommodate noise in the data (Phillips et al., 2002, 2005). The regularized version of REST is developed which we call “rREST.” It is evident that a regularized version of AR is also possible, which we term as “rAR.” AR and REST are just the special cases of rAR and rREST when regularization parameter tends to zero, respectively.

We further investigate the effect of the volume conduction model on rREST. To avoid the “inverse crime” (Kaipio and Somersalo, 2007), the volume conduction model used in the simulation should be different with the one used to generate EEG potentials. We call this “volume conduction model matching problem” for REST that may produce the spurious results in simulation. Although equivalent source models are used for REST in simulation, the volume conduction model matching problem still cannot be neglected (Hu et al., 2018). Within this framework and using extensive simulations, the performances of AR, REST, rAR, and rREST are compared in terms of the relative error of estimation of φ_∞. Additionally, in these simulations, we explore the performances of the model selection criteria for selecting the regularization parameter.

Finally, we assess their performances of all the estimators using real EEG data from 89 subjects with both regularization and volume conduction matching problem tested exhaustively.

Materials and Methods

Unified Reference Estimator

In what follows, we denote scalars with lowercase italic symbols (e.g., x), vectors with lowercase bold (e.g., x), matrices with uppercase bold (e.g., X); unknown parameters will be denoted by Greek letters (e.g., ξ). Furthermore, 1 is the vector of ones; I_Ne is a N_e by N_e identity matrix; N(μ, Σ) is the multivariable Gaussian distribution with mean vector μ and covariance matrix Σ; (·)^T is the transpose of (·); X⁺ is the pseudo-inverse of X; tr(·) is the trace of (·); $\widehat{\text{X}}$ is the estimation of X; ||·||₂, ||·||_F, and ||·||_M are the Euclidean norm, the Frobenius norm, and the Mahalanobis norm, respectively.

General Reference Model

The EEG is always recorded with respect to a time-varying reference. This is usually modeled as a constant subtracted from all electrodes at each instant. In the general case, we can consider that there are two separate reference constants, one for the scalp EEG signal, and another for the sensor noise (if they come from distinctly different source). In this case, the online recorded EEG signal at a given instant is modeled as

$\begin{array}{ll}v=\phi -1\times \rho +\epsilon -1\times \zeta & (1)\end{array}$

where φ is the pure EEG signal with the neutral reference over N_e electrodes, i.e., abovementioned φ_∞, and its distribution is N(0, Σ_φφ); ε is the sensor noise with $N(0,{\sigma }^2{\text{I}}_{N_e})$; ρ and ζ are two reference constants of EEG signal φ and sensor noise ε, respectively.ρ is assumed from a cephalic source, but ζ may come from either cephalic, non-cephalic or the coupled sources. Due to the uncertainty of these constants, the reference of v is an unknown variable. Note that other distributions for φ and ε may be used with our same general framework.

Applying a reference process is just a linear transformation of EEG data. Formally, it is the pre-multiplication of the reference transformation matrix with the EEG data. Thus, supposing a reference transformation matrix H = I − 1f^T (Hu et al., 2018), a referential recording is

$v_r=Hv=H(\phi +\epsilon )-(I-1f^T)\times 1\times (\rho +\zeta )$

Notably, the equation f^T 1 = 1 is satisfied for all the unipolar references, such as monopolar recording references (e.g., Cz, Fz, Oz, etc.), linked mastoids and average reference.

Thus, the general EEG reference model becomes

$\begin{array}{ll}{v}_r=H\phi +e,e=H\epsilon & (2)\end{array}$

where r denotes a specific reference. Note that H is a matrix of the rank as N_e − 1. Thus, the estimate of φ is transformed into an undetermined generalized linear inverse problem.

By means of maximum a posterior estimation (Murphy, 2012), or maximum penalized likelihood estimation (LaRiccia and Eggermont, 2009), we have the objective function.

$\begin{array}{ll}l={(v_r-H\phi )}^T{\sum }_{ee}^{+}(v_r-H\phi )+{\phi }^T{\sum }_{\phi \phi }^{+}\phi & (3)\end{array}$

After finding the partial derivative of (3) with respect to φ, it follows that

$\widehat{\phi }={(H^T{\sum }_{ee}^{+}H+{\sum }_{\phi \phi }^{+})}^{+}{H}^T{\sum }_{ee}^{+}{v}_r$

Referring to the matrix inversion lemma (Hager, 1989; Tarantola, 2005), $\widehat{\phi }$ is re-expressed as

$\begin{array}{ll}\widehat{\phi }={\sum }_{\phi \phi }{H}^T{(H{\sum }_{\phi \phi }{H}^T+{\sum }_{ee})}^{+}{v}_r& (4)\end{array}$

which is taken as the unified Bayesian estimator in reconstructing EEG potentials at infinity.

To derive the explicit expression of (4), in addition to assuming ${\Sigma }_{\text{e}\text{e}}={\sigma }^2\text{H}{\text{H}}^{\text{T}}$, Σ_φφ is assumed to have one of the following two different forms.

Uncorrelated Prior

$\begin{array}{ll}{\sum }_{\phi \phi }={\alpha }^2{I}_{N_e}& (5)\end{array}$

which means that the EEG potentials φ have independent priors across all the channels; α² is the mean of variances of the potentials over each electrode.

Substituting (5), v_r = Hv and ${\Sigma }_{\text{e}\text{e}}={\sigma }^2\text{H}{\text{H}}^{\text{T}}$ into (4), it becomes

$\begin{array}{ll}\widehat{\phi }=H^{+}Hv/(1+{\sigma }^2/{\alpha }^2)& (6)\end{array}$

We show that ${\text{H}}^{+}\text{H}={\text{I}}_{N_e}-11^{\text{T}}/N_e$ which is the average reference transforming matrix in the Appendix. Defining the sensor noise to the scalp EEG signal ratio as $ns{r}_1={\sigma }^2/{\alpha }^2$ and ${\text{H}}_{ar}={\text{I}}_{N_e}-11^{\text{T}}/N_e$, (6) is rewritten as

$\begin{array}{ll}\widehat{\phi }=H_{ar}v/(1+ns{r}_1)& (7)\end{array}$

which we shall call the regularized average reference (rAR). It is obvious that the usual AR is the special case of rAR when nsr₁ = 0.

Correlated Prior

$\begin{array}{ll}{\sum }_{\phi \phi }=K_{\infty }{\sum }_{jj}{K}_{\infty }^T& (8)\end{array}$

which models the EEG potentials across all the channels as correlated due to the effect of volume conduction on neural current sources, that is, we assume that φ = K_∞j; K_∞ is the lead field matrix with infinity reference; j is the primal current density of the neural current sources with $\text{j}~N(\text{0},{\beta }^2{\text{I}}_{N_s})$; N_s is the number of neural current sources; β² is the variance of the multivariate Gaussian signal j.

(4) is transformed by substituting (8) and defining K_r = HK_∞ as

$\begin{array}{ll}\widehat{\phi }=K_{\infty }·{\sum }_{jj}{K}_r^T{(K_r{\sum }_{jj}{K}_r^T+{\sum }_{ee})}^{+}{v}_r& (9)\end{array}$

which is the estimator for reconstructing the EEG potentials at infinity named as the regularized reference electrode standardization technique (rREST). This process can be interpreted as processing in two stages,

$\begin{array}{l}\text{Stage}1:\widehat{j}={\sum }_{jj}{K}_r^T{(K_r{\sum }_{jj}{K}_r^T+{\sum }_{ee})}^{+}{v}_r\\ \text{Stage}2:\widehat{\phi }=K_{\infty }\widehat{j}\end{array}$

the first one of which is solving the inverse problem with lead field K_r that has the same reference as the EEG potentials v_r and the second one of which is taking the forward calculation to reconstruct the EEG potentials with the theoretical neutral infinity reference. In stage 1, $\widehat{\text{j}}$ is the standard form of solving linear inverse problems and the reference problem, simultaneously.

Defining the sensor noise to the brain source signal ratio as $ns{r}_2={\sigma }^2/{\beta }^2$ and plugging ${\Sigma }_{\text{j}\text{j}}={\beta }^2{\text{I}}_{N_s}$, ${\Sigma }_{\text{e}\text{e}}={\sigma }^2\text{H}{\text{H}}^{\text{T}}$ into (9), it becomes

$\begin{array}{ll}\widehat{\phi }=K_{\infty }·K_r^T{(K_r{K}_r^T+ns{r}_2·H{H}^T)}^{+}{v}_r& (10)\end{array}$

which is the solution to reconstruct the EEG potential at infinity through solving the inverse solution by incorporating the identity diagonal structure of Σ_jj. Apparently, REST (Yao, 2001) $\widehat{\phi }={\text{K}}_{\infty }·{\text{K}}_r^{+}{\text{v}}_r$ is the special case of rREST when nsr₂ = 0 in (10).

For clarity, we summarize the general reference model and unified reference estimator in Table 1.

TABLE 1

Table 1. EEG reference model, unified estimator and schemes.

Reference Evaluation

Table 1 shows that both AR and REST are special cases of rAR and rREST if either the sensor noise is supposed to be zero or no regularization is applied. In this section, after transforming the general reference model into the standard ridge regression form, we evaluate the references via statistical model selection criteria.

Standard Regression Form

The objective function (3) of reference estimation is equivalent to the general ridge regression form (Chung et al., 2014)

$\begin{array}{ll}\widehat{\phi }(\lambda )=\underset{\phi }{\arg \min }\{{\Vert v_r-H\phi \Vert }_M^2+\lambda {\Vert L\phi \Vert }_2^2\}& (11)\end{array}$

where λ ≥ 0 is the regularization parameter; L is the regularization matrix. For convenience, we call the regularization of λ and L as “parameter regularization” and “structure regularization,” respectively.

Ridge regression is the name in statistics for Tikhonov regularization (Hoerl and Kennard, 1970). The difference between the general and the standard form of ridge regression is whether the regularization matrix L is identity and the misfit term is the Euclidean norm (Chung et al., 2014). Thus, we redefine φ′ = Lφ to make the regularization matrix being identity, and e′ = D^T U^T e (decompose HH^T = USU^T and S⁺ = DD^T) to transform the Mahalanobis norm of the misfit term as the Euclidean norm. To the end, the standard ridge regression form is

$\begin{array}{ll}{\widehat{\phi }}^{\prime }(\lambda )=\underset{{\phi }^{\prime }}{\arg \min }\{{\Vert v_r^{\prime }-H^{\prime }{\phi }^{\prime }\Vert }_2^2+\lambda {\Vert {\phi }^{\prime }\Vert }_2^2\}& (12)\end{array}$

with v_r′ = D^TU^T_vr and H′ = D^T U^T HL⁺. Then, the posterior mean of φ′ given v_r′ is

$\begin{array}{ll}{\widehat{\phi }}^{\prime }={({H^{\prime }}^T{H}^{\prime }+\lambda I_{N_e})}^{+}{H^{\prime }}^T{v}_r{}^{\prime }& (13)\end{array}$

then, the estimate of φ is $\widehat{\phi }={\text{L}}^{+}{({{\text{H}}^{\prime }}^{\text{T}}{\text{H}}^{\prime }+\lambda {\text{I}}_{N_e})}^{+}{{\text{H}}^{\prime }}^{\text{T}}{\text{v}}_r^{\prime }$ which is equivalent to the formula (10).

Model Selection Criteria

Since ridge regression is a linear estimator ($\widehat{\text{v}}$_r′ = Pv_r′) with $\text{P}={\text{H}}^{\prime }{({{\text{H}}^{\prime }}^{\text{T}}{\text{H}}^{\prime }+\lambda {\text{I}}_{N_e})}^{+}{{\text{H}}^{\prime }}^{\text{T}}$ where P is the projection (“hat”) matrix. The residual sum square error (RSS) is defined as

$\text{RSS}={\sum }_{t=1}^{N_t}\Vert v_{rt}{}^{\prime }-H^{\prime }{\widehat{\phi }}_t{}^{\prime }{\Vert }_2^2$

where v_rt′ and $\widehat{\phi }$_t′ with subscript t denote v_r′ and ${\widehat{\phi }}^{\prime }$ at the $t^{th}(t=1,\cdots ,N_t)$ time sample, respectively; N_t is the number of time samples in the whole EEG recording.

Under the standard ridge regression form (12), we explore three information criteria for the model selection: generalized cross-validation (GCV) (Chung et al., 2014), Akaike information criteria (AIC), and Bayesian information criteria (BIC) (Konishi and Kitagawa, 2008) to compare the reference schemes in Table 1. To apply these, we define the degree of freedom (DF) as

$\text{DF}(\lambda )=tr(P)={\displaystyle {\sum }_{i=1}^{N_e}\frac{s_i}{s_i+\lambda }}$

where {s_i} are the eigenvalues of H′^TH′. Since EEG reference acts as adding or subtracting a time-varying constant over all sensors at each instant, this instantaneous effect results in the dynamical alteration in the temporal domain. To investigate the difference of references, we extend the model selection criteria from single instant to the whole recording, approximately. Predefining N_et = N_e·N_t, GCV, AIC and BIC are expressed as

$\begin{array}{ll}\text{GCV}(\lambda )=\text{RSS}/{(N_{et}-\text{DF})}^2& (14)\end{array}$

$\begin{array}{ll}\text{AIC}(\lambda )=N_e{}_t\log (\text{RSS}/N_{et})+N_t·2·\text{DF}& (15)\end{array}$

$\begin{array}{ll}\text{BIC}(\lambda )=N_e{}_t\log (\text{RSS}/N_{et})+N_t·\text{DF}·\log (N_{et})& (16)\end{array}$

Note that GCV, AIC, and BIC at a single instant are the special cases of (14–16) with N_t = 1, respectively.

Regularization Parameter

The regularization parameter λ balances the goodness of fitting (i.e., likelihood) and the prior constraint on the EEG potentials at infinity. One may try to interactively estimate the hierarchal Bayesian hyperparameter via iteration (MacKay, 1992; Trujillo-Barreto et al., 2004). However, this may work for rREST but poorly for rAR because the noise term will be assimilated into the pure EEG signal in (6) due to the uncorrelated covariance prior. Namely, the objective function of AR is non-convex and it cannot converge to a global or local optimal point. Thus, we adopt a search strategy to explore how DF, GCV, AIC and BIC vary with the values of λ (Phillips et al., 2005). The idea is to plot the DF against λ, as well as GCV, AIC and BIC against DF. The theoretical solutions (7) and (10), indicate that the optimal λ is around nsr₁ for rAR, and approximates to nsr₂ for rREST, respectively. Since volume conduction acts as a lowpass spatiotemporal filter, it results in nsr₂ ≪ nsr₁ (Srinivasan et al., 1998; Stinstra and Peters, 1998; Srinivasan, 1999; Nunez and Srinivasan, 2006). Supposing the intervals of SNR are [35, 10] dB for rREST, and [30, −10] dB for rAR, we generate 1,000 values of λ from 1e-3.5 to 1e-1 for rREST, and from 1e-3 to 10 for rAR, by using sampled logarithm, respectively.

In the simulation, we can evaluate the reference estimators with an “oracle” regularization parameter, namely, one for which the smallest relative error regarding the ground truth. Additionally, the efficacy of the model selection criteria (GCV, AIC, and BIC) for selecting the regularization parameter is evaluated. It will be trickier to find a proper λ with actual EEG data where the ground truth is unknown. The value with which one model selection criteria reaches to a global or local minimum is regarded as the optimal λ chosen by the model selection criteria for actual EEG data.

It has been suggested to avoid regularization when applying REST so as not to lose high-frequency information (Yao, 2001). Instead, a truncation of singular value decomposition (TSVD) was proposed to suppress the effect of sensor noise for REST (Zhai and Yao, 2004). Therefore, we empirically adopt the recommended truncation parameter 0.05 for REST but use the model selection criteria for rREST.

Regularization Matrix

The choice of regularization matrix L depends on the prior covariance structure of the potentials at infinity. Table 1 shows that the prior covariance structure of φ as Σ_φφ = α² I_Ne for AR and rAR, and ${\Sigma }_{\phi \phi }=K_{\infty }{\Sigma }_{jj}{K}_{\infty }^T$, for REST and rREST, respectively. Therefore, the choices of L are:

$\begin{array}{ll}\begin{array}{l}\text{for}\text{AR}\text{and}\text{rAR},L_{ar}=I_{N_e}\\ \text{for}\text{REST}\text{and}\text{rREST},L_{rt}={\left[{(K_{\infty }{K}_{\infty }^{\text{T}})}^{+}\right]}^{1/2}\end{array}& (17)\end{array}$

Several cases of K_∞ are detailed in the next section. The degree of faithful biophysical regularization by L_rt increases from the less realistic approximation of volume conductor to the more realistic one.

Volume Conduction Model

For rREST, we study the volume conduction model matching problem, that is, to what extent, the lead field for rREST may be different from the actual one that generated the simulated EEG data or the real EEG recordings. Here, we evaluate several types of lead fields. The well-known spherical lead field (SLF) is a frequently adopted standard lead field. The most precise volume conduction model is the individual lead field (ILF) matched to the structural Magnetic Resonance Image (sMRI) of each subject. We also evaluate the average lead field (ALF) as a substitute for the individual lead field. Finally, we evaluate the “sparse individual lead field” (sILF) for which we switch off the voxels not used in the simulations. We will use suffixes to distinguish between types of lead fields.

Spherical lead field (SLF)

${\text{K}}_{\infty }^{\text{s}}$, is estimated based on the standard 3-layers concentric spherical head model comprising of brain, skull, and scalp with the conductivities being 1, 0.0125, and 1, respectively. For the spherical head shape, the radii are 1.0, 0.92, and 0.87 for the scalp surface, outer and inner skull surface, separately. The source space consists of 2600 discrete dipole sources evenly and radially distributed on the cortical surface with radius r = 0.86 and 400 discrete dipole sources uniformly located perpendicularly to the transverse plane with Z = −0.076. Here, the values of conductivities, radii, and coordinates are not the actual measurements but the relative ratios of conductivities and radii between the head layers, and the relative coordinates in the unit sphere space (Yao, 2001; Hu et al., 2018).

Individual lead field (ILF)

${\text{K}}_{\infty }^{\text{i}}$, is defined by normalization as

$K_{\infty }^{\text{i}}=K_{\infty }^{\text{iraw}}/{[tr(K_{\infty }^{\text{iraw}}{K}_{\infty }^{\text{iraw}}{}^T)]}^{1/2}$

where ${\text{K}}_{\infty }^{\text{iraw}}$ is the raw individual lead field matched to the $i^{th}(i=1,\cdots ,N_b)$ subject who underwent the EEG recording in Cuban Human Brain Mapping Project (CHBMP) (Uludag et al., 2009; Valdés-Hernández et al., 2010; Hernandez-Gonzalez et al., 2011; Bosch-Bayard et al., 2012). It is estimated by the finite element method based on the segmented cortical surface through CIVET pipeline (Yasser Ad-Dab'bagh, 2006) with sMRI. The cortical surface is formed by 6003 vertices and 11998*3 faces. In total, 6003 dipole sources are located at the vertices and activated perpendicularly to the cortical surface, individually. The normalization allows for comparison across subjects.

Average lead field (ALF)

${\text{K}}_{\infty }^a$, is the average of all the normalized ILFs of N_b subjects as

$K_{\infty }^a=\frac{1}{N_b}{\displaystyle {\sum }_{i=1}^{N_b}{K}_{\infty }^{\text{i}}}$

Sparse individual lead field (sILF)

${\text{K}}_{\infty }^{\text{si}}$, (for use in simulation) is obtained after transforming ILF as follows

$K_{\infty }^{\text{si}}=[K_{\infty }^{\text{iraw}}\circ W_i]/{[tr(K_{\infty }^{\text{iraw}}\circ W_i{W}_i^T\circ K_{\infty }^{\text{iraw}}{}^T)]}^{1/2}$

where ◦ means the matrix elementwise multiplication operation (i.e., Hadamard product); W_i is a matrix that consists of binary weights and has the same size with ${\text{K}}_{\infty }^{\text{iraw}}$; the entries at the columns of un-activated brain sources are zeros and the other columns are full of ones. In the simulation, the position of two patches of sources is incorporated into the covariance of the EEG potentials at infinity for rREST. In place of adopting l₀ norm or l₁ norm to sparse the brain electrical source signal j, we set the entries corresponded to non-activated sources of ILF being zeros to constrain the brain source signal indirectly.

Results

Simulation

EEG Generation

The simulation scheme is based on the forward equation below,

$\begin{array}{ll}\{\begin{array}{l}{v}_r=H\phi +H\epsilon ,\phi =K_{\infty }^{\text{iraw}}j\\ \text{SNR}=10{\log }_{10}({\alpha }^2/{\sigma }^2)\end{array}& (18)\end{array}$

where v_r is the simulated EEG potentials with unipolar reference; without loss of generality, the linear combination vector f = [0, ⋯, 0, 1]^T with the last entry being one and the others being zeros; two patches consisting of 150 dipole sources in j are activated, meeting 4-order bivariate autoregressive model (Figure 3); SNR is the scalp EEG signal to the sensor noise variance ratio in dB unit.

FIGURE 3

Figure 3. Source activities for EEG simulation. (A) Location of two source patches; (B) activities of two source patches.

With the ${\text{K}}_{\infty }^{\text{iraw}}$ of 89 subjects from the CHBMP database, the simulated EEG data of one group is 89 samples *58 channels *5120 instants. Totally, we generated the dataset A: 4 groups where the SNR values are different between groups but the same for all the samples in each group, and the dataset B: one group where the SNR values are different for all the samples. The simulation provides the ground truth of EEG potentials with the neutral reference, thus making it possible to intuitively compare the performances of references in terms of the relative error of potentials. For each data sample, the relative error (RE) of potentials is defined as

$\begin{array}{ll}\text{RE}=\Vert \widehat{\phi }-\phi {\Vert }_F^2/\Vert \phi {\Vert }_F^2& (19)\end{array}$

where φ denotes the ground truth; $\widehat{\phi }$ is the EEG potentials estimated by the references in Table 1.

Relative Error of Reference Estimators

The relative error (RE) is calculated using the simulated dataset A of 4 groups where the SNR values are 20dB, 8dB, 4dB, and 2dB for each group, respectively. Figures 4A–D show the REs of the reference estimators, including the lead fields variants (SLF, ILF, ALF, and sILF) for REST and rREST. Boxplots in black, green, red, and blue, show the REs of AR, rAR, REST and rREST, separately. It is evident from the boxplots (Figures 4A–D), that the REs of regularized references (rAR, rREST) are always less than that of unregularized references (AR, REST). Unpaired t-tests were applied to check the differences between unregularized references (AR, REST) and regularized references (rAR, rREST). Figure 4E lists the statistical significance levels (p-values) between AR and rAR, as well as between REST and rREST with various lead fields tested, separately. Except for the case between AR and rAR with SNR = 20dB, the p-values all reach very small values (<1e-7).

FIGURE 4

Figure 4. Relative error (RE) of reference estimators. (A–D) The boxplots of REs with SNR = 20, 8, 4, and 2 dB, respectively. Volume conduction model tested for REST and rREST are spherical lead field (SLF), individual lead field (ILF), average lead field (ALF), and sparse individual lead field (sILF). (E) The p-values of REs between ordinary references (AR, REST) and regularized references (rAR, rREST) under different SNR and various lead fields, separately.

With regularization, the decreases of REs from REST to rREST are more obvious than the decreases of REs from AR to rAR. Especially, regularization with sILF is much more effective than SLF, ILF, and ALF. This is not surprising since the sparse prior information was incorporated into the covariance structure. By contrast, by the simplest volume conduction model, i.e., SLF, the REs of rREST seems to be even larger than that of AR, and REST performs worst among all the references, when SNR = 20 and 8 dB. Comparing the REs by sILF and that of rREST by SLF with the REs of REST by SLF, we found that structure regularization by precise covariance seems to be more effective than the parameter regularization by selecting the optimal λ which led the least RE among all the tested values of λ. And the REs of rREST with sILF are the least among all the REs of other references. This means that structure regularization in combination with parameter regularization will have the best effect. In addition, injecting higher sensor noise with SNR being from 20 to 2 dB, the REs of rAR increase from less than 15% to higher than 60% accordingly, while the REs of rREST with SLF excluded rise from 4.1 to 40%.

These results indicate that: (1) except for the case of AR and rAR with SNR = 20 dB, AR, rAR, REST, and rREST by using SLF that roughly approximated the actual volume conduction model may be not able to reconstruct the EEG signal at infinity due to the quite large REs; (2) the effects of REST and rREST are volume conduction model dependent; (3) stronger regularization applied, better effect of rREST obtained; (4) for REST and rREST, the REs by using ALF seems to be almost same with the REs by ILF; (5) rAR may not have the effect of denoising. Over all, AR and rAR may be the alternative option when SNR is very large (≥20 dB), while rREST with precise volume conduction model should be the first option to estimate the EEG signal at infinity.

Model Selection for Estimators With Simulated Data

The model selection is analyzed using the simulated dataset B where SNR values uniformly distributed in the interval of [5 20] dB are set for the 89 samples to mimic the different SNRs of subjects in the real EEG recordings. The results summarized in Figure 5 allow determining the optimal reference via the model selection criteria. The plotted DF (degree of freedom), RSS (residual sum square), and the model selection criteria (GCV, AIC, BIC) are the average of them explored individually over the 89 data samples with all the regularization parameters λ (i.e., LMD) tested. The curves in Figure 5A show how the DF and GCV vary with the LMDs as well as how the RSS changes with DF. It is easy to see that the DF of rREST are always smaller than the DF of rAR. This means that rREST adopts the simpler model to reconstruct the EEG signal at infinity but employ the more realistic prior information for regularization than rAR. The lower RSS of rREST than rAR indicate that the EEG signal reconstructed by rREST is closer to the truth compared with the EEG signal restored by rAR. The curves in Figure 5B display how the model selection criteria (GCV, AIC, BIC) vary with the DF. Apparently, the model selection criteria values of rREST are always smaller than them of rAR. The prevalent lower values of model selection criteria provide the evidence to prefer rREST over rAR.

FIGURE 5

Figure 5. Model selection with simulated data. DF, RSS, GCV, AIC, and BIC are the average over the 89 EEG data samples individually simulated with a different SNR among [5 20] dB. (A), DF and GCV against LMD, and RSS varying with DF; (B), model selection criteria (GCV, AIC, and BIC) against DF. DF, degree of freedom; RSS, residual sum square; GCV, generalized cross validation; AIC, Akaike information criteria; BIC, Bayesian information criteria; SLF, spherical lead field; ILF, individual realistic lead field; ALF, the average of realistic lead fields; sILF, sparse individual lead field.

Regularization Parameter

For rREST, it is crucial to pick the best regularization parameter, i.e., the value of λ. Figure 6 displays that, to what extent, the values of λ selected by the model selection criteria (GCV, AIC, and BIC) are close to the truth, that is, the oracle picked by the least RE based on the simulated dataset B. Note that the best λ identified by ground truth and the model selection criteria is in the individual data sample level rather than in the group level due to the averaged model selection criteria curves. Comparing the mean relative error (mRE) and the standard deviations in Figures 6B–D with those in Figure 6A, GCV is easily found as the best one to select the proper regularization parameter due to the almost same mRE and standard deviations to the truth; AIC is worse than GCV since except for the rREST by using sILF, the regularized reference (rAR, rREST) show the same or larger mRE and standard deviations than the ordinary reference (AR, REST); BIC is the worst one to select the proper regularization parameters because all the regularized references (rAR, rREST) present the larger mREs and standard deviations than the ordinary references (AR, REST).

FIGURE 6

Figure 6. Regularization parameter selection. Each square and the error bar are the mean relative error (mRE) and standard deviation over the 89 data samples individually simulated with a different SNR among [5 20] dB. (A), The truth where the best λ is picked by the least RE; (B–D), the results where the best λ is selected by the least GCV, AIC, and BIC values, respectively.

Model Selection for Estimators With Real Data

We take the EEG of 89 subjects from the CHBMP database to evaluate the reference estimators. The EEG recordings were carried out in accordance with the recommendations of Ethics committees of Ministry of Public health and Cuban Neuroscience Center with written informed consent from all subjects. The EEG was acquired with 58 channels, 10–10 electrode placement system, sampling rate 200 Hz, recording period 2.5–5 minutes, and with the resting-state of “eyes-closed-open” intermittently conditioned. Figure 7 displays a real EEG data sample. To compare the performance of references over all subjects, we normalize the EEG data by

$v_r=v_r^{raw}/\Vert v_r^{raw}{\Vert }_F$

The model selection criteria GCV, AIC, and BIC are calculated for each subject with the matched ILF, sILF, the identical ALF and SLF. The mean model selection criteria are averaged over 89 subjects.

FIGURE 7

Figure 7. The illustration of a real EEG data sample. (A) The electrodes position, (B) the waveform of one epoch.

To validate the preference of rREST over rAR, a similar analysis as descried before for simulations regarding performance is now applied to the realistic resting state EEG dataset of 89 subjects from the CHBMP database. The results in Figure 8 show how the DF (degree of freedom) and GCV (generalized cross validation) change with the regularization parameter λ (i.e., LMD) and how the RSS (residual sum square) and model selection criteria (GCV, AIC, and BIC) vary with the DF. The plotted DF, RSS, and model selection criteria are the average of them individually explored over the 89 subjects. Since the reference transformation matrix H and the limits of regularization parameter λ (i.e., LMD) are the same as that in simulation, the curves of DF against LMDs in Figure 8A are identical to them in Figure 5A. The lower RSS and model selection criteria curves in Figure 8B confirm, for real data, our previous findings in simulations, that (1) the EEG signal reconstructed by rREST has lower RSS than that restored by rAR; (2) rREST has the smaller values of GCV, AIC, and BIC than rAR, except for the almost same BIC between rREST and rAR around DF = 28. Since GCV was found to be the best criteria to choose the proper regularization parameter λ in the simulation shown in Figure 6, we suggest adopting GCV to select the value of λ for each EEG recording in practice when the ground truth is unknown. The curves in the middle of Figure 8A shows how GCV varies with the values of λ (i.e., LMDs). For rREST, the global minimum of GCV occurs around λ = 1e-2 or DF = 10. We therefore conjecture that for the group analysis, it is possible to find the best regularization parameter by the lowest GCV which may offer an empirical example for the possible optimal LMD for rREST. However, it will be the safest to choose the regularizer parameter by using GCV for each individual recording as what we have shown in Figure 6. By contrast, GCV of rAR seems to be one constant which is caused by the nonconvex solution of rAR where it is hard to find the proper λ. These results indicate that the preference of rREST over rAR is validated for real EEG. Moreover, for rREST, the best regularization parameter can be picked at the global minimum of GCV curves.

FIGURE 8

Figure 8. Model selection with actual EEG. DF, RSS, GCV, AIC, and BIC are the average of them individually explored over the 89 real EEG recordings. (A) DF and GCV against LMD, and RSS varying with DF; (B), model selection criteria (GCV, AIC, and BIC) against DF. DF, degree of freedom; RSS, residual sum square; GCV, generalized cross validation; AIC, Akaike information criteria; BIC, Bayesian information criteria; SLF, spherical lead field; ILF, individual realistic lead field; ALF, the average of realistic lead fields; sILF, sparse individual lead field.

Discussion

Although the reference electrode standardization technique (REST) was put forward some time ago, its theoretical underpinnings have not been deeply studied, particularly from a mathematical statistics perspective. Prior to REST, the average reference (AR) had been broadly adopted, e.g., in the microstates analysis (Khanna et al., 2015), and as well offered as the final solution to “reference electrode problem” in the inverse solutions (Pascual-Marqui, 1999, 2002, 2007; Pascual-Marqui et al., 2011). Currently, both AR and REST are the main competing estimators (Nunez, 2010). Many comparative studies have been carried out, however, without providing the definitive evidence to prefer one over another (Qin et al., 2010; Chella et al., 2016, 2017; Lei and Liao, 2017). The need to settle this issue has been reinforced recently by the theoretical results of Yao, who demonstrated that the main assumption of AR—the cancelation of brain potentials averaged over the scalp is, in general, false (Yao, 2017). However, it is difficult to decide the reference issue solely by the biophysical interpretations. Empirical verification of the best reference using a full statistical model is also essential.

In this study, we propose to view the estimation of the potentials at infinity and the determination of reference as a linear inverse problem that can be attacked using well known Bayesian techniques. To our surprise, both AR and REST are two special cases with different prior distributions for the covariance of the EEG potentials referenced to infinity. By explicitly introduced the sensor noise term into the reference model, we combined the estimation of the potentials at infinity with denoising. The formulation, based on maximum a posterior estimation, leads to the regularized estimators, rAR and rREST. Finally, recognizing that the reference is a linear inverse problem leads to the use of model selection criteria to examine several issues. It is found, for both simulated and actual data, that (1) regularization is critical to solving the reference problem and denoising simultaneously; (2) the regularized reference (rAR/ rREST) has a better performance than the ordinary reference (AR/REST), respectively; (3) rREST outperforms rAR; (4) to apply rREST to real EEG data, generalized cross-validation is recommended as an effective measure to select the optimal regularization parameter. In our opinion, the definitive argument in favor of rREST is that for 89 resting state EEGs it provides a smaller Generalized Cross Validation. This is the first empirical comparison of references using an information criterion that approximates the Bayesian model evidence.

We have demonstrated that AR is not “the final solution to the reference electrode problem” (Pascual-Marqui et al., 2011), but rather a special case of uncorrelated prior and noise-free of the unified reference estimator. Pascual-Marqui et al derived AR with the assumption of exact noise-free, or say, under the condition that the covariance matrix for the measurement noise is identity. Before the inverse solution, the reference problem is solved by the derived AR which is supposed to achieve the best fit for the reference process. However, in this study, if the reference process is also involved for the measurement noise but not only for the EEG signal, AR cannot be derived. Since the inverse solution does not depend on the reference electrode, both AR and the monopolar reference can be used to transform the lead fields and the EEG recordings with the same reference before the inverse solution. This is described as the stage 1 of the implementation of rREST following the Equation (9).

REST has attracted attention due to its theoretically sound basis (Yao, 2001; Kayser and Tenke, 2010; Nunez, 2010). However, some studies with ordinary REST suggest that it does not uniformly dominate AR (Yao, 2001; Hu et al., 2018). Though ordinary REST was found to be more effective than AR for vertically oriented and shallow dipole sources, this was not so for transverse or deep dipoles. These findings, in our perspective, were due to that the covariance structure for ordinary REST is derived from the two factors, a spherical volume conductor, and limiting sources to the equivalent distributed dipoles layer (Yao and He, 2003), i.e., the sources over the 2D cortical sheet with radial orientation (Yao, 2001). By contrast, we tested here more realistic volume conduction models. Also, our simulations were based on multiple cortical patches. Our results indicate that more realistic volume conductor and source space do make the reference estimator better, and that in fact, REST does dominate AR in all cases.

We emphasize multi-possibilities of source modeling and restate the conception of the generalized inverse problem. As REST is a generalized inverse solution, the multipole sources (Daunizeau et al., 2006) and the general 3D distributed sources at each grid of brain volume (Michel et al., 2004) can be also be adopted for REST as well (Yao, 2000). We have shown that all generalized inverses are not equal and an interesting line of research will be to explore how different source model procedures can be translated into variants of REST.

The volume conduction model matching test showed that REST and rREST is volume conduction model dependent and the importance of the adequate volume conduction model. This is in agreement with the findings of Hu et al. (2018) and Liu et al. (2015). Liu et al reported that a realistic volume conduction model is critical to ordinary REST. Hu et al stated that ordinary REST is volume conduction model dependent but imprecise or slightly perturbated lead fields does not deteriorate ordinary REST much. This is in line with our simulation results that better estimates of both the volume conductor and source lead to better reference estimates. The result that sILF achieves the least relative errors among all the volume conduction models tested for rREST is obvious since the prior sILF matches the forward calculation. The only point of this simulation is to caution that the prior for rREST should be as close as possible to the actual. Of course, this can only be known to an actual approximation. In future, a promising way to account the uncertainty for a correct prior is to employ Bayesian Model Averaging (Trujillo-Barreto et al., 2004).

However, it is computational costly to estimate the individual lead field which requires the subject's sMRI—something not usual in many clinical settings. However, we found that the average lead field achieves the almost same performance as obtained with the exact individual lead field. This validates the proposal that approximate head models without individual MRI can be quite useful (Valdés-Hernández et al., 2009).

A critical point for rREST is to choose the optimal regularization parameter which has been the topic of intense research in statistics (Konishi and Kitagawa, 2008). Our simulations and validation on the real EEG recording suggest that the generalized cross validation criteria could be a simple and sensitive procedure to solve this issue. Once again, in practice, we recommend adopting the generalized cross validation criteria to select the regularization parameter for each individual EEG recording.

Main contributions of this paper are:

(1) We propose that reference estimation is a unified Bayesian linear inverse problem.

(2) This framework explicitly models sensor noise as a part of the EEG generative model.

(3) AR and REST are shown to be the special cases of the linear inverse problem, with a spatially independent prior for AR and a spatially correlated prior for REST due to volume conduction.

(4) Regularized estimators, rAR and rREST, are developed to implement reference estimation and denoising simultaneously.

(5) We adopted the model selection criteria (GCV, AIC, BIC) for not only to select the hyperparameter but also to compare model families. GCV was found to be the most useful indicator.

(6) We propose the average lead field as a practical substitute for the individual lead field to construct near optimal estimators.

Several extensions of this study are being explored:

• Artifact suppression may be incorporated together with reference estimation and denoising. For example, outliers can be eliminated by utilizing a likelihood function designed for robust statistics (Huber and Ronchetti, 2009).

• More biophysical information may be built into the prior distributions to more effectively differentiate the EEG signal from the sensor noise. Particularly, covariance matrices corresponding to different types of structured sparsity source models should be examined (Paz-Linares et al., 2017).

• We have dealt only with spatial priors for the covariance matrix of the EEG. Dynamical priors can easily be incorporated. Temporal autocorrelations may be modeled as state space models and estimated via the Kalman filter (Galka et al., 2004). Alternatively, formulations for the frequency, or time frequency domain are possible.

• The framework may be also extended to event related potentials incorporating prior work from our group in this direction (Carbonell et al., 2004).

Conclusion

We state the EEG reference problem as a unified inverse problem that can be solved via Bayesian techniques. To our best knowledge, this is a novel approach to the problem. This formulation allows us to:

• Adopt regularization methods to estimate the potential referenced to infinity.

• Demonstrate that REST and AR are the special cases of the unified estimator with different EEG spatial covariance priors.

• Simultaneously carry out denoising as part of the reference estimation procedure.

• Use model selection criteria to determine the optimal reference estimator. These results can be summarized as:

◦ Regularized references (rREST or rAR) are superior to the ordinary REST or AR, with rREST having the overall best performance for both simulations and real data.

◦ The optimal choice of volume conductor model is the individual or averaged lead field.

Regularized REST (rREST) may be used in clinical settings, as an improved estimator of EEG potentials referenced to infinity.

Author Contributions

SH developed the theories with PV-S, processed the data, and wrote the complete manuscript. DY made some valuable comments to this study. PV-S posed the scientific question to this study and polished the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors declare no conflict of interest. This research was co-funded by the National Natural Science Foundation of China projects (NSFC Grants No. 61673090, 81330032, and 81601585). We thank Esin Karahan and Pedro A. Valdes-Hernandez for their valuable technical assistance.

References

Baksalary, J. K., Maria Baksalary, O., and Trenkler, G. (2003). A revisitation of formulae for the Moore–Penrose inverse of modified matrices. Linear Algebra Appl. 372, 207–224. doi: 10.1016/S0024-3795(03)00508-1

CrossRef Full Text | Google Scholar

Bigdely-Shamlo, N., Mullen, T., Kothe, C., Su, K.-M., and Robbins, K. A. (2015). The PREP pipeline: standardized preprocessing for large-scale EEG analysis. Front. Neuroinform. 9:16. doi: 10.3389/fninf.2015.00016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosch-Bayard, J., Valdés-Sosa, P. A., Fernandez, T., Otero, G., Pliego Rivero, B., Ricardo-Garcell, J., et al. (2012). 3D Statistical Parametric Mapping of quiet sleep EEG in the first year of life. Neuroimage 59, 3297–3308. doi: 10.1016/j.neuroimage.2011.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbonell, F., Galán, L., Valdés, P., Worsley, K., Biscay, R. J., Diaz-Comas, L., et al. (2004). Random Field–Union Intersection tests for EEG/MEG imaging. Neuroimage 22, 268–276. doi: 10.1016/j.neuroimage.2004.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Chella, F., D'Andrea, A., Basti, A., Pizzella, V., and Marzetti, L. (2017). Non-linear analysis of scalp EEG by using bispectra: the effect of the reference choice. Front. Neurosci. 11:262. doi: 10.3389/fnins.2017.00262

PubMed Abstract | CrossRef Full Text | Google Scholar

Chella, F., Pizzella, V., Zappasodi, F., and Marzetti, L. (2016). Impact of the reference choice on scalp EEG connectivity estimation. J. Neural Eng. 13:36016. doi: 10.1088/1741-2560/13/3/036016

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, J., Español, M. I., and Nguyen, T. (2014). Optimal Regularization Parameters for General-Form Tikhonov Regularization. Available online at: http://arxiv.org/abs/1407.1911

Daunizeau, J., Mattout, J., Clonda, D., Goulard, B., Benali, H., and Lina, J. M. (2006). Bayesian spatio-temporal approach for EEG source reconstruction: conciliating ECD and distributed models. IEEE Trans. Biomed. Eng. 53, 503–516. doi: 10.1109/TBME.2005.869791

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferree, T. C., Luu, P., Russell, G. S., and Tucker, D. M. (2001). Scalp electrode impedance, infection risk, and EEG data quality. Clin. Neurophysiol. 112, 536–544. doi: 10.1016/S1388-2457(00)00533-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Galka, A., Yamashita, O., Ozaki, T., Biscay, R., and Valdés-Sosa, P. (2004). A solution to the dynamical inverse problem of EEG generation using spatiotemporal Kalman filtering. Neuroimage 23, 435–453. doi: 10.1016/j.neuroimage.2004.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, D. (1950). The clinical use of the “average” reference electrode in monopolar recording. Electroencephalogr. Clin. Neurophysiol. 2, 209–212.

PubMed Abstract | Google Scholar

Guruvareddy, A. (2013). Artifact Removal from EEG Signals. Int. J. Comput. Appl. 77, 975–8887. doi: 10.5120/13543-1175

CrossRef Full Text

Hager, W. W. (1989). Updating the inverse of a matrix. SIAM Rev. 31, 221–239. doi: 10.1137/1031049

CrossRef Full Text | Google Scholar

Hernandez-Gonzalez, G., Bringas-Vega, M. L., Galán-Garcia, L., Bosch-Bayard, J., Lorenzo-Ceballos, Y., Melie-Garcia, L., et al. (2011). Multimodal quantitative neuroimaging databases and methods: the cuban human brain mapping project. Clin. EEG Neurosci. 42, 149–159. doi: 10.1177/155005941104200303

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoerl, A. E., and Kennard, R. W. (1970). Ridge regression: biased estimation for nonorthogonal problems. Technometrics 12, 55–67. doi: 10.1080/00401706.1970.10488634

CrossRef Full Text | Google Scholar

Hu, S., Lai, Y., Valdes-Sosa, P. A., Bringas-Vega, M. L., and Yao, D. (2018). How do reference montage and electrodes setup affect the measured scalp EEG potentials? J. Neural Eng. 15:26013. doi: 10.1088/1741-2552/aaa13f

PubMed Abstract | CrossRef Full Text | Google Scholar

Huber, P. J., and Ronchetti, E. M. (2009). Robust Statistics. Hoboken, NJ: John Wiley and Sons, Inc.

Google Scholar

Kaipio, J., and Somersalo, E. (2007). Statistical inverse problems: discretization, model reduction and inverse crimes. J. Comput. Appl. Math. 198, 493–504. doi: 10.1016/j.cam.2005.09.027

CrossRef Full Text | Google Scholar

Kayser, J., and Tenke, C. E. (2010). In search of the Rosetta Stone for scalp EEG: converging on reference-free techniques. Clin. Neurophysiol. 121, 1973–1975. doi: 10.1016/j.clinph.2010.04.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Khanna, A., Pascual-Leone, A., Michel, C. M., and Farzan, F. (2015). Microstates in resting-state EEG: current status and future directions. Neurosci. Biobehav. Rev. 49, 105–113. doi: 10.1016/j.neubiorev.2014.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Konishi, S., and Kitagawa, G. (2008). Information Criteria and Statistical Modeling. New York, NY: Springer.

Google Scholar

Kulaichev, A. P. (2016). Optimal choice of a reference electrode for EEG recording. Moscow Univ. Biol. Sci. Bull. 71, 145–150. doi: 10.3103/S0096392516030068

CrossRef Full Text | Google Scholar

LaRiccia, V. N., and Eggermont, P. P. (2009). Maximum Penalized Likelihood Estimation. New York, NY: Springer.

Google Scholar

Lei, X., and Liao, K. (2017). Understanding the influences of EEG reference: a large-scale brain network perspective. Front. Neurosci. 11:205. doi: 10.3389/fnins.2017.00205

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemm, S., Curio, G., Hlushchuk, Y., and Müller, K.-R. (2006). Enhancing the signal-to-noise ratio of ICA-based extracted ERPs. IEEE Trans. Biomed. Eng. 53, 601–607. doi: 10.1109/TBME.2006.870258

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., and Yao, D. (2007). A new method of spatio-temporal topographic mapping by correlation coefficient of K-means cluster. Brain Topogr. 19, 161–176. doi: 10.1007/s10548-006-0017-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., Balsters, J. H., Baechinger, M., Van Der Groen, O., Wenderoth, N., and Mantini, D. (2015). Estimating a neutral reference for electroencephalographic recordings: the importance of using a high-density montage and a realistic head model. J. Neural Eng. 12:56012. doi: 10.1088/1741-2560/12/5/056012

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKay, D. J. C. (1992). Bayesian Interpolation. Neural Comput. 4, 415–447. doi: 10.1162/neco.1992.4.3.415

CrossRef Full Text | Google Scholar

Marzetti, L., Nolte, G., Perrucci, M. G., Romani, G. L., and Del Gratta, C. (2007). The use of standardized infinity reference in EEG coherency studies. Neuroimage 36, 48–63. doi: 10.1016/j.neuroimage.2007.02.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Michel, C. M., Murray, M. M., Lantz, G., Gonzalez, S., Spinelli, L., and Grave De Peralta, R. (2004). EEG source imaging. Clin. Neurophysiol. 115, 2195–2222. doi: 10.1016/j.clinph.2004.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, K. (2012). Machine Learning: A Probabilistic Perspective. Cambridge, MA: The MIT Press. Available online at: https://books.google.com/books?id=NZP6AQAAQBAJ&printsec=frontcover#v=onepage&q&f=false

Nunez, P. L. (2010). REST: a good idea but not the gold standard. Clin. Neurophysiol. 121, 2177–2180. doi: 10.1016/j.clinph.2010.04.029

CrossRef Full Text | Google Scholar

Nunez, P. L., and Srinivasan, R. (2006). Electric Fields of the Brain, 2nd Edn. New York, NY: Oxford University Press.

Offner, F. F. (1950). The EEG as potential mapping: the value of the average monopolar reference. Electroencephalogr. Clin. Neurophysiol. 2, 213–214. doi: 10.1016/0013-4694(50)90040-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Pascual-Marqui, R. D. (1999). Review of methods for solving the EEG inverse problem. Int. J. Bioelectromagn. 1, 75–86.

Google Scholar

Pascual-Marqui, R. D. (2002). Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find. Exp. Clin. Pharmacol. 24 (Suppl. D), 5–12.

PubMed Abstract | Google Scholar

Pascual-Marqui, R. D. (2007). Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization. Neurosci. Lett. 485, 198–203.

Google Scholar

Pascual-Marqui, R. D., Lehmann, D., Koukkou, M., Kochi, K., Anderer, P., Saletu, B., et al. (2011). Assessing interactions in the brain with exact low-resolution electromagnetic tomography. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 369, 3768–3784. doi: 10.1098/rsta.2011.0081

PubMed Abstract | CrossRef Full Text | Google Scholar

Paz-Linares, D., Vega-Hernández, M., Rojas-López, P. A., Valdés-Hernández, P. A., Martínez-Montes, E., and Valdés-Sosa, P. A. (2017). Spatio temporal EEG source imaging with the hierarchical bayesian elastic net and elitist lasso models. Front. Neurosci. 11:635. doi: 10.3389/fnins.2017.00635

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, C., Mattout, J., Rugg, M. D., Maquet, P., and Friston, K. J. (2005). An empirical Bayesian solution to the source reconstruction problem in EEG. Neuroimage 24, 997–1011. doi: 10.1016/j.neuroimage.2004.10.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, C., Rugg, M. D., and Friston, K. J. (2002). Systematic regularization of linear inverse solutions of the EEG source localization problem. Neuroimage 17, 287–301. doi: 10.1006/nimg.2002.1175

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, Y., Xu, P., and Yao, D. (2010). A comparative study of different references for EEG default mode network: the use of the infinity reference. Clin. Neurophysiol. 121, 1981–1991. doi: 10.1016/j.clinph.2010.03.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Robert, C. P. (2007). The Bayesian Choice. New York, NY: Springer New York.

PubMed Abstract | Google Scholar

Srinivasan, R. (1999). Methods to improve spatial resolution of EEG. Int. J. Bioelectromagn. 1, 102–111.

Google Scholar

Srinivasan, R., Nunez, P. L., and Silberstein, R. B. (1998). Spatial filtering and neocortical dynamics: estimates of EEG coherence. IEEE Trans. Biomed. Eng. 45, 814–826. doi: 10.1109/10.686789

PubMed Abstract | CrossRef Full Text | Google Scholar

Stinstra, J. G., and Peters, M. J. (1998). The volume conductor may act as a temporal filter on the ECG and EEG. Med. Biol. Eng. Comput. 36, 711–716. doi: 10.1007/BF02518873

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarantola, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia, PA: Society for Industrial and Applied Mathematics. doi: 10.1137/1.9780898717921

CrossRef Full Text | Google Scholar

Teplan, M. (2002). Fundamentals of EEG measurement. Meas. Sci. Rev. 2, 1–11.

Google Scholar

Tian, Y., and Yao, D. (2013). Why do we need to use a zero reference? Reference influences on the ERPs of audiovisual effects. Psychophysiology 50, 1282–1290. doi: 10.1111/psyp.12130

PubMed Abstract | CrossRef Full Text | Google Scholar

Trujillo-Barreto, N. J., Aubert-Vázquez, E., and Valdés-Sosa, P. A. (2004). Bayesian model averaging in EEG/MEG imaging. Neuroimage 21, 1300–1319. doi: 10.1016/j.neuroimage.2003.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Uludag, K., Evans, A. C., Della-Maggiore, V., Kochen, S., Amaro, E., Sierra, O., et al. (2009). Latin American brain mapping network (LABMAN). Neuroimage 47, 312–313. doi: 10.1016/j.neuroimage.2009.03.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdés-Hernández, P. A., Ojeda-González, A., Martínez-Montes, E., Lage-Castellanos, A., Virués-Alba, T., Valdés-Urrutia, L., et al. (2010). White matter architecture rather than cortical surface area correlates with the EEG alpha rhythm. Neuroimage 49, 2328–2339. doi: 10.1016/j.neuroimage.2009.10.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdés-Hernández, P. A., von Ellenrieder, N., Ojeda-Gonzalez, A., Kochen, S., Alemán-Gómez, Y., Muravchik, C., et al. (2009). Approximate average head models for EEG source imaging. J. Neurosci. Methods 185, 125–132. doi: 10.1016/j.jneumeth.2009.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yasser Ad-Dab'bagh (2006). “The CIVET image-processing environment: a fully automated comprehensive pipeline for anatomical neuroimaging research,” in Proceedings 12th Annual Meeting Organization for Human Brain Mapping, 2266.

Yao, D. (2000). High-resolution EEG mappings: a spherical harmonic spectra theory and simulation results. Clin. Neurophysiol. 111, 81–92. doi: 10.1016/S1388-2457(99)00205-9

CrossRef Full Text | Google Scholar

Yao, D. (2001). A method to standardize a reference of scalp EEG recordings to a point at infinity. Physiol. Meas. 22, 693–711. doi: 10.1088/0967-3334/22/4/305

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, D. (2017). Is the surface potential integral of a dipole in a volume conductor always zero? A cloud over the average reference of EEG and ERP. Brain Topogr. 30, 161–171. doi: 10.1007/s10548-016-0543-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, D., and He, B. (2003). Equivalent physical models and formulation of equivalent source layer in high-resolution EEG imaging. Phys. Med. Biol. 48, 3475–3483. doi: 10.1088/0031-9155/48/21/002

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, D., Wang, L., Oostenveld, R., Nielsen, K. D., Arendt-Nielsen, L., and Chen, A. C. N. (2005). A comparative study of different references for EEG spectral mapping: the issue of the neutral reference and the use of the infinity reference. Physiol. Meas. 26, 173–184. doi: 10.1088/0967-3334/26/3/003

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhai, Y., and Yao, D. (2004). A study on the reference electrode standardization technique for a realistic head model. Comput. Methods Progr. Biomed. 76, 229–238. doi: 10.1016/j.cmpb.2004.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix: Demonstration that H⁺H = I − 11^T/N_e

Recall f^T 1 = 1, H = I−1f^T (Hu et al., 2018) which reduces the matrix rank by one. Keeping same form with the formula (1.2) in Baksalary et al. (2003), we rewrite

$M=A+b{c}^T,\text{with}M=H,A=I,b=-1,c=f$

Referring to the Theorem 1.1 in Baksalary et al. (2003), it follows that rank(M) = rank(A)−1, since λ = 1+c^T A⁺ b = 1+f^T(−1) = 0 and both b and c belong to the column space of A = I. By utilizing the case (↓) of List 2 in Theorem 2.1 in Baksalary et al. (2003), we take

$M^{+}M=A^{+}-{\delta }^{-1}d{d}^T$

The formulas in (1.3) and (1.4) from Baksalary et al. (2003) are d = A^T b and δ = d^T d. Applying these relations to our problem, it turns to be

$H^{+}H=I-11^T/N_e$

which is the classical average reference.