Electromyogram (EMG) Removal by Adding Sources of EMG (ERASE)—A Novel ICA-Based Algorithm for Removing Myoelectric Artifacts From EEG

Li, Yongcheng; Wang, Po T.; Vaidya, Mukta P.; Flint, Robert D.; Liu, Charles Y.; Slutzky, Marc W.; Do, An H.

doi:10.3389/fnins.2020.597941

ORIGINAL RESEARCH article

Front. Neurosci., 15 January 2021

Sec. Neural Technology

Volume 14 - 2020 | https://doi.org/10.3389/fnins.2020.597941

Electromyogram (EMG) Removal by Adding Sources of EMG (ERASE)—A Novel ICA-Based Algorithm for Removing Myoelectric Artifacts From EEG

$\nYongcheng Li$ Yongcheng Li¹^*

Po T. Wang²

Mukta P. Vaidya^3,4,5

Robert D. Flint^3,4,5

Charles Y. Liu^6,7,8

Marc W. Slutzky^3,4,5

An H. Do¹^*

¹Department of Neurology, University of California, Irvine, Irvine, CA, United States
²Department of Biomedical Engineering, University of California, Irvine, Irvine, CA, United States
³Department of Neurology, Northwestern University, Chicago, IL, United States
⁴Department of Physiology, Northwestern University, Chicago, IL, United States
⁵Department of Physical Medicine and Rehabilitation, Northwestern University, Chicago, IL, United States
⁶Department of Neurosurgery, University of Southern California, Los Angeles, CA, United States
⁷Rancho Los Amigos National Rehabilitation Center, Downey, CA, United States
⁸Neurorestoration Center, University of Southern California, Los Angeles, CA, United States

Electroencephalographic (EEG) recordings are often contaminated by electromyographic (EMG) artifacts, especially when recording during movement. Existing methods to remove EMG artifacts include independent component analysis (ICA), and other high-order statistical methods. However, these methods can not effectively remove most of EMG artifacts. Here, we proposed a modified ICA model for EMG artifacts removal in the EEG, which is called EMG Removal by Adding Sources of EMG (ERASE). In this new approach, additional channels of real EMG from neck and head muscles (reference artifacts) were added as inputs to ICA in order to “force” the most power from EMG artifacts into a few independent components (ICs). The ICs containing EMG artifacts (the “artifact ICs”) were identified and rejected using an automated procedure. ERASE was validated first using both simulated and experimentally-recorded EEG and EMG. Simulation results showed ERASE removed EMG artifacts from EEG significantly more effectively than conventional ICA. Also, it had a low false positive rate and high sensitivity. Subsequently, EEG was collected from 8 healthy participants while they moved their hands to test the realistic efficacy of this approach. Results showed that ERASE successfully removed EMG artifacts (on average, about 75% of EMG artifacts were removed when using real EMGs as reference artifacts) while preserving the expected EEG features related to movement. We also tested the ERASE procedure using simulated EMGs as reference artifacts (about 63% of EMG artifacts removed). Compared to conventional ICA, ERASE removed on average 26% more EMG artifacts from EEG. These findings suggest that ERASE can achieve significant separation of EEG signal and EMG artifacts without a loss of the underlying EEG features. These results indicate that using additional real or simulated EMG sources can increase the effectiveness of ICA in removing EMG artifacts from EEG. Combined with automated artifact IC rejection, ERASE also minimizes potential user bias. Future work will focus on improving ERASE so that it can also be used in real-time applications.

1. Introduction

Electroencephalographic (EEG) signals are often contaminated by surface electromyographic (EMG) signals and can make it difficult to appropriately interpret EEG signals or use them in various neuroengineering applications. There are many approaches to remove EMG artifacts from EEG signals, which can be broadly classified in the following categories (Fatourechi et al., 2007; Muthukumaraswamy, 2013; Urigüen and Garcia-Zapirain, 2015; Minguillon et al., 2017): (1) filtering (e.g., adaptive filtering, Wiener filtering, and Bayes filtering; Morbidi et al., 2007; Guerrero-Mosquera and Vazquez, 2009; Sweeney et al., 2012b), (2) linear regression method (Gratton et al., 1983; Gwin et al., 2010; De Vos et al., 2014), (3) blind source separation (BSS) methods which includes principal component analysis (PCA) (Berg and Scherg, 1991; Ille et al., 2002), independent component analysis (ICA) (Comon, 1994; James and Hesse, 2004; Nolan et al., 2010; Mognon et al., 2011), canonical correlation analysis (CCA) (Safieddine et al., 2012; Sweeney et al., 2012a; Chen et al., 2015; Mowla et al., 2015), sparse component analysis (SCA) (Gribonval and Lesage, 2006), singular spectrum analysis (SSA) (Teixeira et al., 2006), and etc. (4) source decomposition, including discrete wavelet transform (DWT) (Donoho and Johnstone, 1995; Unser and Aldroubi, 1996; Safieddine et al., 2012; Peng et al., 2013), empirical mode decomposition (EMD) (Mijovic et al., 2010; Mourad and Niazy, 2013; Zeng et al., 2013), stationary wavelet transform (SWT) (Zikov et al., 2002; Krishnaveni et al., 2006), and ensemble empirical mode decomposition (EEMD) (Wu and Huang, 2009; Sweeney et al., 2012a; Chen et al., 2014a,b), (5) neural networks (NN) (Nguyen et al., 2012; Burger and van den Heever, 2015) and adaptive neural fuzzy inference systems (ANFIS) (Hu et al., 2015; Suja Priyadharsini et al., 2016). In the setting of a plethora of EMG artifacts removal algorithms, previous surveys (Muthukumaraswamy, 2013; Urigüen and Garcia-Zapirain, 2015) suggest that BSS methods are the most commonly used and outperform other approaches at EMG artifacts removal from EEG. Among the BSS methods mentioned above, Urigüen and Garcia-Zapirain (2015) states that 45% of work in their bibliography use ICA to remove the EMG artifact from EEG. CCA is another widely employed method for the removal of EMG artifacts from EEG in recent years (Chen et al., 2015, 2017). However, it does not outperform ICA at removing EMG artifacts from EEG (McMenamin et al., 2010; Escudero et al., 2011; Safieddine et al., 2012) or at removing ocular and cardiac artifacts (Romero et al., 2008, 2009; Pham et al., 2011; Delorme et al., 2012; Evans et al., 2012; Daly et al., 2013).

Despite the popularity of ICA for EMG artifact rejection in EEG, its use is still affected by several issues (Tran et al., 2004; Delorme et al., 2007; Muthukumaraswamy, 2013). For example, nearly all EEG channels are typically contaminated by EMG, and there is a high spatiotemporal overlap between EMG artifacts and EEG signal (Shackman et al., 2009). Therefore, conventional ICA algorithms are usually unable to separate EMG artifacts from the EEG signal—that is, it is difficult to “force” all of the EMG artifacts into an isolated set of independent components. Hence, post-ICA-treated data may still include residual EMG (Shackman et al., 2009; Olbrich et al., 2011). Also, since the EMGs could have multiple source locations and show a large overlap with the higher frequency components (>20 Hz) of EEG signals, it is difficult to assign a universal operational definition for EMG components (Mammone et al., 2012; Gross et al., 2013). Hence, rejection of EMG artifact components is typically performed manually when no prior knowledge about artifacts is available. This leads to potential over- or under-rejection of components as users attempt to distinguish between neurogenic and myogenic components in common ICA. Further, rejecting components manually is cumbersome/time-consuming and can introduce subjectivity. Additionally, accurate extraction of source signals in the ICA model is another issue since the global optimum in these algorithms is typically affected by the contrast function. Some approaches, which are presented in further detail below, were developed to solve these issues.

Some studies demonstrated that prior knowledge about artifacts or source signals can improve the effectiveness of ICA at removing artifacts (Akhtar et al., 2012; Urigüen and Garcia-Zapirain, 2015). Therefore, constrained ICA (cICA) or ICA with reference (ICA-R) which incorporate prior knowledge about the artifact and/or source signals were developed (James and Hesse, 2004). This is performed by imposing temporal or spatial constraints on the source mixture model. In temporally constrained ICA, prior knowledge about artifacts can be introduced into the ICA model to identify the artifact IC/ICs by solving a constrained optimization problem (James and Gibson, 2003; James and Hesse, 2004; Lu and Rajapakse, 2005; Lin et al., 2007; Romero et al., 2008). Temporally constrained ICA can only be used for the removal of EOG and ECG, but is not highly effective for EMG artifacts (James and Gibson, 2003; Lu and Rajapakse, 2005; Romero et al., 2008). EMG artifacts sources, which are more time-varying and non-stationarity, are too complicated to characterize for optimization constraints. Hence, no prior studies have adequately addressed the separation of spatiotemporal overlap between EMG and EEG by using temporally constrained ICA. Spatially constrained ICA incorporates prior knowledge or assumption of spatial topographies of some source projections acting as a spatial filter, and limit the degree to which some of the columns of the mixing matrix may deviate from the known projections (Ille et al., 2001, 2002; Hesse and James, 2006; Akhtar et al., 2012). As mentioned above, EMG artifacts are time-varying and highly overlapped with EEG. Hence, known spatial topographies of EMG source projections derived from previous data recording or mathematically simulated model are usually inaccurate and difficult to achieve. Although some work on EMG artifact removal utilized spatially constrained ICA with EEG during seizures (Hesse and James, 2006), the performance of spatially constrained ICA on these ictal EEG, where the ground truth of EMG artifacts is actually unknown, has not been fully and rigorously established. Hence, using spatially constrained ICA to remove the EMG artifacts from real EEG is still unsubstantiated. Additionally, running ICA iteratively is required for both types of constrained ICA to reject the artifacts, which is time-consuming and hard to achieve real-time application. Even so, when the prior knowledge about artifacts or source signal is available, some form of spatially constrained ICA is preferred as compared to using ICA alone (Akhtar et al., 2012; Urigüen and Garcia-Zapirain, 2015).

Another important issue associated with removing EMG artifacts from EEG via ICA is automated EMG artifact rejection. Generally, when reference waveforms are available, there is one method to achieve automated rejection based on ICA in previous studies (Urigüen and Garcia-Zapirain, 2015). Specifically, cICA compares spatial and temporal statistical characteristics of ICs to those from background EEG or artifacts. Subsequently, a combination of thresholds for those statistical characteristics is used in an automated algorithm to reject the artifacts (Delorme et al., 2001, 2007; Nolan et al., 2010; Mognon et al., 2011; Daly et al., 2012). This method has been demonstrated to reliably remove EOG and ECG artifacts, particularly since these signals drastically differ from EEG in spatial and temporal statistical characteristics (Wallstrom et al., 2004; Romero et al., 2008, 2009). However, using this method for EMG artifact removal is still inadequate due to the spatiotemporal overlap between EMG artifacts and EEG signals (i.e., EMG artifacts always introduce a large number of unique scalp maps, leaving few ICs available for capturing brain sources). Therefore, it is necessary to develop an automated technique that can more effectively and systematically remove EMG artifacts while not affecting any of the underlying signal features.

In order to improve the effectiveness and reliability of ICA in removing the EMG artifacts from EEG and establish an effective automated artifacts rejection procedure, we introduce a novel method, termed as EMG Removal by Adding Sources of EMG (ERASE) and subject it to rigorous mathematical and experimental validation. ERASE combines the advantages of two types of cICA. It aims to improve upon ICA by adding either real or simulated EMG artifacts as extra “reference” channel signals into the EEG data. We mathematically demonstrated that if the reference EMG artifacts were not independent of the contaminant EMG artifacts in EEG, a larger proportion of the artifacts could be identified by several specific independent components (ICs) after running ICA. Also, we proposed criteria based on the mixing matrix to automatically identify and reject the artifacts components. Our results revealed that ERASE had higher effectiveness in removing the EMG artifacts compared to conventional ICA. In summary, this study developed an effective EMG rejection approach, which can provide more confidence for the utilization of EEG in applications such as physiological studies underlying motor behaviors.

2. Methods

2.1. Description of ERASE ICA Model Based on Added EMG Sources

2.1.1. Model Description

To facilitate more effective removal of EMG artifacts from the EEG data, we combined EMG artifacts (here, either simulated EMG or recorded EMG) with EEG datasets and applied a modified ICA model as follows:

\begin{array}{l} (\begin{matrix} {\hat{X}}_{c_{E}} \\ n_{c_{M}}^{*} \end{matrix}) = A_{c_{E} + c_{M}} \times S_{c_{E} + c_{M}} \\ (\begin{matrix} X_{c_{E}} + b_{c_{E}} \cdot N_{c_{E}} \\ n_{c_{M}}^{*} \end{matrix}) = A_{c_{E} + c_{M}} \times (\begin{matrix} s_{c_{E}} + m_{c_{E}} \\ m_{c_{M}}^{*} \end{matrix}) & (1) \end{array}

where c_E is the number of the EEG channels (c_E dimension); c_M is the number of the reference EMG channels (c_M dimension), ${\hat{X}}_{c_{E}} = X_{c_{E}} + b_{c_{E}} \cdot N_{c_{E}}$ , and X_{c_E} is the uncontaminated EEG data; N_{c_E} is the contaminant EMG artifacts, which usually is the real EMG artifacts inside of EEG; b_{c_E} is the linear coefficients; $n_{c_{M}}^{*}$ is the reference EMG artifacts, which are extra channels containing EMG signals from muscles or simulation; A_{c_E+c_M} is the mixing matrix of dimension (c_E+c_M) × (c_E+c_M); S_{c_E+c_M} is the independent component sources with c_E+c_M dimension, in which s_{c_E} is the sources representing the uncontaminated EEG; m_{c_E} are the sources representing the contaminant EMG artifacts, and the use of m_{c_E} in the formula can help to differ our model from conventional ICA model and explain the conditions where our model is true; $m_{c_{M}}^{*}$ are the reference EMG sources.

Theorem: Given that the reference EMG sources are independent, and $n_{c_{M}}^{*}$ is only dependent with N_{c_E}, then m_{c_E} = 0.

Proof: We considered N_{c_E}, $n_{c_{M}}^{*}$ , m_{c_E}, and $m_{c_{M}}^{*}$ as stochastic variables, and assume that m_{c_E}≠0. In this model, $N_{c_{E}} = a 1 \times m_{c_{E}} + a 2 \times m_{c_{M}}^{*}$ and $n_{c_{M}}^{*} = a 3 \times m_{c_{E}} + a 4 \times m_{c_{M}}^{*}$ (a1, a2, a3 and a4 are the corresponding submatrices of the mixing matrix, of which dimensions are c_E×c_E, c_E×c_M, c_M×c_E, and c_M×c_M, respectively). Since N_{c_E} is dependent with $n_{c_{M}}^{*}$ , dependence between N_{c_E} and $n_{c_{M}}^{*}$ can be defined as follows:

\begin{array}{l} f_{N_{c_{E}} n_{c_{M}}^{*}} (N_{c_{E}}, n_{c_{M}}^{*}) \neq f_{N_{c_{E}}} (N_{c_{E}}) \cdot f_{n_{c_{M}}^{*}} (n_{c_{M}}^{*}) & (2) \end{array}

where × denotes matrix multiplication, · denotes element-wise multiplication, f_{N_{_c_E}} is the probability density function for N_{c_E}, $f_{n_{c_{M}}^{*}}$ is the probability density function for $n_{c_{M}}^{*}$ , and $f_{N_{c_{E}} n_{c_{M}}^{*}}$ is the joint probability density function of N_{c_E} and $n_{c_{M}}^{*}$ .

This expression of N_{c_E} and $n_{c_{M}}^{*}$ can be combined as follows:

\begin{array}{l} f_{N_{c_{E}} n_{c_{M}}^{*}} ((a 1 \times m_{c_{E}} + a 2 \times m_{c_{M}}^{*}), (a 3 \times m_{c_{E}} + a 4 \times m_{c_{M}}^{*})) \\ \neq f_{N_{c_{E}}} (a 1 \times m_{c_{E}} + a 2 \times m_{c_{M}}^{*}) \cdot f_{n_{c_{M}}^{*}} (a 3 \times m_{c_{E}} + a 4 \times m_{c_{M}}^{*}) & (3) \end{array}

We changed the expression above with m_{c_E} and $m_{c_{M}}^{*}$ :

\begin{array}{l} f_{m_{c_{E}} m_{c_{M}}^{*}} (m_{c_{E}}, m_{c_{M}}^{*}) \neq f_{m_{c_{E}}} (m_{c_{E}}) \cdot f_{m_{c_{M}}^{*}} (m_{c_{M}}^{*}) & (4) \end{array}

where f_{m_{_c_E}} is the probability density function of m_{c_E}, $f_{m_{c_{M}}^{*}}$ is the probability density function of $m_{c_{M}}^{*}$ , and $f_{m_{c_{E}} m_{c_{M}}^{*}}$ is the joint probability density function of m_{c_E} and $m_{c_{M}}^{*}$ . From the Equation (4), we know that m_{c_E} is dependent with $m_{c_{M}}^{*}$ . This represents reference EMG sources being dependent with sources representing contaminant EMG artifacts. Since this violates the ICA principle of component independence, m_{c_E} must equal 0.

Two key assumptions are made in our model. The first assumption is the independence among the reference EMG sources and the second assumption is the dependence between the contaminant EMG artifacts and the reference EMG artifacts. This means if the simulated EMG can meet these two assumptions, in situations where real EMG was not collected, or the situation does not allow for EMG recordings, simulated EMG can also act as the reference EMG artifacts.

2.1.2. Rejection Criteria

We define the (c_E+c_M) × (c_E+c_M) mixing matrix A_{c_E+c_M} in equation (1) as below:

\begin{array}{l} A_{c_{E} + c_{M}} = (\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, c_{E} + c_{M}} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, c_{E} + c_{M}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{c_{E}, 1} & a_{c_{E}, 2} & \dots & a_{c_{E}, c_{E} + c_{M}} \\ a_{c_{E} + 1, 1} & a_{c_{E} + 1, 2} & \dots & a_{c_{E} + 1, c_{E} + c_{M}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{c_{E} + c_{M}, 1} & a_{c_{E} + c_{M}, 2} & \dots & a_{c_{E} + c_{M}, c_{E} + c_{M}} \end{matrix}) & (5) \end{array}

where the first c_E rows are the coefficients corresponding to the EEG channels, the last c_M rows are the coefficients corresponding to the reference EMG channels, the coefficients in each column represent the weights of this IC to all EEG/EMG channels.

To develop an automatic method of identifying and rejecting the independent components related to EMG artifacts (referred to as “artifact ICs”) in real EEG data after running ICA, two criteria were defined:

• First, a threshold was established based on the root mean square (RMS) values of coefficients in the mixing matrix rows corresponding to the EMG channels, which was defined as:

\begin{array}{l} R_{m s} = \frac{\sqrt{\frac{\sum_{n = 1}^{c_{E} + c_{M}} a_{c_{E} + 1, n}^{2}}{c_{E} + c_{M}}} + \dots + \sqrt{\frac{\sum_{n = 1}^{c_{E} + c_{M}} a_{c_{E} + c_{M}, n}^{2}}{c_{E} + c_{M}}}}{c_{M}} & (6) \end{array}

where R_ms is the average RMS value across all the reference EMG artifact channels, a_{c_E+1, n} is the n^th coefficient in ${(c_{E} + 1)}^{t h}$ row in the mixing matrix, and a_{c_E+c_M, n} is the n^th coefficient in ${(c_{E} + c_{M})}^{t h}$ row in the mixing matrix. The true threshold was calculated by RMS value times gain. Gain is a constant which was empirically set between 0.4 and 3. The ICs whose absolute value of coefficients in the corresponding EMG rows were above the threshold were defined to be artifact ICs.

• Second, the ICs whose maximal absolute value of coefficients (max(|a_{1, i}|, …, |a_{c_E+c_M, i}|), 1 ≤ i ≤ (c_E+c_M)) corresponds to a hat band electrode, were rejected. Note that hat band electrodes were defined as all the EEG electrodes which were on the outermost circumference of the head, as defined by Yamada and Meng (2012) (Supplementary Figure 8).

In order to find the proper threshold, we changed the gain set with 0.1 intervals so that a threshold set was limited to 0.4–3 times RMS value. The threshold was automatically set at the value which simultaneously minimizes high-frequency band (40–100 Hz) synchronization for all the EEG channels and maximizes μ (8–12 Hz) desynchronization in the EEG channel of interest (e.g., C3/C4 for hand movements). Specifically, the threshold was decided automatically by finding the minimal value from the summated high-frequency band synchronization and μ band desynchronization vector (an example of finding the proper gain was shown in Figure 2).

2.2. Validation With Simulated EEG/EMG Data

To mathematically verify ERASE, simulated EMG and EEG were generated, subjected to ERASE, and its performance was assessed by several metrics as follows. Here, simulated EMG was also used as reference EMG artifacts for the experimental data processing.

2.2.1. Generating Simulated EMG

Simulated EMG was generated as the reference EMG artifacts, using the approach below (Li et al., 2018):

1. The Hodgkin-Huxley model was used to simulate extracellular current/action potentials (AP). For skeletal muscles, the Hodgkin-Huxley model is a widely accepted model for simulating extracellular current (Hodgkin and Huxley, 1952).

2. Single fiber action potentials (SFAP) were generated with a volume conduction model, defined as follows (Duchene and Hogrel, 2000):

\begin{array}{l} V_{E} (z, y) = K {\int_{S_{1}} \frac{\partial e (z)}{\partial z} \cdot \frac{1}{r} d S + \int_{S} d S \int_{- \infty}^{+ \infty} \frac{\partial^{2} e (z)}{\partial z^{2}} \cdot \frac{1}{r} d z \\ - \int_{S_{2}} \frac{\partial e (z)}{\partial z} \cdot \frac{1}{r} d S} & (7) \end{array}

where V_E is the SFAP, e(z) is the extracellular current (from step 1 above), z and y are the axial and radial directions, respectively, S₁ and S₂ are the fiber sections at the fiber ends, and r is the distance between the surface element, and dS is the observation point. The equation above was discretized to generate the SFAP using known parameter values from the literature, including fiber length, endplate position, observation position, etc. (Stegeman and Linssen, 1992; Muthukumaraswamy, 2013).

3. A Gaussian distribution with 0 mean and standard deviation (SD) = 2.5 mm (Stegeman and Linssen, 1992) was used to depict the endplate positions. The voltage propagation velocities were considered as a Gaussian distribution with an average of 4 m/s and SD = 0.125 (Stegeman and Linssen, 1992). A total of 100 SFAPs were first generated and their average served as one activation of the motor unit action potential (MUAP).

4. A Poisson process was employed to model the firing rate of the MUAPs (as defined in Stegeman and Linssen, 1992). The EMG firing rate and amplitude were assumed to increase during the movements. Hence, different firing rates were applied to each state (idle vs. movement). Specifically, two constant firing rates were used for idle (40 spikes/second) and movement (100 sipkes/second), respectively, in both experimental and simulation verification. For each muscle, the new Poisson process with the same initial firing rate (20 spikes/second) was launched to generate the time points of the firing of MUAP.

5. Eight different facial muscles, including bilateral frontalis, temporalis, masseter, and trapezius were simulated for each session (one session denoted one record, which included several trials). Each muscle's simulated EMG was filtered based on its frequency characteristics (spectra) as described in the literature (Muthukumaraswamy, 2013).

This approach was also summarized in Figure 1A. The above approach ensured that the simulated reference EMG was dependent on the contaminant EMG artifacts to some degree. Therefore, the effectiveness in removing EMG artifacts by using the simulated EMG should be similar with that obtained by using the real EMG.

FIGURE 1

Figure 1. Diagram of simulation validation for ERASE. (A) Diagram for generating simulated EMG. (B) Diagram for generating simulated EEG. (C) Diagram for effectiveness test. (D) Diagram for false positive test. (E) Diagram for sensitivity test. (F) Diagram for data processing of real EEG. AP, action potential; SFAP, single fiber action potential; MUAP, motor unit action potential.

These simulated EMG artifacts were then combined with the simulated EEG data (as described below) as separate channels (the “EMG channels”). These separate channels were located at different positions on the edge of brain topographic map (Supplementary Figure 8). The coordinates of these positions corresponded to the approximate muscle locations on the head.

2.2.2. Generating Simulated EEG

Simulated EEG data was generated (in Figure 1B) by using the approach in (Grouiller et al., 2007):

1. Simulated EEG was created using a linear mixture of five Gaussian noises. Each Gaussian noise was bandpass filtered in different frequency bands (1–30, 20–40, 40–80, 80–100, and 100–200 Hz). The amplitude and variance of each Gaussian noise were adjusted to fit the average values of real EEG data.

2. To mimic the spatial correlation between EEG channels due to volume conduction, a smoothing convolution was performed across channels to increase the spatial correlation amongst adjacent channels. The smoothing convolution kernel was a Gaussian function with a standard deviation equal to four channels. We considered the last channels on the list were correlated with the first channels on the list. Note that this resulted in consecutive channels being highly correlated (i.e., channel 1 becomes highly correlated to channel 2, channel 2 to 3, etc).

Simulated EEG data sets were generated with a temporal structure similar to real EEG. More specifically, each simulated dataset contained 32 channels, the maximal amplitude was set to 60 μV with a variance of 30 μV², and the sampling rate was set to 2,000 Hz. A total of 5 min of simulated EEG was generated. Theoretically, the number of EEG channels (still meet the minimum for ICA run) has no effect on our approach, so we can select any number of EEG channels as long as we feel it is appropriate.

2.2.3. Performance Assessment

The performance of ERASE was assessed in two simulated scenarios. In Scenario 1, ERASE was tested across an increasing number of contaminated EEG channels. Here, three types of simulated EMG (frontalis, temporalis, and masseter) were generated (as in section 2.2.1). In the first iteration, 200 sets of simulated EEG (32 channels) were generated (6 channels were contaminated). Specifically, pairs of simulated EEG channels were then contaminated by a single simulated EMG type, and this was repeated until all 3 types were exhausted. Simulated EMG was multiplied by a randomly generated weight factor and added to the simulated EEG. For each simulated EMG, 2 normally distributed pseudorandom numbers were generated (mean: 0 and standard deviation: 1) and normalized to act as weight factors. Subsequently, the simulated EMG without multiplication was combined with the contaminated EEG set as separate channels. This combined dataset was subjected to ICA. This process was repeated with 12, 18, 24, and 30 simulated EEG channels, in which groups of 4, 6, 8, and 10 simulated EEG channels, respectively, were contaminated with a single EMG type.

For Scenario 2, ERASE was tested on simulated EEG across an increasing burden of EMG contamination. Here, 1 set of simulated EMG was generated for each of the following types: frontalis, temporalis, masseter, trapezius, and eye blinks (as in section 2.2.1). Also, 200 sets of simulated EEG were generated as described above. For each set of simulated EEG, 6 randomly chosen simulated EEG channels were contaminated by a single EMG type. More specifically, the simulated EMG was multiplied by a randomly generated weight factor and then added to the simulated EEG. For each simulated EMG, 6 normally distributed pseudorandom numbers were generated (mean: 0 and standard deviation: 1) and normalized to act as weight factors. Likewise, the simulated EMG without multiplication were combined with the contaminated set of EEG as separate channels, thereby acting as reference EMG artifacts (e.g., there were 32 EEG channels, the simulated EMG was the 33rd, 34th,…37th channel, depending on how many simulated EMG channels were used). This combined dataset was subjected to ICA. This process was repeated, each time incrementally adding another simulated EMG type to an additional six randomly chosen simulated EEG channels, until all EMG types were exhausted. Note that to simplify the process, one simulated EEG channel was contaminated with no more than one type of simulated EMG in both scenarios.

Each combined simulated EEG/EMG set was then subjected to the ICA algorithm (FastICA). The performance was assessed by calculating the effectiveness, false positive rate and sensitivity (described below) across all of the above data sets. Note that simulated EEG (section 2.2.2), does not include any physiological brain features and therefore is easier to be distinguished from simulated EMG as compared to real EEG and EMG. Therefore, a threshold was not necessary to define artifact ICs. Instead, for each reference EMG channel, the IC with the highest coefficient was defined as an artifact IC.

2.2.3.1. Effectiveness

To compare how well EMG artifacts were removed between ERASE and the conventional ICA, the effectiveness of both methods were compared. Here, effectiveness was defined as the ratio of the amount of simulated EEG signal in the artifact ICs. Effectiveness was expressed as artifact index (AI), which was defined as the following:

\begin{array}{l} A I = \frac{(| a_{i + 1, l}^{*} | + \dots + | a_{i + j, l}^{*} |) / j}{(| a_{1, l} | + \dots + | a_{i, l} |) / (32 - j)} \\ (6 \leq j \leq 30, i + j = 32, 1 \leq l \leq 37) & (8) \end{array}

where l is the lth mixing matrix column corresponding to artifact ICs (refer to Equation 9), a is the coefficients corresponding to uncontaminated channels (the “non-artifacts coefficients”), a^* is the coefficients corresponding to contaminated channels (“artifacts coefficients”), and i and j are the number of uncontaminated channels, and contaminated channels, respectively. Equation (8) denotes the ratio of the average mixing matrix coefficient values in the contaminated rows and those in the uncontaminated rows in the artifact IC columns (referred to Equation 9). A larger artifact index indicated that this identified artifact IC contained more artifacts, but less simulated EEG signal.

We applied ERASE to two conditions of EMG-artifacts contaminated simulated EEG: with and without separate channels (Figure 1C). Generally, we referred to this latter case as “conventional ICA” condition. The artifact indices calculated for the conditions with and without simulated EMG were statistically compared in order to verify the effectiveness of our model. Note that only criteria 2 from the rejection criteria above was employed in the conventional ICA condition.

2.2.3.2. False Positive Rate

To determine how frequently ERASE would erroneously detect sources that were independent of the reference EMG artifacts, a false positive rate was designed. More specifically, false positive rate determined how often contaminant artifacts (N_{c_E}, Equation 1) that were independent of the reference EMG artifacts ( $n_{c_{M}}^{*}$ , equation 1) were erroneously “pushed” into the artifact ICs. To assess this, N_{c_E}, composed of Gaussian random noise (mean 0, S.D. 30), was used to further contaminate the simulated EEG. Simulated EMG was still used as reference EMG artifacts, $n_{c_{M}}^{*}$ , and combined with the contaminated EEG data sets and acted as separate channels as described above. To simplify the process, only one type of (independent) Gaussian random artifacts was employed for both scenarios. These steps are summarized in Figure 1D.

Here, we define the mixing matrix column corresponding to the artifact IC (the “artifact IC column”) as follows:

\begin{array}{l} \begin{matrix} V = {(a_{1, l}, \dots, a_{i, l}, a_{i + 1, l}^{*}, \dots, a_{i + j, l}^{*}, ã_{i + j + 1, l}, \dots, ã_{i + j + k, l})}^{T} \\ (6 \leq j \leq 30, i + j = 32, 1 \leq l \leq 37, 1 \leq k \leq 5,) \end{matrix} & (9) \end{array}

where a, a^*, i, j, and l are defined in Equation (8), V is the “artifact IC column,” ã denotes the coefficients corresponding to the EMG channels, k are the number of EMG channels. Given that the ICs were normalized to unit variance, the coefficients contained in a given column of the mixing matrix (Equation 5) can be interpreted as relative loads by which this IC contributed to the mixed signals. Therefore, in the artifact IC columns, large coefficient values are usually associated with the channels that were contaminated by artifacts, whereas uncontaminated channels would have relatively smaller coefficients. Satisfying this inequality (10) meant that the corresponding artifacts were detected in these artifact ICs. After running ERASE on the combined simulated EEG data sets, all the artifact ICs, which were decided by the position of the maximal absolute values in the corresponding mixing matrix rows representing the EMG channels, were found. A false positive was formally defined as:

\begin{array}{l} \frac{(| a_{i, l}^{*} | + \dots + | a_{i + j, l}^{*} |)}{j} - \frac{(| a_{1, l} | + \dots + | a_{i, l} |)}{32 - j} \\ > max (| ã_{i + j + 1, l} |, \dots, | ã_{i + j + k, l} |) * 5 % & (10) \end{array}

where max(|ã_{i+j+1, l}|, …, |ã_{i+j+k, l}|) is defined as the maximal artifacts coefficient. If this inequality was satisfied in any artifact ICs, i.e., contaminated channels coefficients exceeded those for the uncontaminated channels by 5% of the maximal artifacts coefficient (denoted as the threshold in Figure 2), a false positive event occurred. The ratio of those false positive events in 200 simulated data sets was the false positive rate.

FIGURE 2

Figure 2. Effectiveness, false positive, and sensitivity of ERASE on the simulated data. (A,B) Artifact indices, i.e., effectiveness in Scenario 1 and 2, respectively. (C,D) False positive for Scenario 1 and 2, respectively. (E,F) Sensitivity results for Scenario 1 and 2, respectively. Scenario 1 was considered as the fact that the number of contaminated EEG was changed for the performance assessment of ERASE. Here, the number of contaminated EEG channels were chosen from 2 to 10 for each simulated EMG artifacts (three EMG artifacts were used here, so numbers are 6–30 in figures). We increased the number of added EMG channels in Scenario 2 for the performance assessment of ERASE. The number of added EMG channel was chosen from 1 to 5. Asterisks indicate a significant difference between two data sets (Wilcoxon Rank-Sum Test), and the significance level = ***p<0.001.

2.2.3.3. Sensitivity

A sensitivity was designed to assess if contaminant EMG artifacts N_{c_E} were accurately “pushed” into the artifact ICs, when contaminant EMG artifacts, N_{c_E}, were dependent on reference EMG artifacts $n_{c_{M}}^{*}$ . Sensitivity was defined as the ratio of events in which the contaminant EMG artifacts N_{c_E} were detected in the artifact ICs after running ERASE. The simulated EMG which served as both contaminant EMG artifacts and reference EMG artifacts here, were used to contaminate the simulated EEG data as described above and also acted as separate channels. These combined simulated EEG data sets were subjected to ICA. These steps are summarized in Figure 1E. Referring to the artifact ICs columns, the sensitivity in our study was calculated by the ratio of events in which the inequality (10) was satisfied in all the artifact ICs columns in corresponding 200 simulated data sets.

2.3. Validation With Real EEG

The ability of ERASE to automatically reject real EMG artifacts from real EEG was assessed as follows.

2.3.1. Experiments

This study was approved by the Institutional Review Boards of the University of California, Irvine. Healthy subjects with no history of neurological conditions were recruited for this study. Subjects were fitted with a 64-channel EEG cap (ActiCap, Brain Products, Gilching, Germany) and asked to perform repetitive fist clenching and unclenching of the dominant hand while their EEG signals were acquired by two, linked NeXus-32 systems (Mind Media, Herten, Netherlands). EMG was recorded from the bilateral frontalis, left temporalis to masseter, right temporalis to masseter and bilateral trapezius using a MP150 system (BIOPAC, Goleta CA), respectively. The subjects were asked to sit in front of a computer screen, which prompted them to alternate between idling (for 5 s) or hand fist-clenching (for 2 s). This was repeated for a total of 10 trials over a 100 s-long session. At least 2 sessions were performed by each subject. The EEG and EMG data were recorded at 2,048 and 4,000 Hz sampling rates, respectively.

2.3.2. Experimental Data Processing

The whole data processing for real EEG was summarized in Figure 1F. For the experimental EEG, both real and simulated EMG artifacts acted as separate channels and were not mixed to any EEG channels. The combined EEG/EMG data were bandpass filtered from 3 to 100 Hz (3rd order, forward-backward filter with no phase distortion). Note that the 100 Hz upper cutoff was chosen since the amplifiers attached to MP150 have a 100 Hz low-pass filter in hardware. Each trial, comprised 1-s idle time followed by 2-s movement, was identified and extracted from the combined EEG. Due to the non-stationarity of EEG, ICA decomposition was just applied to concatenated EEG trial datasets for each session (each run). The FastICA version in the EEGLAB toolbox (Delorme and Makeig, 2004) was used to run ICA on the EEG trials data for all the subjects in the three conditions (simulated EMG, real EMG, and conventional ICA conditions). The artifact ICs were rejected as above. Short-time Fourier transform was applied to EEG trials and the signal power in different frequency bands [μ band: 8–12 Hz, high frequency (HF) band: 40–100 Hz] was compared across all ICA conditions. All the data were z-scored after the time-frequency decomposition, which were separately normalized to the statistics of the EEG during the idling epochs.

The z-scored power of the μ and high-frequency bands during idle and movement was statistically compared for all the EEG channels using a Wilcoxon rank-sum test. The z-scored power of μ/high-frequency band for each channel was then topographically mapped. For channels where there was no significant difference between idle and movement, these values were nulled (set to 0).

Given that the absolute power at any frequency bands was reduced after artifact ICs were rejected from the original signal, we used the z-scored power of high-frequency band to calculate the decrease of high frequency after removal of EMG artifacts for each subject. The percent reduction was defined as below:

\begin{array}{l} P D = \frac{| \sum_{c = 1}^{C} \sum_{n = 1}^{N} P_{z}^{b} (X_{n}^{c}) - \sum_{c = 1}^{C} \sum_{n = 1}^{N} P_{z}^{a} (X_{n}^{c}) |}{\sum_{c = 1}^{C} \sum_{n = 1}^{N} P_{z}^{b} (X_{n}^{c})} \times 100 % & (11) \end{array}

where $P_{z}^{b}$ is the z-scored power of high frequency in baseline, $P_{z}^{a}$ is the z-scored power of high frequency after removal of EMG artifacts, X are the EEG trials data, n is the nth trial, N is the total number of available trials for each subject, c is the cth channel, C is the total number of EEG channels. The rationale for this approach is that the high-frequency signal is dominated by EMG in EEG and any reduction in the high-frequency band power was considered as the reduction in EMG. The Wilcoxon rank-sum test was also employed to examine the difference of the μ band power in the C3/C4 channel and the high-frequency band power in all of the channels during movement between all combinations of the four conditions (baseline, after ERASE with simulated EMG, after ERASE with real EMG, and after conventional ICA).

3. Results

3.1. Simulation Verification

a. Effectiveness: The artifact indices for ERASE and conventional ICA were summarized in Figures 2A,B. Across all parameters, the artifact indices calculated by equation 8 in ERASE were significantly larger than those in the conventional ICA model (Wilcoxon rank-sum test, P < 0.001). This indicates that ERASE has better effectiveness in removing the EMG artifacts compared to the conventional ICA model.

b. False Positive Rate: The artifacts and non-artifacts coefficients for the false positive were summarized in Figures 2C,D. Across all parameters, there was a 0% rate of false positive, as defined by Inequality 10 (details regarding thresholds are in Supplementary Table 1). Furthermore, the values of artifacts and non-artifacts coefficients were small and were not significantly different from one another. This indicates that signals independent of the reference artifacts were not erroneously “pushed” into artifact ICs by ERASE.

c. Sensitivity: The artifacts and non-artifacts coefficients for the sensitivity were summarized in Figures 2E,F. Across all parameters, there was a 100% rate of sensitivity, as defined by Equation (10) (details about threshold were shown in Supplementary Table 2). Furthermore, the values of artifacts and non-artifacts coefficients were significantly different from one another (Wilcoxon rank-sum test, P < 0.001) and the values of non-artifacts coefficients were always significantly lower. This indicates that contaminant EMG artifacts can be accurately identified by ERASE.

3.2. Experimental Verification

A total of eight subjects gave their informed consent to participate in the study. All subjects were healthy male volunteers between 20 and 50 years old. A total of 12 sessions (120 trials) were performed across all subjects.

An example of selecting the proper threshold for ERASE was shown in Figure 3. Table 1 summarized the effect of ERASE on the μ-band and high-frequency band across all subjects. Overall, the high-frequency band power, typically dominated by EMG artifacts, was reduced by 75.31% using ERASE with real EMG, by 63.46% through ERASE with simulated EMG and only by 48.88% with the conventional ICA approach (Table 1). At the group level (Figure 4A), the z-scored power of the high-frequency band during movement was significantly reduced after running ERASE. Furthermore, the high-frequency band power in the real EMG and simulated EMG conditions were both significantly lower than that in the conventional ICA (Figure 4A). However, there was no difference between running ERASE with either the real EMG or simulated EMG. At the group level, there was no change in the μ band (only for C3/C4 channel) among the four conditions based on all trials (Figure 4B). This indicates that the expected μ-band desynchronization phenomenon during the fist-clenching task was not adversely affected by ICA.

FIGURE 3

Figure 3. Representative z-scored μ and high-frequency band power from Subject 1 after running ICA (as in section 2.3.2), using simulated EMG and real EMG. (A) Average of z-scored μ power in the C3 channel in the real EMG condition. (B) Average of z-scored high-frequency band power for all the EEG channels in the real EMG condition. (C) Values added by z-scored high-frequency band and μ power in the real EMG condition. (D) Average of z-scored μ power in the C3 channel in the simulated EMG condition. (E) Average of z-scored high-frequency band power for all the EEG channels in the simulated EMG condition. (F) Values added by z-scored high-frequency band and μ power in the simulated EMG condition. After running ICA on the combined EEG/EMG data, the artifact ICs were rejected based on the criteria depicted above for each threshold. After the joint consideration, the final threshold was 1.5 for real EMG condition, and 0.4 for simulated EMG, which were denoted by arrows.

TABLE 1

Table 1. Change of HF band power in different conditions for each subject [mean ± standard deviation (S.D)].

FIGURE 4

Figure 4. Comparison of z-scored power of μ/high-frequency band during movement in different conditions. (A) z-scored power of high-frequency band during the movement for all the EEG channels under different conditions. (B) z-scored power of μ band during movement in the C3/C4 channel under different conditions. Data were from all the subjects with a total of 120 trials. ***Significant differences between the two datasets (p < 0.001).

A representative example of significant changes in z-scored power of μ and high-frequency bands during idle and movement is shown in Figure 5. The z-scored power of high-frequency band during movement was reduced after running ERASE, and the power in both ERASE conditions was smaller than that in the conventional ICA (Figures 5E–H, Supplementary Figures 1E–H–7E–H). A representative example of time series also showed the same findings (Supplementary Figure 9). In Figures 5A–D, the μ desynchronization during right hand movement was well-preserved after running ICA in all the conditions and was localized to the C3 channel and surrounding electrodes. This localized μ desynchronization was also present around the C3 channel for Subjects 3–7 (Supplementary Figures 2A–D–6A–D) and around the C4 channels for Subjects 2 and 8 (Supplementary Figures 1A–D–7A–D). The channels exhibiting μ desynchronization were always contralateral to the hand movement except Subjects 2 and 8 (Supplementary Figures 1, 7). Combined with the findings above, this indicates that ERASE did not disturb the spatial distribution of the expected brain features underlying the motor task of interest.

FIGURE 5

Figure 5. Brain topography maps for Subject 1 displaying the z-scored power of μ band and high-frequency band in different conditions [baseline, after ERASE with real EMG and with simulated EMG and after conventional ICA, (A–D) for μ band and (E–H) for high-frequency band] on the Subject 1. Channels whose z-scored power of μ/high-frequency band were not significantly different between idle and movement states (Wilcoxon rank-sum test) were nulled (values were set to zero). The significance threshold was P < 0.01 for μ band and P < 0.05 for high-frequency band.

4. Discussion

Here, we proposed a modified ICA model that combined reference EMG artifacts with EEG data to facilitate an enhanced automated removal of EMG artifacts. We tested and validated this method using both simulated and actual EEG during hand movement. We found that it had high sensitivity at detecting EMG artifacts and an extremely low false positive rate (Figures 2C–F, Supplementary Tables 1, 2). With simulated data, ERASE effectively removed a large proportion of the EMG artifacts (Figures 2A,B). It also removed EMG artifacts in real EEG recordings, while preserving the expected μ desynchronization associated with movement (Figure 5 and Supplementary Figures 1–7). This may also indicate that our approach can remove any potential confounding overlap between EMG and EEG and thereby improves the confidence that low-frequency brain features extracted from ERASE are mostly EMG-free (which cannot be achieved by a simple low pass filter).

These results suggest that using reference EMG artifacts can force ICA to “learn” and detect the EMG artifacts by forcing the contaminant EMG within EEG into a minimal number of ICs. We established an operational definition (rejection criteria) for identifying EMG artifacts components to enable automated component rejection and thereby minimizing user bias. Compared to conventional ICA, results showed that ERASE removed on average 26% more EMG artifacts from EEG data than conventional ICA (Figure 4 and Table 1), which indicated that ERASE improved the ICA algorithm. Although this approach may require slightly more preparation time to record real EMG, it is still possible to use it in situations where real EMG recordings were either not possible or not available by substituting it with simulated EMG (Figure 4). The process to generate simulated EMG mimics a typical situation where a subject will likely generate increased EMG activity during movements. Although this process cannot perfectly emulate the time-varied EMG signals, it simulates the firing rate, amplitude and spectrum of each muscle to ensure the statistical parameters of EMG signals in different states to some extent. Those statistical parameters are exactly required for successfully running ERASE/ICA, and as such, may explain why the effectiveness at removing EMG artifacts by using simulated EMG was found to be similar to that obtained by using the real EMG. The advantages and novelty of ERASE are discussed in further detail below.

First, ERASE directly introduces reference EMG artifacts into the ICA model as prior knowledge to more accurately maximize the separation between EMG and EEG ICs as well as to minimize the computational complexity of removing EMG with respect to existing cICA approaches. On the other hand, previously reported forms of reference ICA, such as spatially cICA (Hesse and James, 2006; Akhtar et al., 2012), involves first performing an initial run of ICA on a previous EEG data segment, followed by manually selecting ICs believed to represent EMG sources as EMG reference. Subsequently, ICA is run iteratively on current EEG data to find all the ICs which have a strong correlation with pre-defined EMG reference. By comparison, ERASE real EMG (where available) or simulated EMG as the reference. Since such references are expected to provide a more reliable representation of the ground truth for EMG sources, ERASE is likely more reliable and systematic than other forms of ICA. In addition, ERASE utilizes properties of mixing matrix for identification of artifact ICs to avoid the complicated computation of optimization problem, which is employed in temporally cICA (James and Gibson, 2003; James and Hesse, 2004; Lu and Rajapakse, 2005). When combined with the use of a simulated or simultaneously recorded real EMG used in ERASE, it is not necessary to perform iterative run for the ICA process. Our results (Figures 4, 5 and Table 1) showed that ERASE can remove most of EMG artifacts and preserve expected brain features just after a single ICA run.

Unlike previously reported versions of ICA, ERASE does not require manual intervention. More specifically, ICA is typically run directly on EEG, and the operator manually rejects the artifact ICs. The effectiveness of this step depends substantially on the experience of researchers and may be biased due to subjectivity. Hence, ERASE can minimize the biases of researchers and improve the efficiency of artifacts rejection. As mentioned in the introduction, automated rejection is not necessarily unique to ERASE (Delorme et al., 2001, 2007; Nicolaou and Nasuto, 2007; Nolan et al., 2010; Mognon et al., 2011; Daly et al., 2012, 2013; Wu et al., 2018; Vaidya et al., 2019), given that other methods, such as cICA can also involve automatic IC rejection when prior knowledge of EMG signals is available (Hesse and James, 2006; Akhtar et al., 2012; Urigüen and Garcia-Zapirain, 2015). Previously reported EMG artifacts removal methods also proposed automated rejection techniques, in which some classifiers were built to classify the ICs into EMG sources and EEG sources based on ICs statistical features (Nolan et al., 2010; Gabsteiger et al., 2014; Wu et al., 2018). However, one unique aspect of ERASE compared to these prior reports is that rejection criteria are based on physiological features of both EEG and EMG for automated EMG artifacts rejection procedure (section 2.1.2), which makes ERASE more focused on preserving relevant EEG phenomenon. In addition, although we used μ band as the physiological feature of interest due to its known modulation during human motor behaviors, this automated rejection technique can be adjusted to other bands of interest in other EEG studies (e.g., β band from pre-frontal areas in cognitive studies).

ERASE was validated with both simulated and behavioral EEG data, whereas the physiological information and properties of EEG are completely overlooked in other EMG artifact removal studies using ICA, cICA, or other BSS methods. Researchers typically employ simulated or synthetic EEG data to validate corresponding artifacts removal algorithms. Based on various metrics (discussed in the paragraph below), all of these algorithms have declared themselves effective at removing EMG artifacts, but have not answered a critical question as to whether the information encoded in EEG is retained after artifact removal. Here, we show that EEG μ band modulation that typically underlies hand movement is preserved or even enhanced after ERASE. Combined with the reduction of EMG elsewhere in the EEG, these findings unequivocally demonstrate that meaningful EEG is retained. Such a demonstration has been generally absent from the validation of other ICA approaches. Note that the high-frequency band in EEG is assumed to be EMG artifact since the skull filters out all neurogenic high-frequency information (Whitham et al., 2007; Dannhauer et al., 2011; Lanfer et al., 2012). Therefore, any attenuation of the high-frequency band in EEG can be considered as removal of EMG artifacts. Also, μ desynchronization in EEG and ECoG signals is a well-known modulation underlying motor behavior (Pfurtscheller and Da Silva, 1999; Miller et al., 2007; Schalk et al., 2007; Jiang et al., 2020). Therefore, these two bands are used jointly as an indicator that our approach is effective at removing EMG artifacts while preserving the brain physiological features underlying human movement.

ICA was selected as the basis of ERASE due to the fact that ICA has typically been shown to have superior performance to most other artifact removal methods. For example, CCA, a popular algorithm for artifact removal, does not outperform ICA at removing EMG artifacts from EEG (McMenamin et al., 2010; Escudero et al., 2011; Safieddine et al., 2012; Urigüen and Garcia-Zapirain, 2015) as well as at removing ECG and EOG artifacts (Romero et al., 2008, 2009; Pham et al., 2011; Delorme et al., 2012; Evans et al., 2012; Daly et al., 2013). Some BSS methods and source decomposition methods have been combined for removal of EMG artifacts (e.g., EEMD-CCA Chen et al., 2015, 2017; Mucarquer et al., 2019, as well as EEMD-ICA Mijovic et al., 2010). Both EMD and EEMD are single-channel techniques, so EEMD-CCA and EEMD-ICA are only tested on the fact that a few channels of EEG recording are involved. Since running ICA or CCA sometimes is time-consuming, EEMD-CCA or EEMD-ICA probably is typically less than ideal. While the ideas of using EMG data to make the EMG rejection method more systematic were presented before (i.e., cICA mentioned above), to the best of our knowledge, there were no prior attempts to fundamentally improve the separation of EMG artifacts from EEG by adding extra reference artifacts to the ICA process as in current study. The closest to this was our initial preliminary report on this technique (Li et al., 2018), followed by Richer et al. (2019). In Richer et al. (2019), extra EMG channels were also used and combined with EEG to “force” ICA algorithm to detect more EMG artifact. However, the EMG artifact ICs were arbitrarily rejected manually, there was no rigorous mathematical proof or experimental verification for the approach and there was no assessment to ensure that meaningful EEG data was preserved during this process.

There are no uniformly accepted performance metrics for artifacts removal algorithm in practical EEG. In most studies (Crespo-Garcia et al., 2008; McMENAMIN et al., 2009; Vos et al., 2010), the performance of EMG artifact removal is typically assessed, in part, by visual inspection. Although highly subjective, it may give an indication with respect to whether the algorithm has improved the quality of the EEG signal or has distorted one or more time intervals or frequency bands. This is also performed in our study and the results are showed in Supplementary Figure 9. For simulated and synthetic EEG, there are many metrics employed for the assessment of the performance of artifact removal approaches, since the “ground truth” in such scenarios are explicitly known (De Clercq et al., 2006; Mijovic et al., 2010; Akhtar et al., 2012; Albera et al., 2012; Safieddine et al., 2012; Zhang et al., 2012). The most widely used metrics for performance assessment are listed in the literature (Urigüen and Garcia-Zapirain, 2015) [e.g., the relative root mean squared error (RRMSE), the signal to artifact ratio (SAR), etc.]. In our work, we used standard statistical metrics (false positive, sensitivity, effectiveness) to validate our approach in simulated EEG. Meanwhile, time-series validation results are shown in Supplementary Figures 10–12 in the Appendix. It should also be emphasized that since researchers in each study used different data sets, a head-to-head performance comparison across various artifact removal approaches in EEG studies is difficult. Therefore, we further validated our novel approach by using and assessing both the elimination of EMG and the preservation of the EEG features underlying motor behaviors. We propose that such an evaluation is employed for validation of effectiveness of any new artifact rejection approaches in future studies, as it is otherwise not possible to know whether the methodology aggressively removed EMG as well as erroneously eliminating the signals of interest if real brain signals underlying behavioral experiments are not collected. In order to further evaluate ERASE, we employed ERASE on EEG recorded from traumatic brain injury (TBI) patients with hemicraniectomy, and ERASE can show a better effectiveness compared to conventional ICA (Li et al., 2020).

The main limitations of ERASE are that it is still impossible to remove all EMG from EEG. The main reason originates from the assumptions of our theoretical model. Namely, our model requires that the reference EMG artifacts and contaminant EMG artifacts are completely dependent. However, contaminant EMG artifacts in EEG data cannot have complete dependence on the real EMG sources due to several reasons. For example, the signal propagation path between the true source and recording electrodes may significantly distort observed EMG. It is also not possible to include all the reference EMG artifacts which could be contributing to EEG, because some head and neck muscles are not easily recorded at the surface. In addition, another assumption in our theoretical derivation is that the reference EMG artifacts are independent of EEG. Also, the reference EMG electrodes are close to potential EEG sources (such as EMG reference electrodes located over the frontalis and temporalis muscles), so it is difficult to ensure that the reference EMG artifacts channels in this study are fully independent of EEG. Hence, the identified artifact ICs may still contain some contribution from EEG. However, EEG from these areas is not highly involved in motor tasks in our studies. Finally, this was only tested using EEG from 8 healthy subjects, and future more rigorous experimental validation can be performed on an even larger cohort, with EEG underlying other types of behaviors (e.g., cognitive tasks), or in those with neurological conditions.

5. Conclusion

Here, we proposed a modified ICA model that can automatically remove EMG artifacts by combining reference EMG artifacts with EEG. This new approach can more effectively remove EMG artifacts from EEG while preserving the expected brain features underlying motor behavior. Also, the approach proposed in this work is automated, which minimizes experimenter bias and speeds up analysis. The utilization of the simulated EMG as the reference EMG source potentially extends the application of this approach. The EEG recovered by our approach can provide more confidence for further neuroscience analysis. Meanwhile, future work will focus on testing ERASE on EEG from patient populations, and adapting it for real-time applications, such as in the BCI system.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The studies involving human participants were reviewed and approved by Institutional Review Boards of the University of California, Irvine. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

YL, AD, and MS: study design and manuscript preparation. PW: data collection and manuscript review. RF, MV, and CL: data collection. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the National Institutes of Health, R01NS094748.

Conflict of Interest

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2020.597941/full#supplementary-material

We showed the results of z-scored power of μ band (8–12 Hz) and high frequency band (40–100 Hz) in different conditions (baseline, after ERASE with real EMG and simulated EMG and after conventional ICA) for Subject 2–8 (Supplementary Figures 1–7), the 2D image of the electrode locations (Supplementary Figure 8), time series band-pass filtered with 40–100 Hz for one trial in different conditions (Supplementary Figure 9), time series validation results by using simulated EEG with simulated EMG (Supplementary Figure 10) and real EMG (Supplementary Figure 11), and performance comparison between ERASE and conventional ICA by a common metric. Also, Supplementary Tables 1, 2 list the average thresholds for false positive and sensitivity test, respectively. Supplementary Table 3 lists the average artifact indices for the simulated EMG and conventional ICA conditions.

References

Akhtar, M. T., Mitsuhashi, W., and James, C. J. (2012). Employing spatially constrained ICA and wavelet denoising, for automatic removal of artifacts from multichannel EEG data. Signal Process. 92, 401–416. doi: 10.1016/j.sigpro.2011.08.005