Refinement of High-Gamma EEG Features From TBI Patients With Hemicraniectomy Using an ICA Informed by Simulated Myoelectric Artifacts.

Recent studies have shown the ability to record high-γ signals (80-160 Hz) in electroencephalogram (EEG) from traumatic brain injury (TBI) patients who have had hemicraniectomies. However, extraction of the movement-related high-γ remains challenging due to a confounding bandwidth overlap with surface electromyogram (EMG) artifacts related to facial and head movements. In our previous work, we described an augmented independent component analysis (ICA) approach for removal of EMG artifacts from EEG, and referred to as EMG Reduction by Adding Sources of EMG (ERASE). Here, we tested this algorithm on EEG recorded from six TBI patients with hemicraniectomies while they performed a thumb flexion task. ERASE removed a mean of 52 ± 12% (mean ± S.E.M) (maximum 73%) of EMG artifacts. In contrast, conventional ICA removed a mean of 27 ± 19% (mean ± S.E.M) of EMG artifacts from EEG. In particular, high-γ synchronization was significantly improved in the contralateral hand motor cortex area within the hemicraniectomy site after ERASE was applied. A more sophisticated measure of high-γ complexity is the fractal dimension (FD). Here, we computed the FD of EEG high-γ on each channel. Relative FD of high-γ was defined as that the FD in move state was subtracted by FD in idle state. We found relative FD of high-γ over hemicraniectomy after applying ERASE were strongly correlated to the amplitude of finger flexion force. Results showed that significant correlation coefficients across the electrodes related to thumb flexion averaged ~0.76, while the coefficients across the homologous electrodes in non-hemicraniectomy areas were nearly 0. After conventional ICA, a correlation between relative FD of high-γ and force remained high in both hemicraniectomy areas (up to 0.86) and non-hemicraniectomy areas (up to 0.81). Across all subjects, an average of 83% of electrodes significantly correlated with force was located in the hemicraniectomy areas after applying ERASE. After conventional ICA, only 19% of electrodes with significant correlations were located in the hemicraniectomy. These results indicated that the new approach isolated electrophysiological features during finger motor activation while selectively removing confounding EMG artifacts. This approach removed EMG artifacts that can contaminate high-gamma activity recorded over the hemicraniectomy.

1. Introduction

Conventional EEG has a poor signal-to-noise ratio in the high-γ band due to the low-pass filter characteristics of the skull and scalp (Nunez and Srinivasan, 2006; Luck, 2014). Meanwhile, movement and force are strongly encoded in the high-γ band activity from brain signals (Crone et al., 1998; Mehring et al., 2003; Pfurtscheller et al., 2003; Miller et al., 2007a; Schalk et al., 2007; Flint et al., 2012a,b, 2014, 2016). Traumatic brain injury (TBI) patients with hemicraniectomy may be a useful model for human electrophysiology with high bandwidth and spatiotemporal resolution (Voytek et al., 2010; Vaidya et al., 2019). In particular, substantial high-γ band power can be detected in these patients' electroencephalogram (EEG) due to the absence of the skull in the hemicraniectomy area (referred to as hEEG) (Dannhauer et al., 2011; Lanfer et al., 2012). However, the extraction of the high-γ band features from hEEG in TBI patients remains challenging due to the large confounding spectral overlap between the high-γ band of EEG and EMG, which is primarily caused by facial and head movement and has broad bandwidth (Duchene and Hogrel, 2000; Goncharova et al., 2003; Fatourechi et al., 2007; Dalal et al., 2011).

To surmount this challenge, we tested the ability of our novel approach, referred to as EMG Reduction by Adding Sources of EMG (ERASE) (detailed in our previous work: Li et al., 2018, 2020), to minimize the EMG artifacts in hEEG. ERASE is a modified ICA model that can automatically remove EMG artifacts by combining reference EMG artifacts with EEG. In this new approach, real or simulated EMG from neck and head muscles were defined as reference EMG artifacts, and combined with EEG as extra channels. ERASE was validated using both simulated and experimentally-recorded EEG and EMG from healthy subjects, and shown to remove EMG artifacts while preserving relevant electrophysiological features underlying motor behaviors (Li et al., 2018, 2020). Simulation results showed ERASE removed EMG artifacts from EEG significantly more effectively than conventional ICA and had a low false positive rate and high sensitivity (Li et al., 2020). By comparison, there have been several EMG artifact removal approaches previously reported, but a critical difference with ERASE was that these approaches did not experimentally test their ability to remove EMG while preserving relevant electrophysiological features underlying human behavior. Such methods include ones based on independent component analysis (ICA) (James and Hesse, 2004; Nolan et al., 2010), constrained ICA (cICA) (Lu and Rajapakse, 2005; Romero et al., 2008), canonical correlation analysis (CCA) (Safieddine et al., 2012; Mowla et al., 2015), empirical mode decomposition (EMD) (Mourad and Niazy, 2013; Zeng et al., 2013), ensemble empirical mode decomposition (EEMD) (Wu and Huang, 2009; Chen et al., 2014), EEMD-CCA (Chen et al., 2017; Mucarquer et al., 2019), as well as EEMD-ICA (Mijovic et al., 2010). In this study, we applied the ERASE technique to EEG from 6 TBI patients while they performed isometric finger flexion, and found the residual high-γ signal in hEEG was correlated to force level.

2. Methods

2.1. Experiments

This study was approved by the Institutional Review Boards of the University of California, Irvine, Northwestern University and the Rancho Los Amigos National Rehabilitation Center (RLANRC). TBI patients with hemicraniectomy and mild to moderate weakness on the contralesional hand were recruited for our study. Subjects were fitted with a 128-electrode EEG cap (ActiCap, Brain Products, Gilching, Germany) and asked to perform a thumb flexion task on the contralesional side while the EEG signals were acquired at 2,000 Hz (Neuroport, Blackrock Microsystems, Salt Lake City, UT). The thumb flexion force was simultaneously measured by a load cell sensor. The subjects were instructed to flex their thumbs on the affected side to apply varying levels of force on the load cell sensor to move a computer cursor to acquire targets in a 1D, random-force target task. Each flexion event was set to be 2 s long (target displayed for 2 s) with a 3–5 s interval between each event. Subjects nominally performed 20 thumb flexions in each 120-s run at RLANRC and 53 thumb flexions in each 300 s-long run at Northwestern University. All of the trails were used in the analysis described below.

2.2. Experimental Data Processing

The EEG from TBI patients with hemicraniectomy were subjected to following processing approach. Due to TBI patients'limited ability to tolerate testing, real EMG was not collected for use in the ERASE algorithm. Instead, only simulated EMG was used as reference EMG artifacts in this study. The approach to generate simulated EMG was detailed in our prior work (Li et al., 2018, 2020), and we summarized the 5-step approach as below:

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

Single fiber action potentials (SFAP) were generated with a volume conduction model described in Stegeman and Linssen (1992), Duchene and Hogrel (2000), and Li et al. (2018).

A total of 100 SFAPs were first generated and their average served as one activation of the motor unit action potential (MUAP).

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 hand/finger movements. Hence, firing rates which were proportional to thumb flexion force in different states (idle vs. movement) were applied to different states. For each muscle, the new Poisson process with the same initial firing rate (20 spikes/s) was launched to generate the time points of the firing of MUAP. Firing rate will change with thumb flexion force but be limited to a maximum of 5 times of initial firing rate (i.e., 100 spikes/s).

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 ensured that the simulated reference EMG was dependent on the contaminant EMG artifacts to some degree.

EEG from TBI patients with hemicraniectomy and simulated EMG data generated by the approach above were combined (simulated EMG acted as separate electrodes and were not mixed with any EEG signal), and the combined EEG/EMG data were subjected to a 3–200 Hz 3rd-order bandpass filter. All trials were identified and extracted from the EEG/EMG data. Each trial was defined as 1-s idle time (remaining 2–4 s of idle time discarded) followed by 2-s movement. The extracted trials were concatenated (which were referred to as baseline condition). Since EEG always includes other unexplained noise, leading to long-term EEG non-stationary, ICA decomposition was applied to concatenated EEG trials [FastICA algorithm, EEGLAB toolbox (Delorme and Makeig, 2004)] which contained the entire broad bandwidth. After running ICA on this combined EEG/EMG data, an automated artifact independent components (ICs) rejection process was applied (Li et al., 2018, 2020) as follows. The ICs containing EMG artifacts (the “artifact ICs”) were identified and rejected using an automated procedure which optimized two rejection criteria based on mixing matrix. These two rejection criteria were described as below: first, ICs whose maximal absolute value of coefficients corresponds to a hat band electrode were identified as artifact ICs and rejected. Second, ICs whose absolute value of coefficients in the corresponding EMG rows were above a defined threshold, were rejected. Here, threshold was calculated from the root mean square (RMS) values of coefficients in the mixing matrix rows corresponding to the EMG channels. The rationale for why adding reference EMG artifacts was effective at removing EMG artifacts was detailed in Li et al. (2020). Note that for comparison purposes, conventional ICA was also applied to concatenated EEG trials without extra simulated EMG channels, and artifact ICs were rejected manually by the experimenter. Short-time Fourier transform was applied to EEG trials under three conditions (baseline, ERASE with simulated EMG, conventional ICA) and data were z-scored to the power in idle time in each trial after the time-frequency decomposition. Average z-scored signal power across all the trials in different frequency bands (μ band: 8–12 Hz, high-γ band: 80–160 Hz) were compared across all conditions (baseline, ERASE with simulated EMG, conventional ICA). To assess the effectiveness of EMG artifact removal, we first assumed that signals within the high-γ band from non-hemicraniectomy were EMG artifacts due to the skull essentially filtering out all neurogenic high-γ signals. Therefore, we calculated the percent reduction (PR) for high-γ band across all the electrodes in non-hemicraniectomy areas by using the equation below:

where is the z-scored power of high-γ band in baseline, is the z-scored power of high-γ band after removal of EMG artifacts, X are the EEG trials data, i is the ith trial, N is the total number of available trials for each subject, c is the cth channel in non-hemicraniectomy areas, C is the total number of EEG channels in non-hemicraniectomy areas. For each subject, percent reduction was calculated after running ERASE and conventional ICA, respectively.

2.3. Signal-to-Noise Ratio (SNR) Calculation

In previous studies, high-γ band of brain activities was demonstrated to synchronize to movement (Flint et al., 2014, 2016; Wang et al., 2017; McCrimmon et al., 2018; Branco et al., 2019) while μ band was desynchronized (Miller et al., 2007a; Schalk et al., 2007; Jiang et al., 2020). To evaluate if these electrophysiological features and information content underlying movement were retained after ERASE, the signal-to-noise ratio (SNR) between movement and idle was calculated for electrodes expected to be related to thumb movements, including C3, C5, C1, FCC5h, FCC3h, CCP5h, and CCP3h for subjects with left-sided hemicraniectomy, or C4, C2, C6, FCC6h, FCC4h, CCP4h, and CCP6h for those with right-sided hemicraniectomy. These electrodes will be subsequently referred to as hand motor electrodes (see Supplementary Figure 11). Note that the homologous electrodes on the non-hemicraniectomy side will be referred to as the contralesional electrodes. Generally, most of hand motor electrodes were located in the hemicraniectomy areas. Here, we treated the high-γ power during movement as signal and the baseline high-γ power during idle as noise. The SNR for each trial before and after applying ICA was calculated as in (McCrimmon et al., 2018). The segmented trials were concatenated and the resulting time series was referred to as X(t). For each electrode, the EEG was subjected to bandpass filtered from 80 to 160 Hz (4th order Butterworth filters) to extract the high-γ band, referred to as Xγ(t). The resulting signal was then squared to obtain the instantaneous high-γ band power, . was low-pass filtered (4 Hz, 4th order Butterworth filter) to create the envelope signal, Pγ. The average power of high-γ band during movement and idling time were calculated, denoted as and , respectively. The SNR was defined as . This processing was also applied to μ band. We calculated the SNR at both of μ and high-γ bands for all the available trials of each subject. A Wilcoxon rank-sum test with Holm Bonferroni correction was used to compare the SNR between all pair combinations of conditions (baseline, ERASE with simulated EMG, conventional ICA) for both bands.

2.4. EEG High-γ Verification

To measure the success of ERASE, we assessed the extent to which high-γ movement-related information content was retained after it has been applied. Previous studies have demonstrated that high-γ band power increased during muscle activation or movement (Crone et al., 1998; Miller et al., 2007b; Zhuang et al., 2010; Flint et al., 2014; Wang et al., 2017; McCrimmon et al., 2018) and that fractal dimension (FD) can be used as a measure of brain activity by quantifying the complexity of EEG (Acharya et al., 2005; Subha et al., 2010). Also, the FD of EEG signal has been shown to change with the level of force (Liu et al., 2005). Since relative FD (change in FD between idle and move states) can be an indicator of cortical activity modulation during movement (Liu et al., 2005), the correlation between thumb flexion force and the high-γ relative FD will be used to assess if the recovered high-γ in hEEG still retains encoding of movement. The correlation between the high-γ relative FD and the thumb flexion force was calculated as follows. First, the FD of the EEG high-γ was calculated as follows for the move and idle epochs in each trial, respectively (Katz, 1988; Esteller et al., 2001; Liu et al., 2005):

where N is the total number of time points to be analyzed for each epoch (2,000 for idle epoch, 4,000 for movement epoch), L is the sum of the Euclidean distances between successive data vectors, and d is the Euclidean distance between the first data vector and the vector that provides the farthest distance. Data vector was composed of successive time points, which were . Any adjacent data vectors included no overlapped time points. Since the FD value was dependent on quantization units in this algorithm (Katz, 1988), 1 ms was chosen as the quantization unit of time (i.e., one data vector was composed of two time points, hence, multiple data vectors were included in each epoch), and 1μV as that of the EEG potential. Here, we calculated the Euclidean distances for all the data vectors in the linear space.Next, the relative FD, defined as that FD value during idle epoch subtracted from that during move epoch, was calculated for each trial.

For each subject, the mean force during the move epoch for each trial was calculated. Then the maximum and minimum of these mean force were found. Since the resolution of high-γ in hEEG was not sufficient to precisely decode continuous force, the mean force was evenly discretized into 10 force levels from minimum to maximum. Subsequently, the relative FD values were averaged over the trials at each force level for each subject. The correlation coefficient (R, |R| denoted the absolute value of correlation coefficient) between force level and relative FD values for all electrodes was calculated with Pearson Correlation. The significance of correlation was calculated by considering the correlation coefficients as a t distribution. The t-value can be calculated as below (Soper et al., 1917):

where R is correlation coefficient, N is the sample number. In our work, the correlation coefficients with no significance (P > 0.05) were set as zero, and the ones with significance denoted as significant R-value (or significant |R|-value for absolute value).

3. Results

A total of 58 sessions (including 1,771 trials) were collected from 6 TBI patients with hemicraniectomy (subject demographics summarized in Table 1).

Table 1

Subject demographics.

Subject	Age	Sex	HA side	Contralesional hand
S1	23	Female	Left	Right hand
S2	34	Male	Left	Right hand
S3	30	Male	Right	Left hand
S4	40	Male	Right	Left hand
S5	29	Male	Right	Left hand
S6	56	Female	Left	Right hand

HA, hemicraniectomy areas.

3.1. Brain Features Evaluation

After processing of EEG as described in section 2.2, representative examples showed that the z-scored EEG high-γ power in the non-hemicraniectomy areas (NHAs) was reduced while μ desynchronization (negative SNR value) and high-γ synchronization (increased SNR) in the motor-related areas [in the hemicraniectomy areas (HAs)] were apparent after ERASE (Figure 1 and Supplementary Figures 1–5). In baseline, the average Z-scored high-γ power was not different between HA and non-HA across all subjects. Only after ERASE did the average Z-scored high-γ power during movement across all the electrodes in the HAs became significantly larger than that across all the electrodes in the NHAs (P-value < 0.05, Wilcoxon rank-sum test, Table 2). The z-scored high-γ power during movement in the NHA was reduced by an average of 52.03 ± 12.08% (mean ± S.E.M) across all subjects (Table 2), indicating that artifacts were removed. On the other hand, the corresponding reduction in non-HA for conventional ICA condition averaged 26.50 ± 18.91% (mean ± S.E.M) across all subjects (Table 2).

Figure 1

Brain topography maps of z-scored μ and high-γ power before and after artifacts rejection with ERASE and conventional ICA on the Subject 6. Only electrodes whose z-scored power of μ and high-γ during idle time and movement were significantly different were shown (Wilcoxon rank-sum test). P-value for significant difference was 0.05. The dots outlined the position of the electrodes, and details of electrode position can refer to Supplementary Figure 11. Colors denoted the z-scored power of μ or high-γ in corresponding electrodes. The purple outline in each subfigure denoted the HA. A-C: μ band desynchronization under different conditions (baseline, ERASE with simulated EMG, conventional ICA). D-F: high-γ band synchronization under different conditions.

Table 2

z-scored high-γ power in different conditions for each subject.

Subject	S1	S2	S3	S4	S5	S6
Number of trials	261	97	864	230	209	110
Baseline
Mean z-scored in NHA	0.26	0.13	0.21	0.24	0.26	0.15
Mean z-scored in HA	0.25	0.12	0.21	0.23	0.35	0.13
P-value	0.35	0.23	0.47	0.42	0.67	0.62
ERASE with simulated EMG
Mean z-scored in NHA	0.16	0.06	0.12	0.11	0.12	0.04
Mean z-scored in HA	0.24	0.10	0.19	0.19	0.34	0.12
P-value	0.02	0.03	0.04	0.04	0.03	0.03
Conventional ICA
Mean z-scored in NHA	0.10	0.11	0.16	0.18	0.24	0.11
Mean z-scored in HA	0.14	0.12	0.16	0.15	0.2	0.12
P-value	0.12	0.63	0.57	0.17	0.96	0.61

Mean z-score denoted the average z-scored high-γ power during movement over all the available trials. The mean z-scored high-γ power in the HAs got averaged across all the electrodes in the HAs. Meantime, The mean z-scored high-γ power in the NHAs got averaged across all the electrodes in the NHAs.

HA, Hemicraniectomy area; NHA, Non-HA.

P-value: comparison of z-scored high-γ power between HAs and NHAs for each subject (Wilcoxon rank-sum test).

The high-γ synchronization described above was seen in the C3 (for subjects with left hemicraniectomy, right thumb flexion)/C4 (subjects with right hemicraniectomy, left thumb flexion) electrode in 4 out of 6 subjects in baseline. When considering all hand motor electrodes, high-γ synchronization was observed in 5 out of 6 subjects in baseline (Figures 2A,B). The high-γ synchronization at the C3/C4 electrode was significantly larger after ERASE in 5 out of 6 subjects compared to baseline condition. However, the C3/C4 high-γ synchronization after conventional ICA was never significantly larger than those in baseline (Figure 2A). The high-γ synchronization from hand motor electrodes was significantly larger after ERASE compared to those both in baseline and after conventional ICA for all the subjects. However, except in Subject 1, the high-γ synchronization from hand motor electrodes in the conventional ICA condition was smaller than those in baseline (Figure 2B). These findings suggested that high-γ remaining in the HAs after ERASE was movement-related EEG, and ERASE was more effective at preserving the EEG features.

Figure 2

Mean SNR [± standard deviation (S.D.)] across all available trials for each subject in different conditions (baseline, after ERASE and after conventional ICA). S1 was the abbreviation for Subject 1, and so on. The asterisks indicated the significant differences between the two datasets, and the significance level = ***p < 0.001, level = **p < 0.01, and level = *p < 0.05 (Wilcoxon rank-sum test). (A) SNR of high-γ band. Data were from the C3/C4 electrode. (B) SNR of high-γ band. Data were from hand motor electrodes. (C) SNR of μ band. Data were from the C3/C4 electrode. (D) SNR of μ band. Data were from hand motor electrodes.

In order to further verify these high-γ activities were neurogenic and related to movement, the μ desynchronization was calculated as described in section 2.3. The μ desynchronization was present in 5 out of 6 subjects in baseline at C3/C4. After applying ERASE, it was present in all subjects at C3/C4, and the magnitude of desynchronization increased significantly in 4 subjects and was not significantly different in the remaining 2 (Figure 2C). In contrast, after conventional ICA, the μ desynchronization was no longer present in 1 out of 6 subjects. Furthermore, the magnitude of desynchronization was significantly reduced in 4 subjects (Figure 2C). Across all of the hand motor electrodes, μ desynchronization was present in all subjects in baseline. After applying ERASE, the magnitude of desynchronization was significantly larger in 5 out of 6 subjects compared to baseline condition (Figure 2D). After conventional ICA, the magnitude of desynchronization was significantly reduced in all subjects. The magnitude of desynchronization was always significantly larger after applying ERASE compared to that of the conventional ICA (Figure 2D).

3.2. EEG High-γ Verification

A strong correlation between thumb flexion force and the relative FD of high-γ was primarily found in the HAs for all the subjects after ERASE (Figure 3). After conventional ICA, either there were no strong correlations in the HAs (Figure 3 and Supplementary Figures 1, 4), or a strong correlation was seen in both HAs and non-HAs (Figure 3 and Supplementary Figures 2, 3, 5, 6).

Figure 3

Correlation between the relative FD of high-γ and the amplitude of thumb flexion force in different conditions (baseline, after ERASE and after conventional ICA). Only the significant correlation coefficients were showed. P-value for significant difference was 0.05. The dots outlined the position of the electrodes, and details of electrode position can refer to Supplementary Figure 11. Colors denoted the z-scored power of μ or high-γ in corresponding electrodes. The purple outline in each subfigure denoted the HA. (S1–S6) denoted Subject 1–6, respectively.

The correlations between relative FD and thumb flexion force for each condition are summarized in Table 3. In baseline, 27.34% of electrodes with significant correlation were located within HAs across all subjects. Subsequently, 83.33 and 18.99% of electrodes with significant correlation were located within HAs after applying ERASE and the conventional ICA, respectively (Table 3, Figure 4, and Supplementary Figures 6–10).

Table 3

Average of significant |R|-values in hand motor electrodes and contralesional electrodes in different conditions.

Subject	S1	S2	S3	S4	S5	S6
Baseline
|R|-value in hand motor	0.61	0.76	0	0.72	0.53	0
SCE number in hand motor	4	1	0	1	3	0
|R|-value in contralesional	0.73	0.81	0.77	0	0.68	0.93
SCE number in contralesional	1	1	1	0	1	1
Total number of SCE	31	8	14	4	8	21
Proportion of SCE in HA (%)	35.48	37.50	14.29	25	37.5	14.29
ERASE with simulated EMG
|R|-value in hand motor	0.62	0.83	0.80	0.74	0.69	0.90
SCE number in hand motor	3	2	2	1	1	1
|R|-value in contralesional	0.68	0	0	0	0	0
SCE number in contralesional	1	0	0	0	0	0
Total number of SCE	12	5	3	2	4	2
Proportion of SCE in HA (%)	83.33	100	66.67	100	50	100
Conventional ICA
|R|-value in hand motor	0	0.76	0.75	0	0.80	0.86
SCE number in hand motor	0	1	2	0	1	1
|R|-value in contralesional	0.73	0	0	0.81	0	0
SCE number in contralesional	1	0	0	1	0	0
Total number of SCE	10	10	8	14	11	13
Proportion of SCE in HA (%)	0	30	50	0	18.19	15.38

SCE was significant correlation electrodes.

Hand motor denoted hand motor electrodes, which included C3, C5, C1, FCC5h, FCC3h, CCP5h, and CCP3h for right hand, or C4, C2, C6, FCC6h, FCC4h, CCP4h, and CCP6h for left hand. Contralesional denoted the contralesional electrodes, which were the homologous motor electrodes in the non-hemicraniectomy area side.

Figure 4

Bar graphs for showing the electrodes with significant correlation in three conditions (baseline, after ERASE with simulated EMG and after running conventional ICA). Here, the correlation coefficients were absolute values from 0 to 1. Data were from Subject 1. Blue bars denoted the electrodes with significant correlation in the NHAs, and red bars were the ones in HAs. (A) The electrodes with significant correlation in baseline. (B) The electrodes with significant correlation after ERASE with simulated EMG. (C) The electrodes with significant correlation after conventional ICA.

In baseline, the average significant |R|-value in hand motor electrodes was 0.436 across all subjects (subjects with no significant correlation were treated as though the significant |R|-value was 0). This was lower than the significant |R|-values on the contralesional electrodes area (0.654 from Table 3). The average significant |R|-values within the hand motor electrodes were 0.762 and 0.528 across all subjects for ERASE and conventional ICA conditions, respectively (Table 3). By comparison, ERASE removed the presence of any correlation in the contralesional electrodes areas except in 1 subject (only one electrode was still correlated, average significant |R| across all the subjects was 0.11). After conventional ICA, the average significant |R|-value within contralesional electrodes across all the subjects was 0.26 (Table 3).

4. Discussion

Our novel approach, ERASE, was demonstrated to be an effective tool at removing EMG artifacts from EEG during thumb flexion by TBI patients with HAs while preserving the expected underlying EEG features. Specifically, both μ desynchronization and high-γ synchronization during movement in TBI patients were found in anatomically expected areas, while μ desynchronization was not seen on the contralesional side. This indicated that the high-γ synchronization in hEEG can be more confidently interpreted as cortical activity and not EMG (Figure 1, Supplementary Figures 1–5). Also, the SNR of high-γ and μ were both significantly improved compared to conventional ICA and baseline conditions (Figure 2). The relative FD of high-γ was strongly correlated to the thumb flexion force primarily within the HAs after applying our new approach (Figures 3, 4 and Supplementary Figures 6–10). These results indicated that ERASE can help to isolate cortical sources of high-γ by significantly reducing the presence of EMG artifacts.

Information about the primary kinetics appeared to be preserved in the high-γ signal after ERASE. At baseline, high correlation between high-γ relative FD and thumb flexion force was similarly present in both HA and NHAs. After applying ERASE, hand motor electrodes located in the HAs typically had the strongest positive correlation across all the subjects, indicating that the new algorithm selectively removed EMG artifacts from the NHAs. This implied that there were fewer artifacts in the HA after ERASE. Given that ERASE helps to isolate electrophysiological features during hand motor movements while removing confounding EMG artifacts, the justification for using TBI patients with hemicraniectomy as a model for studying high-bandwidth EEG signals is bolstered.

Our recent study on EEG from TBI patients with hemicraniectomy demonstrated that the high-γ over HAs was capable of effectively decoding the thumb flexion force (Vaidya et al., 2019). In this study, we demonstrated that high-γ over HAs may contain information about thumb flexion force, as demonstrated by a strong correlation between relative FD of high-γ and flexion force. Anecdotally, neither this study nor our prior work showed a strong direct correlation between hEEG high-γ power and thumb flexion force as seen in ECoG signals (Li et al., 2018). Most likely, the presence of scalp still interferes with the high-γ signal enough that such robust features are not preserved when compared to subdurally (Pistohl et al., 2012; Wang et al., 2017) and epidurally (Flint et al., 2016) recorded ECoG signals. Alternatively, the removal of EMG may also be overly aggressive. More specifically, while EMG is effectively removed, a portion of neurogenic high-γ component was also simultaneously removed as well. This highlights a potential tradeoff between EMG removal and neural information preservation in ERASE. In the future, this tradeoff may potentially be minimized by selecting the rejection threshold in a manner that jointly optimizes high-γ encoding as well as EMG reduction. In another study on TBI patients with hemicraniectomy (Voytek et al., 2010), both μ desynchronization and high-γ synchronization underlying movement can be found in hEEG. In our work, we not only found that this modulation was associated with movements (Figures 1, 2, Supplementary Figures 1–5), but also showed that there was a strong correlation between the relative FD of high- and the amplitude of thumb flexion force after employing ERASE (Figures 3, 4, Supplementary Figures 6–10).

In the approach to generate simulated EMG, we assumed that simulated EMG activity is increased during movements and reduced during idling. This simulates a typical situation where a subject will likely generate increased EMG activity during movements. Although the approach to generate the EMG signals cannot emulate the time-varied EMG signals rigorously for each time point, it simulated the firing rate, amplitude and spectrum of each muscle to ensure the statistical parameters of EMG signals in different states to some extent. Those statistics are exactly required for running ERASE/ICA. Therefore, the effectiveness at removing EMG artifacts by using the simulated EMG should be similar with that obtained by using the real EMG. We demonstrated that using simulated EMG as reference artifacts can achieve similar effectiveness at removing EMG artifacts (while preserving the brain features) as real EMG did in Li et al. (2020).

Even after ERASE with simulated EMG, the overall reduction of high-γ power was 52% (Table 2), and there was still one subject with EEG electrode in the NHA with a strong correlation between the high-γ relative FD and thumb flexion force (Table 3). These findings indicated that not all of the EMG artifacts were easily removed. This is likely due to the fact that simulated EMG artifacts cannot precisely mimic real EMG artifacts since many components of EMG are still difficult to simulate (e.g., the neural sources of real EMG artifacts). We hypothesize that using real EMG as the reference artifacts may lead to better EMG artifact rejection. However, since real EMG is not always available or possible to collect, future work will involve developing an algorithm that further improves upon the EMG removal. This may involve improving the learning algorithms, such as adaptive system recognition and system recognition with artificial neural network, to characterize the EMG artifacts. In addition, future work will involve making the technique more computationally efficient, such that it can be used in real time for neurorehabilitation applications, such as in BCI systems.

5. Conclusion

Our new approach, ERASE, as described in our prior work (refer to Li et al., 2020) can also be applied to effectively remove EMG artifacts from hEEG. In particular, we have demonstrated the first approach that ERASE can potentially remove the confounding overlap between EMG and high-γ signals, and preserve the expected brain signal features underlying motor behavior. The retained high-γ activities demonstrated the expected increase during thumb flexion in contralateral hand motor cortex area (within the hemicraniectomy site). Moreover, the relative FD of the EEG high-γ from the HAs after applying the new approach was demonstrated to be strongly correlated to the amplitude of thumb flexion force. Therefore, this approach may allow researchers to confidently use the resulting high-γ signals for subsequent analysis or practical applications.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, 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, Northwestern University, Rancho Los Amigos National Rehabilitation Center (RLANRC). 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. MV, CL, and RF: data collection. All authors contributed to the article and approved the submitted version.

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.