Electromyogram (EMG) Removal by Adding Sources of EMG (ERASE) - A novel ICA-based algorithm for removing myoelectric artifacts from EEG - Part 2

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 2000 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 seconds long (target displayed for 2 seconds) with a 3-5 seconds interval between each event. Subjects nominally performed 20 thumb flexions in each 120-s run at RLANRC and 53 thumb flexions in each 300s-long run at Northwestern University. All of the trails were used in the analysis described below.

Only simulated EMG was used as reference EMG artifacts, as described in Section 2.2.1 in Part 1 in Supplementary Materials. The EEG from TBI patients with hemicraniectomy were subjected to the same processing approach described in Section 2.3.2 in Part 1 in Supplementary Materials. Briefly, the EEG and EMG data 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 3${}^{\text{rd}}$-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. 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 [20]) which contained the entire broad bandwidth. The automated artifacts IC rejection process described in Section 2.1.2 in Part 1 was applied. Short-time Fourier transform was applied to EEG trials and the signal power in different frequency bands ($\mu$ band: 8 to 12 Hz, high-$\gamma$ band: 80 to 160 Hz) were compared across all conditions (baseline, ERASE with simulated EMG, conventional ICA). All the data were z-scored after the time-frequency decomposition. Also, we calculated the percent reduction for high-$\gamma$ band across all electrodes in non-hemicraniectomy areas by using equation (11) in Part 1.

To evaluate if the 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 Fig. 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-$\gamma$ power during movement as signal and the baseline high-$\gamma$ power during idle as noise. The SNR for each trial before and after applying ICA was calculated as in [21]. The segmented trials were concatenated as described in Section 2.3.2 in Part 1 in Supplementary and the resulting time series was referred to as $X(t)$. For each electrode, the EEG was subjected to bandpass filtered from 80-160 Hz (4th order Butterworth filters) to extract the high-$\gamma$ band, referred to as $X_{\gamma}(t)$. The resulting signal was then squared to obtain the instantaneous high-$\gamma$ band power, $X_{\gamma}^{2}(t)$. $X_{\gamma}^{2}(t)$ was low-pass filtered (4 Hz, 4th order Butterworth filter) to create the envelope signal, $P_{\gamma}$. The average power of high-$\gamma$ band during movement and idling time were calculated, denoted as $\overline{P_{m,\gamma}}$ and $\overline{P_{i,\gamma}}$, respectively. The SNR was defined as $10\times\log_{10}{(\overline{P_{m,\gamma}}/\overline{P_{i,\gamma}})}$. This processing was also applied to $\mu$ band. We calculated the SNR at both of $\mu$ and high-$\gamma$ bands for all the available trials of each subject. A Wilcoxon rank-sum test was used to compare the SNR between all pair combinations of conditions (baseline, ERASE with simulated EMG, conventional ICA) for both bands.

To measure the success of ERASE, we assessed the extent to which high-$\gamma$ movement-related information content was retained after it has been applied. Previous studies have demonstrated that high-$\gamma$ band power increased during muscle activation or movement [4, 22, 21, 23, 24, 25] and that fractal dimension (FD) can be used as a measure of brain activity by quantifying the complexity of EEG [26, 27]. Also, the FD of EEG signal has been shown to change with the level of force [28]. Since relative FD (change in FD between idle and move states) can be an indicator of cortical activity modulation during movement [28], the correlation between thumb flexion force and the high-$\gamma$ relative FD will be used to assess if the recovered high-$\gamma$ in hEEG still retains encoding of movement. The correlation between the high-$\gamma$ relative FD and the thumb flexion force was calculated as follows. First, the FD of the EEG high-$\gamma$ was calculated as follows for the move and idle epochs in each trial, respectively [28, 29, 30]:

where N is the total number of time points to be analyzed for each epoch (2000 for idle epoch, 4000 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 less than $\frac{N}{2}$. Any adjacent data vectors included no overlapped time points. Since the FD value was dependent on quantization units in this algorithm [29], 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$\mu$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-$\gamma$ 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 [31]:

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).

After processing of EEG as described in Section 2.3.2 in part 1, representative examples showed that the z-scored EEG high-$\gamma$ power in the non-hemicraniectomy areas (NHAs) was reduced while $\mu$ desynchronization (negative SNR value) and high-$\gamma$ synchronization (increased SNR) in the motor-related areas (in the hemicraniectomy areas (HAs)) were apparent after ERASE (Fig. 1 and Supplementary Figs. 1 to 5). In baseline, the average Z-scored high-$\gamma$ power was not different between HA and non-HA across all subjects. Only after ERASE did the average Z-scored high-$\gamma$ 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-$\gamma$ power during movement in the NHA was reduced by an average of 52.03 $\pm$ 12.08% (mean $\pm$ 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 $\pm$ 18.91% (mean $\pm$ S.E.M) across all subjects (Table 2).

The high-$\gamma$ 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-$\gamma$ synchronization was observed in 5 out of 6 subjects in baseline (Fig. 2 C and D). The high-$\gamma$ 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-$\gamma$ synchronization after conventional ICA was never significantly larger than those in baseline (Fig. 2C). The high-$\gamma$ 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-$\gamma$ synchronization from hand motor electrodes in the conventional ICA condition was smaller than those in baseline (Fig. 2D). These findings suggested that high-$\gamma$ remaining in the HAs after ERASE was movement-related EEG, and ERASE was more effective at preserving the EEG features.

In order to further verify these high-$\gamma$ activities were neurogenic and related to movement, the $\mu$ desynchronization was calculated as described in Section 2.3. The $\mu$ 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 (Fig. 2 A). In contrast, after conventional ICA, the $\mu$ desynchronization was no longer present in 1 out of 6 subjects. Furthermore, the magnitude of desynchronization was significantly reduced in 4 subjects (Fig. 2 A). Across all of the hand motor electrodes, $\mu$ 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 (Fig. 2 B). 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 (Fig. 2 B).

A strong correlation between thumb flexion force and the relative FD of high-$\gamma$ was primarily found in the HAs for all the subjects after ERASE (Fig. 3). After conventional ICA, either there were no strong correlations in the HAs (Fig. 3 S1 and S4), or a strong correlation was seen in both HAs and non-HAs (Fig. 3 S2, S3, S5, and S6).

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, Fig. 4 and Supplementary Figs. 6-10).

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).