Scalp electroencephalograms over ipsilateral sensorimotor cortex reflect contraction patterns of unilateral finger muscles

A variety of neural substrates are implicated in the initiation, coordination, and stabilization of voluntary movements underpinned by adaptive contraction and relaxation of agonist and antagonist muscles. To achieve such flexible and purposeful control of the human body, brain systems exhibit extensive modulation during the transition from resting state to motor execution and to maintain proper joint impedance. However, the neural structures contributing to such sensorimotor control under unconstrained and naturalistic conditions are not fully characterized. To elucidate which brain regions are implicated in generating and coordinating voluntary movements, we employed a physiologically inspired, two-stage method to decode relaxation and three patterns of contraction in unilateral finger muscles (i.e., extension, flexion, and co-contraction) from high-density scalp electroencephalograms (EEG). The decoder consisted of two parts employed in series. The first discriminated between relaxation and contraction. If the EEG data were discriminated as contraction, the second stage then discriminated among the three contraction patterns. Despite the difficulty in dissociating detailed contraction patterns of muscles within a limb from scalp EEG signals, the decoder performance was higher than chance-level by 2-fold in the four-class classification. Moreover, weighted features in the trained decoders revealed EEG features differentially contributing to decoding performance. During the first stage, consistent with previous reports, weighted features were localized around sensorimotor cortex (SM1) contralateral to the activated fingers, while those during the second stage were localized around ipsilateral SM1. The loci of these weighted features suggested that the coordination of unilateral finger muscles induced different signaling patterns in ipsilateral SM1 contributing to motor control. Weighted EEG features enabled a deeper understanding of human sensorimotor processing as well as of a more naturalistic control of brain-computer interfaces.


Introduction
During the execution of physical movements, agonist muscles contract and antagonist muscles relax to generate joint torque, while both simultaneously contract (co-contract) to achieve joint stability and limb movement accuracy to maintain posture or hold objects ( Burdet et al., 2001 ;Gomi and Osu, 1998 ;Gribble et al., 2003 ). Studies exploiting invasive recordings or functional imaging techniques revealed that flexible and dynamic coordination of skeletal muscles requires the involvement of various neural substrates, including the cerebral cortex, basal ganglia, thalamus, cerebellum, and spinal cord ( Bostan and Strick, 2018 ;Gallivan et al., 2011 ;Leo et al., 2016 ;Weiler et al., 2019 ). For instance, at the level of the spinal cord, inhibitory circuits for antagonist muscles contingent on agonist activation (i.e., reciprocal inhibition) are pathways ( Edelman et al., 2019 ;Hochberg et al., 2012 ;Kato et al., 2019 ). As decoding the presence of movement attempts under naturalistic conditions is the crucial aspect of a non-invasive BCI, scalp electroencephalogram (EEG) signals have been utilized as the input of such decoders ( Pfurtscheller et al., 2006 ;Soekadar et al., 2016 ). In EEG-based BCI, the sensorimotor rhythm (SMR), an EEG component at 8-30 Hz around the SM1 contralateral to the targeted hand, is a prevalent feature because it is attenuated during execution of voluntary movement or kinesthetic motor imagery ( Neuper and Pfurtscheller, 2001 ;Pfurtscheller, 1992 ). Such power attenuation putatively originates from desynchronization of the thalamo-cortical loop (i.e., event-related desynchronization; SMR-ERD; Hindriks and van Putten, 2013 ;Pfurtscheller and Lopes Da Silva, 1999 ). EEG-triggered transcranial magnetic stimulation (TMS) studies confirmed that the SMR-ERD is accompanied by the disinhibition of intracortical inhibition of SM1 ( Naros et al., 2016 ;Takemi et al., 2013 ). Furthermore, as the locus of SMR-ERD reflects somatotopic organization during voluntary movement ( Pfurtscheller et al., 2006 ), it has been utilized as a reliable biomarker to probe attempted movement in BCIs for clinical treatment of post-stroke hemiplegia ( Biasiucci et al., 2018 ;Mukaino et al., 2014 ;Ramos-Murguialday et al., 2013 ).
Numerous studies have focused on methods to enhance decoder performance for efficient BCIs and practical applications ( Lotte et al., 2018 ;F. 2007 ;Wierzga ł a et al., 2018 ). A promising approach is to leverage the physiological dynamics of brain activity to realize effective decoder performance ( Kamitani and Tong, 2005 ), yet such decoder designs are rarely proposed in the context of BCIs. Here we employed a physiologically inspired, two-stage method to decode relaxation and three patterns of contraction in unilateral finger muscles (i.e., extension, flexion, and co-contraction) using high-density EEG. The decoder consisted of two parts employed in series. The first discriminated between relaxation and contraction. If EEG data at this stage discriminated a contraction, the second stage discriminated among the three contraction patterns. We hypothesized that such a design would efficiently decode movement attempts because cortical dynamics are dynamically modulated during the transition from relaxation to motor execution compared with the differences among contraction patterns ( Miller et al., 2007( Miller et al., , 2009. Moreover, to identify EEG features contributing to the decoding in each stage and to explain the physiological correlates of decoder performance, we measured changes in decoder performance induced by randomizing the EEG features.

Participants
Nineteen neurologically healthy participants (all right-handed, aged 22.1 ± 0.9, 18 men, 1 woman) were explained the aim of the experiments and gave informed written consent. The experiment was approved by the local ethics committee of the Keio University (IRB approval number: 30-78) and performed in accordance with the Declaration of Helsinki. Anonymized data is available on reasonable request, as data protection is required in the ethics agreement.

Experimental setup
Participants were seated on a comfortable chair with armrests in an unshielded experimental room. A display was placed about one meter in front of the chair to show visual feedback of electromyogram (EMG) signals and task instructions (See also 2.4 Experimental Design). EMG signals were recorded from the extensor digitorum communis (EDC) and the flexor digitorum superficialis (FDS) muscles of the right forearm. These two muscles were selected for visual feedback, as they are respectively protagonist muscles of finger extension and flexion. In addition to these finger muscles, EMG signals from the biceps brachii (BB) and triceps brachii (TB) muscles were also recorded to confirm the relax-ation of proximal arm muscle contraction. Additionally, we confirmed the absence of muscle contraction by measuring EMG of the left EDC and FDS muscles during experiments in three randomly selected participants. All surface EMG signals were acquired using Ag/AgCl electrodes ( = 10 mm). The Electrical Geodesics system with 128-channel HydroCel Geodesic Sensor Net (GES 400, Physio16; Electrical Geodesics, Inc.) was used to acquire scalp EEG and EMG signals. EEG electrode impedance levels were kept below 50 k Ω throughout experiments and a 100 Hz online low pass filter was applied to recorded data at a sampling rate of 1000 Hz. Channels placed at CPz and Cz in the extended 10-20 system served as ground and reference, respectively.

Experimental design
The participants conducted four types of tasks: relaxation (Relax), finger extension (Ext), finger flexion (Flx), and co-contraction (Cocont). For the Relax task, participants were asked to keep their eyes open without any overt muscle contractions. For the three contraction tasks (Contraction), participants were asked to isometrically contract the extensor muscles (Ext), the flexor muscles (Flx), or both of these two muscles (Cocont). Additionally, participants were asked to maintain their left hand relaxed through the experiment and keep it on armrest.
At the beginning of the experiment, maximum voluntary contractions (MVC) of the EDC and FDS muscles were measured at once for each motor task to determine the target force level, which was then set at approximately 30% MVC. The targets and their force levels were displayed during motor execution periods. Online-processed EMG signals were used for visual feedback (See 2.4. EMG and EEG data processing). As shown in Fig. 1 A, visual feedback of integrated EMG (iEMG) signals was employed to control force levels generated during experiments. Participants rehearsed the four tasks long enough to reach the displayed targets within 1 s with the proximal muscles relaxed. It took approximately 10 min to get used to the tasks. After rehearsal, participants successfully executed directed tasks as shown in Fig. 1 B.
In the main experiment, participants executed one task, pseudorandomly selected from Relax and the three contraction tasks, in each trial ( Fig. 1 C). One trial started at the "Rest" period, and participants were instructed to keep their eyes open. Subsequently, one of the four tasks was directed at the "Ready" period. At the "Task" period, participants then executed the task for 5 s under the feedback of EMG signals acquired from EDC and FDS muscles. Participants executed 40 trials for each task (i.e., each participant conducted 160 trials).

EMG and EEG data processing
All online and offline analyses of EMG and EEG signals were conducted using MATLAB 2016b (The MathWorks, Inc, Natick, MA, USA). The schematic flow of all offline analyses is shown in Fig. 2 A. Surface EMG signals were filtered with a second-order Butterworth bandpass filter (5-450 Hz) and a 50-Hz notch filter. Filtered EMG signals were rectified and smoothed with 300-ms sliding windows and 33% overlap. The iEMG signals were used for online feedback and offline preprocessing.
Scalp EEG signals were processed in the following steps for offline analyses. EEG electrodes around the eyes and ears were excluded from further analyses to minimize the effect of artifacts related to face movements; therefore, EEG signals recorded from 78 channels were used following analysis ( Debnath et al., 2019 ;Pichiorri et al., 2015 ). A thirdorder Butterworth bandpass filter (3-70 Hz) and 50-Hz notch filter was applied. The common average reference was used to enhance the signalto-noise ratio and to extract localized activities ( McFarland et al., 1997 ). Subsequently, a short-time Fourier transform (STFT) was applied with 1-s sliding windows and 90% overlap. To check the spectro-temporal patterns and their spatial patterns, event-related spectral perturbations (ERSP) were computed using the "Rest" period as the reference and the averaged value of ERSP was used to conduct non-parametric permutation test to compare relaxation and three contraction tasks ( Maris and Oostenveld, 2007 ).
We also conducted cortical source estimation of the ERSPs using sLoreta ( Pascual-Marqui, 2002 ). We used band-power of EEG signals between 1 s and 5 s during each trial to calculate the cross spectra for an inverse solution of sLoreta. A non-parametric permutation test was conducted to visualize the difference between the resting and task periods ( Nichols and Holmes, 2002 ). The frequency of interest was set to the frequency bands of the SMR (i.e., 8-30 Hz) because both components of the SMR in the alpha (8-13 Hz) and beta (14-30 Hz) bands have distinct functional roles. Furthermore, not only alpha and beta, but also the upper and lower frequency components within each band are functionally dissociable ( Hacker et al., 2017 ;Ritter et al., 2009 ;Tsuchimoto et al., 2017 ). Thus, the frequency bands of the SMR were subdivided as follows: low-alpha; 8-10 Hz, high-alpha; 11-13 Hz, low-beta; 14-18 Hz, middle-beta; 19-24 Hz, high-beta; 25-30 Hz ( Babiloni et al., 2011( Babiloni et al., , 2008Iturrate et al., 2018 ;Rangaswamy et al., 2002 ). As we conducted individual statistical test for these 5 bands, we corrected the alpha-level with Bonferroni correction and set as p < 0.01.

Construction of decoders using random forest algorithm
The magnitudes of calculated Fourier coefficients were averaged to subdivide the SMR frequency into the five frequency bands just described. The magnitude of each frequency band was calculated for each window and used as the datasets for the construction of decoders after a third-order Butterworth lowpass filter (0.2 Hz) was applied for each time course of frequency band power. Time points between 1 s to 4 s from the beginning of "Task " period were used as dataset with each time window as individual observation (i.e., no temporal information was considered in the decoder construction). The 0.2 Hz lowpass filter was applied to reduce the flicker of band-power within a trial and improve the signal-to-noise ratio of voluntarily controllable SMR (   The architecture of the two-stage decoder. Note that each decoder was independently constructed and the two then combined after the constructions. the contraction patterns were the signal strengths within the 5 bands acquired from 78 electrodes (i.e., 390 features). All EEG features were used to construct the Decoder-All decoder because the spatial patterns of EEG signals acquired from the whole-head possibly contain information about finger movements. To test this assumption, we also constructed a control decoder (Decoder-Ctrl), which used only features around the contralateral SM1, where motor output is conveyed to the hand muscles ( Bradnam et al., 2012 ;Soteropoulos et al., 2011 ) and is therefore a common region in decoding the presence of right finger movement. Another decoder (Decoder-Ctr) was also constructed using EEG features derived only from the C3 electrode according to the international 10/20 convention.
We constructed decoders using a consecutive two-stage method to dissociate the decoding task into two parts ( Fig. 2 B): one that decoded relaxation/contraction and the second that decoded the three contraction patterns because EEG signals are dynamically modulated during the transition from relaxation to motor execution compared with the difference between contraction patterns. Using this physiologically inspired approach, Hayashi et al. (2019) reconstructed produced force from scalp EEG signals by dividing their decoders into two components (i.e., one to decode the presence/absence of contraction and one to decode the force production). Decoder-1 (Dec-1) discriminated between relaxation and contraction at the first stage of decoding, and Decoder-2 (Dec-2) then discriminated the type of contraction pattern for EEG datasets decoded as contraction by Dec-1. Dec-2 was trained with the EEG datasets acquired only during the contraction task. This two-stage procedure was also essential to separately identify features that contributed to the decoding of relaxation or contraction and those to the decoding of muscle contraction patterns (See also 2.7. Calculation of feature-importance scores).
Both Dec-1 and Dec-2 were constructed using the Random Forests classifier (RF) ( Breiman, 2001 ). We chose this algorithm because the decoders so constructed would be interpretable by visualizing the weighted features in the trained models, and because their ability to generalize given a dataset from minimal training would be suitable for EEG-based decoding ( Lotte et al., 2018 ;Zhang et al., 2017 ). RF generates an ensemble of individual decision trees that define the decision boundary in a limited range of feature space to fit subsampled datasets, and the ensemble then determines the prediction based on the majority. In our study, the number of decision trees in the model was set at 100, and a curvature algorithm was exploited in the construction of the decision trees ( Loh, 2002 ). Such an ensemble classifier decreases the variance of models through combining classifiers constructed from different subsampled datasets. Thus, RF is less sensitive to the effects of multicollinearity and nonlinearity of EEG-datasets compared with other algorithms, such as logistic regression or linear discriminant analyses.
There are a couple of other standard algorithms that are also interpretable and can be generalized (e.g., logistic regression, partial least squares discriminant analysis); however, they were not suitable for the present study because the assumption that the acquired datasets are linearly separable is uncertain. Typical options for EEG-decoding (e.g., support vector machine, filter bank common spatial pattern) were not suitable, as the low interpretability of the models blurs the contributions of EEG features to decoding the tasks ( Ang et al., 2008 ;Liao et al., 2014 ;Pfurtscheller et al., 2006 ).

Evaluation of decoding performance
Decoder performance was evaluated by repeated 10-fold crossvalidation. The data obtained from 160 trials were partitioned into ten subsets containing equal numbers of datasets with four "movement " types: relaxation, extension, flexion, and co-contraction. In the first stage, datasets were decoded by Dec-1 as relaxation or contraction. If a trial was decoded as a contraction, it was further decoded by Dec-2 according to contraction type. The types resulting in the highest probability were used to evaluate decoding performance. The confusion matrix for each fold was calculated and averaged after the repetitive crossvalidation test was carried out. The F-measure was used for performance evaluation because the harmonic mean of precision and recall reflects the comprehensive performance of multi-class decoders ( Chinchor and Nancy, 1992 ;Roy et al., 2019 ). The convergence of the F-measure was confirmed in the range of three significant figures following ten repetitions. Performances of both the Decoder-All and Decoder-Ctrl were also evaluated and averaged F-measure values derived from all participants were compared via the Student's paired t-test.

Calculation of feature-importance scores
To assess the contribution of individual EEG features to decoding performance, we calculated feature-importance scores. The calculation protocol was based on the "out-of-bag" predictor importance ( Breiman, 2001 ). For datasets not used as training data, the original accuracy and its change due to random permutations of each feature across datasets were calculated. Because performance was evaluated by 10 repetitions of 10-fold cross-validation, the feature-importance scores were calculated 100 times. z-score normalized scores were averaged to calculate average importance scores and the same process was applied to scores acquired from all participants. Because the featureimportance score reflects the difference in the decoder performance between the actual and randomized dataset, we identified clusters of features that significantly contributed to the decoding performance. The feature-importance score averaged across participants were assessed using the permutation test according to the following procedure.
The number of spatially neighboring features (a cluster) that scored more than 1.96 ( p < 0.05) were counted for each frequency band based on the adjacency of electrodes. Randomization of scores was applied to measure the probability at least one cluster consisted of the same number of electrodes appeared. The maximum size of the cluster in the randomized data set was counted 5000 times iteratively in each frequency band, and then the permutation p -value were calculated to test the statistical significance of the original cluster (i.e. the proportion of cluster sizes larger than the original in the randomized data; ( Nichols and Holms, 2002 )).
The same procedure was applied to scores calculated for both Dec-1 and Dec-2, and a Bonferroni correction was applied to correct for multiple comparisons (i.e. p = 0.05/5 was set as the significance threshold as the tests were individually conducted for the five frequency bands). Note that the temporal or frequency adjacency was not considered in the non-parametric permutation test as the decoders were trained with band-averaged features but not with the temporal information.

Temporal analysis of EEG signals
Typical EEG signals recorded over the contralateral SM1 (i.e., C3 channel) are shown in Fig. 3 A. SMR was observed during the "Rest" period and the power attenuation occurred during the "Task" period ( Fig. 3 B). Prominent power spectrum density peaks of the EEG signals recorded in channel C3 ( Fig. 3 C) were observed around the high-alpha band (11-13 Hz). The group time-frequency analysis from channel C3 ( Fig. 3 D) showed that task-related modulation of EEG amplitude occurred during the three contraction tasks. Non-parametric permutation test revealed statistically significant modulation of frequency power in 8-30 Hz ( Fig. 3 E).

Spatial analysis of EEG signals
Topographic maps derived from whole-head scalp EEG ( Fig. 4 A) showed that SMR-ERDs were observed only around the bilateral SM1 while those in contralateral hemisphere only exceeded statistical significance ( Fig. 4 B). sLoreta analysis in the high-alpha band-power (11-13 Hz) ( Fig. 4 C), confirmed that the difference between relaxation and contraction tasks were found in the contralateral SM1 and were spatially overlapping for all three contraction tasks (maximum at X = − 40, Y = − 25, Z = 40, best matches contralateral postcentral gyrus). Differences at the significance level of Bonferroni corrected ( p < 0.01) were only confirmed in the high-alpha band.

Decoding performance
The confusion matrices of Decoder-All and Decoder-Ctrl, averaged across the participants are shown in Fig. 5 A. Row labels indicate truth and column labels indicate prediction. Averaged F-measure values for Decoder-All and Decoder-Ctrl were 0.79 and 0.66 for Dec-1 and 0.44 and 0.36 for Dec-2. Decoder-All decoded the datasets more precisely than Decoder-Ctrl by approximately 1.2-fold. Comparisons of F-measure values in the two conditions ( Fig. 5 B) confirmed significant differences between Decoder-All and Decoder-Ctrl (Cohen's d = 1.83. p < 0.001, paired t -test). The effect size is considered to be large ( Cohen, 1988 ).

Spatial patterns of feature-importance score
The number of electrodes in the original cluster and p -values for the permutation tests of feature importance scores are shown in Table 1 . Topographic maps of clusters exceeded statistical significance ( p < 0.05/5, Bonferroni correction for multiple testing) were shown in Figs. 6 A,  Fig. 3. Results of temporal analysis for EEG signals. (A) A typical example of filtered EEG signals recorded over the sensorimotor cortex (SM1) contralateral to the activated fingers (left hemisphere), during a trial. The sensorimotor rhythm (SMR) was recorded during the rest period indicated with the black bar. (B) Arch-shaped SMR measured from 2 to 3 s at the trial shown in Fig. 3 A. (C) The power spectrum density (PSD) elicited from each participant are indicated by gray lines and averaged data across all participants are indicated by the black line. The frequency bands exploited for decoder construction are indicated above as black bars. (D) Results of time-frequency analyses with short-time Fourier transform, elicited from the EEG signals over the contralateral SM1 (i.e. C3 channel) averaged across all participants. Prominent decreases in signal strength were observed during the three contraction tasks. In contrast, signal strength was constant during the relaxation task. (E) Results of non-parametric permutation test for the time-frequency representation in C3 channel. Significant differences between the relaxation and three contraction tasks were confirmed in 8-30 Hz. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) B for Dec-1 and Dec-2, respectively (Note that only frequency-band whose cluster exceeded significance threshold were visualized). Together, these data show that the cluster of electrodes recording significantly altered EEG activity in high-alpha (11-13 Hz) and low-beta (14-18 Hz) for Dec-1 and in high-alpha for Dec-2 exceeded statistical significance. The spatial patterns derived from Dec-1 showed that EEG features in the high-alpha and low-beta bands acquired around the contralateral SM1 contributed to the decoding of relaxation or contraction. By contrast, during Dec-2, EEG features in the high-alpha band around the ipsilateral SM1 contributed to the decoding of the three different contraction patterns. Note that there were no obvious contributions to decoding contraction types of EEG features in the contralateral SM1, which contributed to the decoding of relaxation or contraction by Dec-1.

Discussion
We decoded the contraction patterns of finger muscles from EEG features using an RF algorithm in two-stages. Our results indicated that (1) Decoder-Ctrl, constructed with EEG features derived only from electrodes around the contralateral SM1, was not sufficient to dissociate the various contraction patterns despite previous functional magnetic resonance image (fMRI) studies exhibit contralateral SM1 represent the property of finger movements ( Berlot et al., 2019 ;Diedrichsen et al., 2013 ); (2) Decoder-All, constructed with EEG features derived from the whole head, showed significantly better performance than Decoder-Ctrl; and (3) analysis of feature-importance scores revealed that the EEG features derived from electrodes around the ipsilateral SM1 contributed to decoding of contraction patterns while those around contralateral SM1 contributed to decoding the distinction between relaxation and contraction.

Comparison of Decoder performances
Decoder-All outperformed Decoder-Ctrl constructed with the same procedure yet using features from the electrode around contralateral SM1. In spite of the spatial limitations of EEG signals and their relatively lower signal-to-noise ratio, high-density EEG signals have been exploited to decode the properties of finger movements ( Alazrai et al., 2019 ;Hayashi et al., 2019 ;Iturrate et al., 2018 ;Liao et al., 2014 ;Quandt et al., 2012 ). For instance, using a similar high-density EEG system, Hayashi et al. (2019) reconstructed produced forces of unilat- Fig. 4. Results of spatial analysis for EEG signals. (A) Topographic map of EEG signaling changes during the task period. During contraction tasks, power attenuations of EEG signals were observed only around electrodes over the bilateral SM1. (B) Results of non-parametric permutation test for topographic maps of EEG signals during the task period in 8-30 Hz. Statistical significance was observed around contralateral SM1 but not in ipsilateral. (C) sLoreta analyses of statistical non-parametric mapping for estimated cortical sources of band power in the high-alpha band (11-13 Hz), which exceeded statistical threshold corrected for multiple comparisons. Significant differences between resting and task periods were observed for contraction tasks. Differences were found around the contralateral SM1 cortices for all three contraction tasks ( p < 0.05, Bonferroni correction). Derived masks superimposed on a standard brain template were visualized by MRIcroGL ( https://www. mccauslandcenter.sc.edu/mricrogl/home ). Each mask represents voxels that significantly exceeded resting activity ( p < 0.05, Bonferroni correction). eral extensor muscles and Liao et al. (2014) classified activated finger digits into two classes. However, to the best of our knowledge, there are no reports of the decoding of simple muscle contraction patterns at the same force level but with different effectors in multi-class classification tasks, which is important for naturalistic control of BCIs with high degree of freedom. The performance of Decoder-All averaged across the participants was approximately two times higher than chance-level performance. This is comparable with previous reports that examined multi-class classification of movements within a limb ( Alazrai et al., 2019 ;Quandt et al., 2012 ). To decode extensor and flexor activities, current source estimation (CSE) is thought to be necessary to extract spatially enhanced cortical activity from scalp EEG signals ( Edelman et al.,  Fig. 6. Topographic maps of feature-importance score. (A) The results of permutation tests for the spatial extent of feature-importance scores elicited from Dec-1, which dissociates relaxation and contraction. Statistically significant clusters in high-alpha and low-beta band are shown. (B) The result of permutation tests for the spatial extent of feature-importance scores elicited from Dec-2, which dissociates contraction types (i.e. extension/flexion/co-contraction). A statistically significant cluster in the high-alpha band is shown. Results indicate that decoding of relaxation/contraction was mainly realized by EEG features over SM1 contralateral to the examined fingers, while the decoding properties of the three contraction tasks were mainly realized by EEG features over the ipsilateral SM1.
2019 ; Mattout et al., 2005 ;Pascual-Marqui, 2002 ). For instance, with a CSE method using fMRI signals during execution as a prior information source, Yoshimura et al. (2012) successfully decoded extensor and flexor activities using 32 EEG electrodes. Because the high-density EEG system augments the spatial frequency of EEG, modeling of cortical dynamics can be achieved by two-stage decoding without prior assumption of CSE algorithms ( Robinson et al., 2017 ).

Contribution of features around contralateral SM1 to performance of Dec -1
The topographic maps derived from the feature-importance scores of Dec-1 revealed that the electrode cluster around the contralateral SM1 in the high-alpha (11-13 Hz) and low-beta (14-18 Hz) bands contributed to the decoding of relaxation/contraction. The contribution of contralateral SM1 is consistent with previous reports, as unilateral finger movements induce stronger SMR-ERD in the contralateral SM1 than in the ipsilateral ( Neuper and Pfurtscheller, 2001 ). This laterality was also confirmed by our sLoreta analysis, as voxels significantly exceeding resting level were located around the contralateral SM1 in the high-alpha band. The sLoreta analysis in the middle and high beta-band did not reach statistical significance due to the inter-subject variability of reactive frequency beta band ( Espenhahn et al., 2017 ), while significant differences were observed in time-frequency analysis using ERSP values normalized in a trial by trial manner.
For SMR signals in the low-beta band, significant differences were not found by the sLoreta analysis, yet the time-frequency map derived from the C3 channel ( Fig. 3 E) indicates that EEG signals in low-beta were partially attenuated during the contraction tasks. Espenhahn et al. (2017) reported that such responses in the beta band exhibit high intra-individual reliability. Because such EEG features with low variability contributed to single-trial decoding of relaxation and contraction, the importance of features in the low-beta band over the contralateral SM1 were emphasized. As such task-related modulation of SMR amplitude is a robustly observed phenomenon ( Debnath et al., 2019 ;Jurkiewicz et al., 2006 ;Miller et al., 2010 ;Pfurtscheller and Lopes Da Silva, 1999 ;van Wijk et al., 2012 ), this result corroborated that the observation of important EEG features using the RF algorithm can be utilized to find task-related changes in EEG features.

Contribution of EEG features around ipsilateral SM1 to performance of Dec -2
Our results indicated that EEG features in the high-alpha band elicited from the ipsilateral SM1 mainly contributed to the performance of Dec-2, which decoded three contraction patterns of hand muscles. Even though the innervation of the hand muscles by the ipsilateral hemisphere is limited ( Gandevia and Burke, 1988 ;Soteropoulos et al., 2011 ), substantial number of neurophysiological TMS studies reported involvement of ipsilateral M1 (iM1) in unilateral voluntary movements ( Davare et al., 2007 ;Duque et al., 2008 ;Hübers et al., 2008 ;Ki či ć et al., 2008 ;Morishita et al., 2011 ;Uehara et al., 2013 ). Davare et al. (2007) reported that single-pulse TMS applied to the iM1 induced delay in proper recruitment of agonist and antagonist muscles during fine motor tasks. As our datasets were acquired under the feedback of iEMG signals and participants were directed to control their force to match the target on display, constant regulation of agonist and antagonist muscles was required. Thereby, our results suggest that the contraction pattern-specific changes in EEG-features around the ipsilateral SM1 reflect such modulatory contributions to motor control, which are likely realized by the cortico-cortical interactions between bilateral M1 via the corpus callosum ( Ferbert et al., 1992 ;Verstynen et al., 2005 ). This explanation is corroborated by reports that cathodal transcranial direct current stimulation applied over ipsilateral M1 improved behavioral results in unilateral motor tasks that require muscle-selective activation of agonist and antagonist muscles ( McCambridge et al., 2011 ;Uehara et al., 2015 ). Moreover, such modulation of M1 excitability is volitionally induced through regulating intracortical inhibitory circuits in contralateral (cM1), as the power attenuation of SMR induced by motor imagery reflects the disinhibition of the short-latency intracortical inhibition in a muscle-selective manner ( Takemi et al., 2018 ). Not only intracortical disinhibition, but also the strength of short-interval interhemispheric inhibition from cM1 to iM1 can be probed with the phase synchronicity and amplitude of SMR signals in the frequency range of the alpha-band, recorded over bilateral M1 ( Stefanou et al., 2018 ). These studies indicate SMR contains information about both the intra-and inter-hemispheric inhibitory activity of M1, that plays an essential role to control contralateral corticomotor excitability in a muscle-selective manner. Therefore, the weighted EEG features around ipsilateral SM1 putatively contributed to decoding of the contraction patterns of finger muscles because they represent the activities of such inhibitory circuits. Future studies testing the relationships between bilateral SMR signals and interhemispheric inhibition during motor tasks would elucidate the characteristics of iM1 activity represented in EEG signals.
One might say that contraction pattern-specific activations of ipsilateral uncrossed pathway contributed to the decoding since uncrossed corticospinal pathways could influence the EEG signals around ipsilateral SM1. Indeed, the ipsilateral SMR-ERD reflects the excitability of uncrossed corticospinal pathway ( Hasegawa et al., 2017 ). However, such influence is less likely in the present study because the distal mus-cles (EDC and FDS) are exclusively innervated by the contralateral SM1 ( Soteropoulos et al., 2011 ).

Implications for BCI
In the present study, we constructed decoders and evaluated them in offline analyses; however, to apply them to BCIs, these analyses need to be implemented online. Generally, there are two issues to be addressed in achieving online processing. The first is the computational requirement for online processing and the second is the ability of the constructed decoders to generalize to the inter-day variability of recorded EEG signals or to the brain state of users. The former issue is readily resolved because the RF classifier we adopted uses an ensemble of decision trees, which does not require much processing power because the output is determined by comparing feature values with tuned thresholds node by node. Hence, as far as the computational power is concerned, the online decoding of finger movements is possible once the decoder is constructed offline. Regarding the latter issue, it is necessary to confirm whether the constructed decoders can be applied to new data obtained from the same subject on different days. However, as previous BCI studies reported, the calibration conducted at the beginning of a BCI run is an effective method to achieve inter-day generalization by constructing robust decoders. Moreover, the closed-loop feedback induces learning in BCI users ( Benabid et al., 2019 ;Perdikis et al., 2018 ;Shenoy and Carmena, 2014 ).
Along with the processing speed and accuracy, the degree of freedom would be an important factor for naturalistic control of BCIs ( Schwarz et al., 2019 ). Thus, studies that test whether it is feasible to expand the two-stage decoding method into datasets during both unimanual and bimanual movements are warranted.

Conclusions
In the present study, we tested whether the contraction patterns of agonist and antagonist muscles of the same finger could be decoded from whole-head high-density EEG signals using a two-stage decoding method, and we examined which EEG features contributed to decoder performance. Decoding performance was significantly higher using features obtained from multiple electrodes than using features derived only from the SM1 contralateral to the finger. The features derived from the contralateral SM1 played an important role in decoding Relax/Contraction (i.e., to predict the presence of finger movement). By contrast, the SM1 ipsilateral to the finger muscles was important for decoding the type of contraction: Ext/Flx/Cocont (i.e. to predict the detailed activities of extensor and flexor muscles in the fingers). Our results indicate that the detailed state of the sensorimotor system can be predicted by the SMR recorded from bilateral SM1 by decoders treating multivariate EEG-features, despite conventional univariate approach did not effectively yield the information regarding the coordination of finger muscles. This finding should play a significant role in deciding the therapeutic cortical target for novel clinical treatments (e.g., neurofeedback, repetitive transcranial magnetic stimulation, and transcranial electrical stimulation), especially for motor deficits of stroke patients characterized by spasticity, which induces abnormal co-activation of extensor and flexor muscles caused by a maladaptive sensorimotor system ( Grefkes et al., 2008 ;Li, 2017 ;Li and Francisco, 2015 ;Pool et al., 2018 ).

Declaration of Competing Interest
J.U. is one of the founders and the Representative Director of the University Startup Company, Connect Inc. for the research, development, and sales of rehabilitation devices including brain-computer interface. J.U. has received a salary from Connect Inc., and has held the shares of Connect Inc. This company does not have any relationship with the present study.