Bridging M/EEG Source Imaging and Independent Component Analysis Frameworks Using Biologically Inspired Sparsity Priors

2.1 . The Measurement Model

In ESI, the electromagnetic activity of the brain sources is defined on a dense grid of known cortical locations (the source space). Typically, a vector of $N_y$ voltage (EEG) or magnetic field (MEG) measurements at sample $k$, ${\mathbf{y}}_k\in {\mathbb{R}}^{N_y}$, relates to the activity of $N_g$ sources, ${\mathbf{g}}_k\in {\mathbb{R}}^{3N_g}$, through the following instantaneous linear equation (Dale & Sereno, 1993),

where ${\mathbf{g}}_k$ is the vector of current source density along the three orthogonal directions and ${\mathbf{e}}_k\in {\mathbb{R}}^{N_y}$ represents the measurement noise vector. The source is projected onto the sensor space through the electric or magnetic lead field matrix $\mathbf{L}=[{\mathbf{l}}_{x_1},{\mathbf{l}}_{y_1},{\mathbf{l}}_{z_1}\dots ,{\mathbf{l}}_{x_{N_g}},{\mathbf{l}}_{y_{N_g}},{\mathbf{l}}_{z_{N_g}}]\in {\mathbb{R}}^{N_y\times 3N_g}$ ($N_y\ll N_g$), where each column ${\mathbf{l}}_i$ denotes the scalp projection of the ith unitary current dipole along a canonical axis. The lead field matrix is usually precomputed for a given model of the head derived from a subject-specific MRI (Hallez et al., 2007). Alternatively, if an individual MRI is not available, an approximated lead field matrix obtained from a high-resolution template can be used (Huang, Parra, & Haufe, 2016). The inverse problem of the M/EEG can be stated as the estimation of a source configuration ${\widehat{\mathbf{g}}}_k$ that is likely to produce the scalp topography ${\mathbf{y}}_k$. To simplify the exposition, from now on, we focus on EEG. However, our results hold for MEG as well, keeping in mind that in the formulas, we swap the respective voltage by magnetic field measurements and the electric by the magnetic lead field matrix.

As we mentioned in section 1, the noise term ${\mathbf{e}}_k$ is usually assumed to be white and gaussian. This simplification is acceptable when the measurements are not affected by nongaussian pseudo-random artifacts generated by eye blinks, lateral eye movements, facial and neck muscle activity, and body movement, among others. With this in mind, we start building our bridge between ESI and ICA by augmenting equation 2.1 with a term that accounts for the linear contribution of artifacts as follows,

where $\mathbf{\nu }\in {\mathbb{R}}^{N_{\nu }}$ is a vector of $N_{\nu }$ artifact sources and $\mathbf{A}=[{\mathbf{a}}_1,\dots ,{\mathbf{a}}_{N_{\nu }}]\in {\mathbb{R}}^{N_y\times N_{\nu }}$ is a dictionary of artifact scalp projections (see Figure 1).

The columns of $\mathbf{A}$ that correspond to muscle activity may be obtained based on a detailed electromechanical model of the body (Böl, Weikert, & Weichert, 2011). Another approach is to expand the lead field matrix to model putative EMG sources located all over the scalp (Janani et al., 2017). EOG sources can be modeled by including the projections of two current dipoles situated right behind the eyes (Fujiwara et al., 2009). Our approach to obtaining $\mathbf{A}$ is rather empirical, inspired by ICA's success for identifying artifact scalp projections.

In section 2.7, we explain how we built the artifact dictionary $\mathbf{A}$ using a set of stereotypical artifact scalp projections obtained by ICA, but for the sake of building up the theory, let's assume that $\mathbf{A}$ is known by the time the data are collected. With this in mind, we rewrite equation 2.2 in a compact manner as follows,

where $\mathbf{H}\triangleq [\mathbf{L},\mathbf{A}]$ is a known observation operator and ${\mathbf{x}}_k\triangleq {[{\mathbf{g}}_k^T,{\mathbf{\nu }}_k^T]}^T$ is the augmented vector of unknown (latent) brain and artifact sources (see Figure 1). Note that, structurally, the standard generative model in equation 2.1 and the augmented one in equation 2.3 are identical. However, they differ in that in equation 2.3, we are explicitly modeling the instantaneous spatial contribution of nonbrain sources to the scalp topography ${\mathbf{y}}_k$. The assumption of gaussian measurement noise yields the following likelihood function,

where $\mathbf{R}$ is the sensor noise covariance, which can be estimated from empty room measurements of MEG data, from prestimulus data (Engemann & Gramfort, 2015; Wipf, Owen, Attias, Sekihara, & Nagarajan, 2010), or learned adaptively (Ojeda et al., 2018). Here we assume that, after accounting for most artifactual sources, common-mode sensor noise is homoskedastic and spatially uncorrelated $\mathbf{R}=\lambda \mathbf{I}$.

2.2 . Block-Sparse Source Imaging

Equation 2.3 does not have a unique solution, so to obtain approximated source maps with biological interpretation, we need to introduce constraints. Since the neural generators of the EEG are the electrical currents produced by distributed neural masses that become locally synchronized in space and time (Nunez & Srinivasan, 2006), it makes sense to constrain source maps to be globally sparse (seeking to explain the observed scalp topographies by a few spots of cortical activity) and locally correlated (so that we obtain spatially smooth maps as opposed to maps formed by scattered, isolated sources). Artifactual sources, on the other hand, can be assumed to be spatially uncorrelated from one another and true brain sources.

We use the assumptions above to furnish the following source prior,

where $\mathbf{Q}$ is the source noise covariance, which we assume has the following block-diagonal structure,

where the upper block corresponds to the covariance of groups of brain sources and has also the following block-diagonal structure,

and ${\mathbf{\Sigma }}_{\nu }=\text{diag}({\gamma }_{N_{ROI}+1},\dots ,{\gamma }_{N_{ROI}+N_{\nu }})$ is the covariance of the artifactual sources. ${\mathbf{\gamma }}_k\in {\mathbb{R}}^{N_{ROI}+N_{\nu }}$ denotes a nonnegative scale vector that encodes the sparsity profile of the group of sources.

In general, the block structure of ${\mathbf{\Sigma }}_g$ could be learned in a data-driven way (Zhang & Rao, 2013), but we do not pursue that route here and assume it is given by the 68 regions of interest (ROI, $N_{ROI}=68$) covering the entire cortical surface of the Desikan-Killiany atlas (Desikan et al., 2006). In this setting, the matrices ${\mathbf{C}}_i\in {\mathbb{R}}^{N_i\times N_i}$ model the intra-ROI source covariances and can be precomputed based on source distances, taking into account the local folding of the cortex as described in Ojeda et al. (2018).

Next, we specify the priors on the hyperparameter $\lambda$ and $\mathbf{\gamma }$ (hyperpriors). Assuming that $\lambda$ and ${\gamma }_i$ are independent yields the factorization,

And since they are scale hyperparameters, we follow the popular choice of assuming gamma hyperpriors with noninformative scale and shape parameters on a log-scale of ${\lambda }^{-1}$ and ${\gamma }_i^{-1}$ (Tipping, 2001). This choice of hyperprior has the effect of assigning a high probability to low values of ${\gamma }_i$, which has been shown to shrink the irrelevant components of ${\mathbf{x}}_k$ to zero (Wipf & Nagarajan, 2008), thereby leading to a sparsifying behavior known as automatic relevance determination (ARD) (Neal, 1996; MacKay, 1992).

We use the Bayes theorem to express the posterior distribution of the sources as follows,

It is well known that a model with gaussian likelihood, $p\left({\mathbf{y}}_k\right|{\mathbf{x}}_k,\lambda )$, and prior, $p\left({\mathbf{x}}_k\right|\mathbf{\gamma })$, like ours, also has a gaussian posterior with the following mean and covariance,

where ${\mathbf{H}}^{†}$ represents the pseudoinverse of the observation operator $\mathbf{H}$ and $\mathbf{S}$ is the model data covariance,

2.3 . Sparse Model Learning

Once the measurements are gathered, the maximum a posteriori estimate of the sources and their covariance can be obtained analytically using equation 2.8. However, these formulas depend on the specific values taken by the hyperparameters that control our model, $\lambda$ and $\mathbf{\gamma }$. Doing a comprehensive survey of the algorithms available to learn these hyperparameters is beyond the scope of this letter, so in this section, we instead focus on the general idea and assumptions that go into this task.

The denominator of equation 2.7 is called the evidence and is a probabilistic measure of how much the data match our modeling assumptions. For a linear gaussian model like ours, it is readily expressed as follows,

Since we used noninformative hyperpriors (see section 2.2), we could optimize equation 2.10 with respect to $\lambda$ and $\mathbf{\gamma }$ to obtain source estimates conditioned on the optimal model. However, pointwise hyperparameter estimates tend to be noisy and computationally expensive. One way to obtain robust estimates is to learn these hyperparameters based on an ensemble of measurements (Cotter, Rao, & Kreutz-Delgado, 2005), $\mathbf{Y}=[{\mathbf{y}}_k,\dots ,{\mathbf{y}}_{k+n}]$. Furthermore, if we ignore the autocorrelation of the data and assume that each single measurement vector in $\mathbf{Y}$ is independent and identically distributed (i.i.d.), we can conveniently approximate the evidence of the ensemble as follows,

which is usually transformed into the so-called type II maximum likelihood (ML-II) cost function (Barber, 2012),

Equation 2.12 embodies a trade-off between model complexity and accuracy. The complexity term can be seen as the volume of an ellipsoid defined by $\mathbf{S}$. The volume of the ellipsoid is reduced by shrinking the axes that correspond to sources that are not needed to explain the data. The second term is a squared Mahalanobis distance that measures the accuracy of our model. Finally, the ML-II cost function is used to learn the optimal hyperparameters by solving the following optimization problem:

One of the key contributions of this letter is the realization that the augmented inverse problem of the M/EEG can be readily solved with existing off-the-shelf algorithms that tackle the optimization problem above (Faul & Tipping, 2002; Wipf & Nagarajan, 2008; Wipf et al., 2010; Zhang & Rao, 2013). In this letter, we solve equation 2.13 with the fast two-stage block-SBL algorithm proposed by Ojeda et al. (2018). Algorithm 1 shows the summary of SBL in pseudocode.

2.4 . Independent Component Analysis

In the analysis presented above, the matrix $\mathbf{H}$ is prespecified. In this section, we analyze the generative model of equation 2.3 from the ICA viewpoint. ICA seeks to estimate the source time series (called activations in the ICA literature) ${\mathbf{x}}_k$ from the data time series ${\mathbf{y}}_k$ without knowing the gain (mixing) matrix $\mathbf{H}$. Here we assume that the latent sources are instantaneously independent, which yields the following prior distribution,

To simplify the exposition, we assume the same number of sensors and sources, $N_y=N_x$; interested readers can find the case $N_y<N_x$ in Le, Karpenko, Ngiam, and Ng (2011) and Lewicki and Sejnowski (1998). From these premises, the objective of the algorithm is to learn the unmixing matrix ${\widehat{\mathbf{H}}}^{-1}$ such that we can estimate the sources with ${\widehat{\mathbf{x}}}_k={\widehat{\mathbf{H}}}^{-1}{\mathbf{y}}_k$. The unmixing matrix ${\widehat{\mathbf{H}}}^{-1}$ can be learned up to a permutation and rescaling factor, which has the inconvenience that the order of the learned components can change depending on the starting point of the algorithm and data. As in section 2.3, we can use a data block $\mathbf{Y}$ to write an approximation to the block likelihood function,

under the i.i.d. data assumption. We obtain each factor in equation 2.15 by integrating out the sources as follows,

As noted by MacKay (2008), assuming that the data are collected in the noiseless limit, $\lambda \to 0$, transforms the gaussian likelihood $p\left({\mathbf{y}}_k\right|{\mathbf{x}}_k,\mathbf{H},\lambda )$ into a Dirac delta function, in which case equation 2.16 leads to the Infomax algorithm of Bell and Sejnowski (1995). The learning algorithm essentially consists in finding the gradient of the log likelihood, $logp\left(\mathbf{Y}\right|\mathbf{H},\lambda )$, with respect to $\mathbf{H}$ and updating $\mathbf{H}$ on every iteration such that the likelihood of the data increases. As pointed out by Comon (1994), the ICA model is uniquely identifiable only if at most one component of ${\mathbf{x}}_{\mathbf{k}}$ is gaussian. Therefore, the prior densities $p_i\left(x_{i,k}\right)$ are usually assumed to exhibit heavier tails than the gaussian and, in particular, the prior $p_i\left(x_{i,k}\right)\propto {cosh}^{-1}{x}_{i,k}$ yields the popular contrast function $tanh\left({\mathbf{H}}^{-1}{\mathbf{y}}_k\right)$. It is worth noting that unlike the block-sparsity prior introduced in section 2.2, this one is motivated not by a biological consideration but by a mathematical necessity.

2.5 . Link Between Block-Sparse Source Imaging and ICA

In this section we show that ESI subject to the block-sparsity constraints laid out earlier necessarily yields groups of sources that are independent (or as independent as possible). We start by rewriting the biologically motivated source prior of equation 2.5 as follows,

where each factor is a gaussian pdf and i indexes a group of sources or an artifact component. To write equation 2.17 as the ICA prior in equation 2.14, we need to integrate out the hyperparameter ${\gamma }_i$ from each factor as follows:

Given our choice of hyperprior on ${\gamma }_i$ renders each marginalized prior $p_i\left({\mathbf{x}}_{i,k}\right)$ a heavy-tailed Student $t$-distribution (Tipping, 2001), which is known to induce independence. In other words, if we use the biologically inspired block-sparsity prior of equation 2.5, the ICA likelihood function of equation 2.16 becomes exactly the evidence in equation 2.10.

Note that we have not said anything very profound here, but we are merely pointing out that these two frameworks can be related in a principled way by choosing sparsity-inducing priors. The main difference between the two frameworks is that in SBL, we use biological information to constrain the $\mathbf{H}$ dictionary and source activity, while in ICA, both are learned from data.

2.6 . Recursive SBL

In this final theoretical section, we propose a simple modification to SBL to eliminate the i.i.d. data assumption while keeping the algorithm adaptive and computationally tractable. We start by introducing a spatiotemporal constraints in the form of a state (source) transition,

where the vector function $F$ models how the source activity evolves from one sample to the next and ${\mathbf{w}}_k$ is a perturbation vector. Several linear (Yang, Aminoff, Tarr, & Kass, 2016; Fukushima et al., 2012; Fukushima, Yamashita, Knösche, & Sato, 2015; Lamus, Hämäläinen, Temereanca, Brown, & Purdon, 2012; Cheung, Riedner, Tononi, & Van Veen, 2010; Galka, Yamashita, Ozaki, Biscay, & Valdés-Sosa, 2004) and nonlinear (Olier, Trujillo-Barreto, & El-Deredy, 2013; Giraldo, den Dekker, & Castellanos-Dominguez, 2010; Valdes-Sosa, Sanchez-Bornot et al., 2009; Daunizeau, Kiebel, & Friston, 2009) models have been proposed to express the neural dynamics; for simplicity, here we use the time-invariant linear model proposed by Yamashita, Galka, Ozaki, Biscay, and Valdes-Sosa (2004). Yamashita's model considers nearby source interactions only and describes the evolution of the ith source $g_{i,k}$ by the following first-order autoregressive model,

where the constants $\alpha$ and $\beta$ are set so that the system is observable (Galka et al., 2004), $\mathcal{N}\left(i\right)$ contains the indices of the direct neighbors of dipole, $i$ and $n_i$ is its total number of neighbors. Dipole neighborhoods can be extracted from the tessellation of the cortical surface. In the absence of any other obvious transition model for artifact sources, we propose a simple random walk, which together with equation 2.19 yields the following linear transition function,

where the constant damping factor $\zeta =0.99$ has the job of stabilizing the random walk.

The state perturbation process ${\mathbf{w}}_k$ encompasses modeling errors as well as random inputs coming from distant sources, which we assume to be gaussian and serially uncorrelated, ${\mathbf{w}}_k\sim N\left({\mathbf{w}}_k\right|\mathbf{0},{\mathbf{Q}}_{\mathbf{k}})$. These assumptions yield the following conditional source prior,

where ${\mathbf{P}}_{k-1}$ is the state covariance at $k-1$ and we assume $p\left({\mathbf{x}}_0\right|{\mathbf{x}}_{-1},\mathbf{\gamma })=N\left({\mathbf{x}}_0\right|\mathbf{0},\mathbf{Q})$.

Swapping out the static prior of equation 2.5 by the one in equation 2.20 transforms our generative model into a linear gaussian dynamic systems (LGDS). In this type of system, the source time series can be estimated optimally using the Kalman filter (Kalman, 1960). We first remove the dependence on the previous state, ${\mathbf{x}}_{k-1}$, by computing the predictive density,

where $p\left({\mathbf{x}}_{k-1}\right|{\mathbf{y}}_{k-1},\lambda ,\mathbf{\gamma })$ is the source posterior already determined in the previous time step, $k-1$. In an LGDS, the predictive density is also gaussian with the following mean and covariance (Ozaki, 2012),

where we use the notation $n|m$ (with $n\ge m$) to indicate quantities estimated in the step $n$ using data up to the $m$ sample. Equation 2.21 is known as the time update. The posterior is also a gaussian, so we readily write down its mean and covariance,

where $\mathbf{K}={\widehat{\mathbf{P}}}_{k|k-1}{\mathbf{H}}^T{\mathbf{S}}^{-1}$ represents the Kalman gain, the analytical data covariance is $\mathbf{S}=\widehat{\lambda }{\mathbf{I}}_{N_y}+\mathbf{H}{\widehat{\mathbf{P}}}_{k|k-1}{\mathbf{H}}^T$, and ${\tilde{\mathbf{y}}}_k$ is the innovation (one-step ahead prediction error). Equation 2.22 is known as the measurement update, and plugging equation 2.21 into it yields the following recursive formula,

In our case, the source space is usually very large (thousands or tens of thousands of sources), making the direct use of the Kalman filter impractical due to the necessity to estimate the $N_x\times N_x$ state covariance matrix ${\mathbf{P}}_k$. To make the filter tractable, Yamashita et al. (2004) proposed ignoring the contribution of the neurodynamic model in the evolution of the state covariance, $\mathbf{F}{\widehat{\mathbf{P}}}_{k-1|k-1}{\mathbf{F}}^T\to 0$, in equation 2.21. With this approximation, the Kalman gain becomes the pseudoinverse ${\mathbf{H}}^{†}$ in equation 2.8 and the model data covariance is calculated as defined in equation 2.9, which we plug into the main recursive formula of the Kalman filter in equation 2.23. We continue estimating the hyperparameters $\lambda$ and $\mathbf{\gamma }$ through the SBL algorithm outlined in section 2.3 with the difference that now the evidence is a function of the innovations:

We call this new approach recursive SBL (RSBL) and we summarize it in algorithm 2.

It is important to note that unlike most common use cases of the Kalman filter, in ESI, we do not have access to a perfect description of the brain dynamics, let alone a computationally tractable one. Thus, we usually rely on approximate neurodynamic models, of which equation 2.19 is an example. We also remark that our model yields a transition function that provides a simple local spatiotemporal filtering effect that only reflects nearby source interactions. However, the actual transition may be far more complicated. Since modeling errors are unavoidable, we argue that it makes sense from a practical standpoint to downplay the role of the neurodynamic model in the update of the state covariance. To compensate for this additional approximation, we exploit the block structure of the state vector induced by the hyperparameter vector $\mathbf{\gamma }$ to constrain the correction of the measurement update to regions of the state-space with biological significance.

Our approach differs from the one proposed by Yamashita et al. (2004) in two ways. First, they regularize the same amount everywhere in the source space (equivalent to considering a single $\gamma$ parameter for the whole cortex), while we do so in groups of nearby sources encouraging sparsity in the ROI domain. Second, they do not reestimate the hyperparameters and we do.

2.7 . Computation of Dictionaries L and A

As we mentioned in section 2.1, a database of stereotypical artifact scalp projections can be obtained by running ICA on a large collection of EEG data sets. The database could include data from different studies, populations, and montages so that a “universal” artifact dictionary can be compiled offline. Note that the key idea is to use the precomputed artifact dictionary to approximate the scalp projection of stereotypical artifact components of new subjects without actually running ICA on their data. Toward that end, ideally, each subject would have its companion structural MRI and digitized channel locations fully describing the anatomical support on which EEG data were collected.

If only sensor positions were available, the procedure proposed by Darvas, Ermer, Mosher, and Leahy (2006) could be used to obtain individualized head models using a single template head. Most EEG studies, however, do not include MRI data or subject-specific sensor positions, just the sensor labels. Here we propose a co-registration procedure that requires sensor positions only, either measured with a digitizer or pulled from a standard montage file, and a template head model. We use a four-layer (scalp, outer skull, inner skull, and cortex) head model derived from the Colin27 MRI template with fiducial landmarks (nasion, inion, vertex, left and right preauricular points) and 339 sensors located on the surface of the scalp. The sensors are placed and named according to a superset of the 10/20 system. We define the brain source space as 8003 orientation-free dipoles scattered uniformly over the template's cortical layer, resulting in $N_g=8003\times 3$ current source elements (8003 dipoles and their orientations along $x$-, $y$-, and $z$-axes).

Although several algorithms can be used to automatically identify artifact scalp projections derived from ICA (Pion-Tonachini et al., 2017; Winkler, Haufe, & Tangermann, 2011), as proof of principle, here we rely on the expertise of the authors. Since inspecting each IC in a large database is a cumbersome task, we propose the following simplification. We first co-register each subject in the database with the head model template to map ICs defined on different sensor spaces to a common space. We collect the mapped ICs in a matrix of 339 (number of channels of the template) by the number of total ICs. We reduce dimensionality by applying the k-means algorithm to the columns of the IC matrix. After inspecting the cluster centroids, we label them as Brain, EOG, EMG, or Unknown (cluster of scalp maps of unknown origin). Then we store the indices of the EOG and EMG cluster centroid nearest neighbors for further use in the creation of the artifact dictionary $\mathbf{A}$.

We use the fiducial points marked on the participant's head at the start of data collection to estimate an affine transformation from the individual space to the space of the template. If fiducial points are not available, a common set of sensors between the template and the individual montage can be used to estimate the transformation. We use the affine transformation to map the montage of the participant onto the skin surface of the template; this is what we call in Figure 2 “individualized head model.” We compute the lead field matrix $\mathbf{L}$ using the boundary element method implemented in the OpenMEEG package (Gramfort, Papadopoulo, Olivi, & Clerc, 2010). Although this procedure can incur significant errors due to differences in subject and template head shapes and errors in the sensor placement, in many cases, this is the best one can do given the lack of accurate spatial information. A promising alternative to explore is the co-registration procedure recently proposed by Duque-Muñoz et al. (2019) for optically pumped magnetometer sensors in which the sensor placement orientation is learned from data.

To produce an individualized $\mathbf{A}$ matrix, we linearly warp the EOG and EMG cluster centroids to the space of the individualized head model so that each artifact scalp projection has the same column length as L and can be appended to it. Finally, we divide each column of $\mathbf{H}$ by its norm so that the relative contribution to each source is determined by the amplitude of the source activation vector ${\mathbf{x}}_k$ only.

2.8 . Data

To construct the artifact dictionary, we used data from two different studies made public under the umbrella of the BNCI Horizon 2020 project² (Brunner et al., 2015). The data sets were selected due to their well-curated metadata and their relatively high number of sensors.

3.1 . Performance on Simulated Data

The objective of this experiment was to unmix the EEG signal and recover all simulated brain and artifact sources under different common-mode noise conditions. To that end, we simulated 2 seconds of evoked EEG data contaminated by lateral and vertical EOG activity (see Figure 4). By construction, we only simulated brain activity in the anterior (ACC) and posterior (PCC) cingulate cortex and designed their time courses to mimic the distribution of scalp voltages observed during ErrP studies (Gehring, Liu, Orr, & Carp, 2011). ErrPs are characterized by a negativity (Ne component) in the frontocentral channels within 100 msec after an erroneous response, followed by a positivity (Pe component) in the centroparietal channels around 250 msec post-error. The timing of these components can be reliably identified by inspecting the EEG activity of frontocentral channels (Fz, FCz, Cz) time-locked to the erroneous responses.

We simulated Ne and Pe spatial profiles with the column vectors ${\mathbf{g}}_{\text{A}}$ and ${\mathbf{g}}_{\text{P}}$, each with $8003\times 3$ elements. ${\mathbf{g}}_{\text{A}}$ had 294 sources (98 dipoles $\times$ 3 orientations) in the ACC set to one and zero elsewhere, and ${\mathbf{g}}_{\text{P}}$ had 390 sources (130 dipoles $\times$ 3 orientations) in the PCC set to one and zero elsewhere. The temporal profiles, ${\mathbf{t}}_{Ne}$ and ${\mathbf{t}}_{Pe}$ were simulated using gamma functions. We simulated the cortical activity by multiplying the spatial and temporal profiles as follows,

The two artifact sources simulated a lateral eye movement (from $-$750 msec to 250 msec) and eye blink (from 0 msec to 400 msec) events. All other artifact sources were set to zero. The ratios between the maximum amplitude of artifacts, ${\text{EOG}}_{\mathrm{h}}$ and ${\text{EOG}}_{\mathrm{v}}$, to cortical source activity were set to 5 and 10, respectively. The simulated source activity was generated by concatenating the cortical ($\mathbf{G}$) and artifact source matrices. We simulated the sensor space with a 64-channel montage placed according to the extended 10/20 system and computed $\mathbf{H}$ as described in section 2.7. Scalp data were generated by projecting the sources onto the sensor space through $\mathbf{H}$ and adding white gaussian noise. The variance of the noise was calculated with respect to the variance of the sensor data before adding EOG artifacts.

Table 1 summarizes the performance of the RSBL algorithm for different levels of noise. Column 1 shows the SNR in dB units, and column 2 shows the equivalent EEG signal to noise amplitude ratio (SNAR). The columns headed by ACC, PCC, ${\text{EOG}}_{\mathrm{v}}$, and ${\text{EOG}}_{\mathrm{h}}$ report the correlation value between their respective simulated and estimated sources; high values indicate good estimation accuracy. The columns headed by $\setminus$ACC and $\setminus$PCC report the mean correlation between simulated ACC and PCC sources and every other estimated cortical source; low correlation values indicate low source leakage. We note that for extremely aggressive noise conditions (i.e., SNR values between $-$10 dB and 0 dB), the algorithm failed to consistently reconstruct all the simulated sources accurately. However, once the amplitude of the EEG was approximately one and a half higher than the amplitude of the noise, we obtained correlations in the 0.63 to 0.99 range. In all cases, the leakage was low, as indicated by correlation values below 0.04.

Figure 4 illustrates the data and sources for a simulation with SNR of 6 dB. The traces in panel A represent the raw, cleaned, and simulated EEG time series for a subset of channels. The cleaned EEG traces were obtained by projecting out the estimated EOG source activity. The panels on the right show the ground truth and estimated source activity in the ACC (B), PCC (C) areas, as well as EOG artifacts (D). Note that the estimated source time series are sparse in space and time (see also Figure 5). In particular, panels B and C demonstrate that source values that randomly oscillate at the noise level are ignored by the algorithm. That does not mean that those cortical areas are silent, but that the postsynaptic potentials produced by the collective firing of pyramidal cells in those regions are not coherent enough to create signals that can be reliably measured by the scalp sensors.

Figure 5 shows a snapshot of the scalp and source maps at the peak of the Ne and Pe components displayed in Figure 4. We linearly extrapolated the sensor values to the portion of head covered by the EEG cap, so that the topographic maps could be rendered on the 3D surface. In the top row, we see that although all the simulated sources within an ROI have the same amplitude, in the recovered maps (bottom row), the nontrivial source values are not identical (though correlated). The panels in the bottom row show the EEG scalp maps with the artifact activity projected out. Note that even when the scalp data are severely affected by the eye blink artifact at the peak of the Ne component (t $=$ 44 msec), the algorithm correctly estimated the orientation and location of its generators in the ACC.

3.2 . Source Reconstruction of Ongoing Real Error-Related EEG Data

In this experiment we used real data to demonstrate the performance of the RSBL algorithm. To that end, we selected a subject from the testing set that belonged to the ErrP study (data set 1 above). In the experimental protocol used to collect those data, the subjects were tasked with moving a cursor toward a target location either using a keyboard or mental commands (Chavarriaga & del R. Millán, 2010). The authors of that paper found consistent evoked potentials produced by errors induced by the computer. In other words, subjects elicited error potentials after watching the computer execute the wrong moves; this is sometimes referred to as feedback error–related negativity/positivity (Ward et al., 2013). These error potentials were characterized by two frontocentral positive peaks around 200 msec and 320 msec after the feedback, a frontocentral negativity around 250 msec, and a broader frontocentral negative deflection about 450 msec.

Figure 6 summarizes the results of applying the RSBL algorithm to one trial selected at random. We checked beforehand that the selected trial wasn't obviously affected by artifacts. Panel A shows the EEG signal of three frontocentral sensors (Fz, FCz, and Cz). Panel B shows the source magnitude time series averaged within the ACC and PCC regions, panel C shows the scalp and cortical maps at the latency of the Ne and Pe components, and panel D shows the maximum intensity projection of each source map in panel C in the sagittal plane. Panels B and C reveal that the error components observed at the scalp level are mostly generated by a complicated interplay of sources in the ACC and PCC regions. However, the sagittal maximum intensity projections indicate that other cortical sources also contribute to the observed EEG signal, although a lesser amount. This is common in single-trial analyses, and usually the sources that are unrelated to the experimental condition being studied clear out after trial averaging. We note that the result that both ACC and PCC sources contribute to Ne and Pe scalp topographies is in agreement with other imaging studies (Buzzell et al., 2017; Ojeda et al., 2018).

3.3 . EOG Artifact Removal on Real Data

In this example, we applied the RSBL algorithm to a trial contaminated by a lateral eye movement and eye blink activity. We used data from the same subject selected for the experiment in the previous section, but this time we analyzed an epoch with no error-related activity so that we could appreciate the artifact rejection performance of the algorithm minimizing task-related confounds. Figure 7 summarizes the data and results.

Figure 7A shows a subset of the raw and reconstructed (cleaned) EEG traces in black and red, respectively. There is a lateral eye movement artifact between $-$1250 msec and $-$800 msec and an eye blink between $-$250 msec and 1000 msec. Panel B shows the estimated EOG${}_v$ and EOG${}_h$ artifact source activity in blue and orange, respectively. We note that these artifact sources were active at the latencies where the EEG is affected and mostly zero elsewhere. This is a desired feature of the algorithms because this way, it only “corrects” the affected segments while leaving clean data unchanged. In panel C, we show the maximum intensity projection maps of the source map at the latency of the peak eye blink artifact, 0 msec. Each column displays a different projection. The rows display source estimates without and with artifact modeling enabled. Panel D shows the time series of the log evidence for generative models with artifact modeling on and off. In panel D, both traces differ mostly only when artifacts occur, where higher log evidence of the blue trace indicates that source estimation benefits from modeling artifacts. This is not a striking finding, but it illustrates the practical utility of the evidence metric for data modeling. Note that although in our implementation, we optimize the ML-II cost function (see equation 2.12), here we put this value back in log evidence units to facilitate its interpretation (i.e., the higher it is, the better our model captures the data).

Visual inspection of the bottom row of panel C reveals that some residual eye blink activity may have leaked into the frontal pole. In practice, it is extremely hard to totally remove artifactual activity because (1) the use of a lead field matrix derived from a template head model may misfit the anatomy of the subject, introducing errors in the $\mathbf{L}$ dictionary, (2) errors in the sensor locations can cause the EEG topography to shift with respect to the expected brain and artifact source projections, (3) EMG scalp projections are difficult to characterize due to their variability, as opposed to EOG projections that are more stereotyped, and (4) unmodeled artifactual activity, such as muscle projections toward the back of the head that were largely ignored in this study, may be suboptimally accounted for. Despite all these issues, Figure 7 demonstrates that RSBL can yield reasonably robust source estimates in the presence of high-amplitude artifacts. Furthermore, panel D suggests that we could use dips in the log evidence to inform subsequent processing stages of artifactual events that were not successfully dealt with.