2.1 Experimental Procedure and Data Acquisition
We recruited a total of 44 healthy subjects for this study, and assigned them randomly to three groups—sham (control), vibrotactile stimulation, and tDCS. This experiment was conducted in a quiet space, and the subjects were seated in a comfortable chair approximately one meter from a 24-inch screen. The subjects were asked to place both of their hands on a desk to prevent hand movement while they performed the motor imagery tasks. We used OpenViBE software  and custom-built MATLAB scripts to acquire EEG, process online signals, and present the motor imagery task. As shown in Figure 1, EEG were recorded first for one minute during the eyes-open resting state. Thereafter, the subjects performed a left- and right-hand motor imagery task that consisted of two offline and two online blocks each. The subjects received different stimuli while performing the four offline motor imagery blocks according to their assigned groups. After the stimulation session, they performed the same motor imagery task that they performed before the stimulation session, and eyes-open resting state EEG were recorded before and after the stimulation session. EEG were recorded from 19 wired dry electrode channels, referenced by the left and right earlobes at 300Hz (DSI-24, Wearable Sensing, USA). This experiment was approved by the Institutional Review Board at Gwangju Institute of Science and Technology (20210806-HR-62-03-02), and all subjects were informed about the experimental procedure and signed informed consent.
In the motor imagery task, the subjects performed left- and right-hand motor imagery rather than actual movement. The subjects were asked to imagine hand movements in a kinesthetic rather than visual way, as a previous study reported that kinesthetic motor imagery (KMI) and visual motor imagery (VMI) showed different brain activations . Each trial began with a fixation cross for the first 1.5 seconds. During the next 4 seconds, an orange-colored circle appeared in the center of the screen, and an instruction bar appeared on the left or right side. The subjects were instructed to keep imagining left- and right-hand movements according to the direction of the bar until the blank screen, which remained blank for 1.5 seconds. This trial was repeated 40 times and included shuffled left- and right-hand trials for each block. During the offline phase (two blocks), there was no motor imagery feedback. After the offline phase, the subjects’ EEG data were used to train common spatial pattern (CSP)-based spatial filter and Fisher’s linear discriminant analysis (FLDA) classifier. Thereafter, the subjects performed the online phase that consisted of two motor imagery blocks that showed visual feedback for 1.5 seconds by moving the circle to classified directions.
For CSP-FLDA, an epoch was extracted from each trial at [1000-3500] ms to the stimulus onset, and band-pass filtered with cutoff frequencies of 8 and 30Hz using the 4th-order Butterworth filter. However, we note that hyper-parameters for online classifiers were changed after the first ten subjects. Specifically, for the first ten subjects, CSP-FLDA was trained by 19 electrode channels and the first and last two CSP filters were selected. However, electrodes around the eyes and occipital areas often yielded notable noise with high variation attributable to non-brain activity, such as motion artifacts, eye movements, and bad contacts. As a result, for the remaining subjects, nine electrode channels around the central area, including the F3, Fz, F4, C3, Cz, C4, P3, Pz, and P4 channels, were used to train CSP, and the first and last CSP filters were selected to train FLDA. In this study, we calculated offline BCI performance to compare BCI performance with the same hyper-parameter over all subjects, as described in section 2.3.
2.2 Stimulation Session Design
In this study, we compared the stimulation effects on brain activity modulation and relevant BCI performance. Depending upon the stimulation group, each subject performed a stimulation session in the middle of two motor imagery sessions. The stimulation groups consisted of the vibrotactile stimulation group, tDCS group, and sham stimulation group. Figure 4.2 illustrates the stimulation session for each group. All stimulation group subjects performed four blocks of offline motor imagery task while receiving stimulation.
With respect to vibrotactile stimulation (Figure 4.2A), two vibration motors (Model 310-113, Precision Microdrives, England) activated by a Arduino Due board were used to deliver vibrotactile stimulation to the left and right index fingertips. Each motor is 10 mm in diameter with a 1.34G vibration amplitude. The Arduino board was triggered through serial communication with MATLAB scripts for closed loop vibrotactile stimulation. Recent studies have shown that EEG-guided stimulation can enhance the stimulation effects in vibrotactile stimulation and transcranial magnetic stimulation (TMS) by targeting motor cortex excitability states [20, 35]. Researchers have observed that stimulation triggered by the EEG alpha falling phase outperformed continuous stimulation with respect to the motor evoked potential (MEP)  amplitude and motor imagery BCI performance . In this respect, our study applied vibrotactile stimulation according to the EEG alpha (8 ~ 13Hz) phase at the electrode channels on the left (C3) and right motor cortex (C4) as in . Specifically, during the vibrotactile stimulation session, the subjects were instructed to place their left and right index fingertips on each vibration motor. For each trial, the epoch was extracted from the contralateral EEG channel as often as 500ms every 50ms to calculate the alpha phase. The epoch extracted was band-pass filtered in the [8-13]Hz frequency range using a 10th-order elliptical infinite impulse response (IIR) filter, and the phase was calculated using the Fast Fourier Transform (FFT)-based phase tracking algorithm, as proposed and adopted previously [20, 36]; among the FFT amplitudes, the dominant alpha frequency component between 8 and 13Hz and the corresponding phase were used to obtain a simple sine function to predict the upcoming phase. When the phase predicted was falling, the vibration was delivered through the left or right vibration motor for 100ms according to the motor imagery class, and the inter-stimulation interval was set to 100ms. Thus, the C4 channel alpha phase was extracted for the left-hand motor imagery trials, and the left vibration motor was activated when the phase predicted was falling and the converse. After the stimulation session, the subjects performed the same motor imagery task that they did before the stimulation session.
For brain stimulation (Figure 2B), high-definition transcranial direct current stimulation (HD-tDCS) was delivered to the motor cortex using one anode electrode and four neighbouring cathode electrodes (Starstim8, Neuroelectrics, Spain) after taking off the EEG cap. An anode electrode was placed to stimulate the contralateral motor cortex of the non-dominant hand; because all subjects in the tDCS group were right-handed, the anode electrode was placed on C4 and the cathode electrodes were placed on FC2, FC6, CP2, and CP6. To minimize pain attributable to tDCS, sufficient gel was injected, and the impedance level was maintained below 2K throughout the stimulation session. The stimulation intensity reached a target intensity of 2mA during a 30-second ramping period. The stimulation session continued until the subjects performed four blocks of the offline motor imagery task, although we did not record EEG during the tDCS session because of electrical interference. The subjects were informed that they could stop the stimulation at any time if they felt severe pain. In addition, the sham stimulation group (control group) subjects performed the same task as the tDCS group subjects, but they received sham, rather than real stimulation. For the sham stimulation group, the stimulation turned off after the ramping period, and the subjects did not know whether they belonged to the tDCS or sham stimulation group. After the stimulation session, the subjects washed their hair and were re-equipped with the EEG cap. Then they performed the same motor imagery task as they did before the stimulation session.
2.3 BCI Performance Evaluation
This study used offline BCI performance because online BCI classification parameters were not the same for all subjects. To evaluate the subjects’ BCI performance, the EEG data were band-pass filtered with [8-30]Hz using the 4th-order Butterworth filter, and any 60Hz line noise that remained was filtered out with the band-stop filter from [58-62]Hz. Epochs were extracted from [500-3500] ms to the stimulus onset, which is the time window obtained heuristically. Except for the electrodes near the eyes, those with an amplitude greater than ±100µV were removed, and subjects who had more than 30% of bad trials among all trials were eliminated from the analysis. Finally, we chose a region of interest (ROI), the nine electrode channels (F3, Fz, F4, C3, Cz, C4, P3, Pz, and P4) for motor imagery BCI classification using visual inspection. To evaluate BCI performance, the first two motor imagery blocks (offline phase) were used for training, and the remaining two blocks were used for testing. The Riemannian minimum distance metric (MDM)  was used to extract features and evaluate BCI performance.
In this study, the subjects performed the motor imagery task before and after the stimulation session to evaluate stimulation effects on brain activity during motor imagery. Moreover, we divided each stimulation group into low- and high-performing groups according to their pre-stimulation BCI performance as previous studies have observed that low and high BCI performers showed different neurophysiological characteristics [9, 30]. Further, our previous study suggested that low and high BCI performers should be treated differently because low BCI performers’ features decreased the ability to generalize the cross-subject BCI model significantly, while subject selection may have increased cross-subject BCI performance . In addition, we assumed that poor BCI performers may have greater potential to change their brain activity through learning or external stimulation, for either better or worse, compared to high BCI performers because high BCI performers’ brain activity may be stable and optimal already. To divide the subjects into low and high BCI performance groups, a previous study used the median value of the BCI performance , as the median value can divide groups with the same size and allow analysis to be performed on sub-groups of the same size. However, the distribution and sample size affect the median value. Accordingly, if the BCI performance data collected are biased toward high or low performance, the divided groups do not represent true high or low performance groups. Instead, we used the statistical random probability introduced in  and used in our previous study  to divide low and high BCI performers, as it can produce a random threshold based upon the number of trials and the statistical significance. With 40 test trials and α=0.05 significance level, we obtained 60.69% as a threshold to divide low and high BCI performers. As a result, we divided each stimulation group’s subjects into low and high BCI performers based upon their pre-stimulation performance and investigated the stimulation effects in the subgroups.
2.4 Pre-stimulus Band Power for Motor Imagery
The pre-stimulus band power was calculated from the online phase of the motor imagery task before and after the stimulation session. Maeder et al. investigated the relation between the pre-stimulus sensory motor rhythm (SMR) band power and motor imagery BCI performance and observed that higher pre-stimulus SMR trials yielded significantly better performance compared to lower trials . In this respect, we compared the pre-stimulus band power before and after the stimulation session for each group to investigate the stimulus effects on the pre-stimulus band power. The online motor imagery task phases before and after the stimulation session were used to calculate and compare the pre-stimulus band power. To perform pre-processing, the epochs were extracted first from the EEG data from [-1000-4000] ms to the stimulus onset, band-pass filtered with [1-40]Hz, and any epochs larger than ±100µV were removed, except for the electrodes near the eyes, and the same bad subject criterion (>30% bad trials) was applied as in the BCI performance evaluation.
After extracting the epochs and eliminating the bad trials, we calculated the pre-stimulus band power as the logarithm of the band power during 1000ms preceding the stimulus onset from the epochs extracted after the band-pass filtering with alpha (8-13Hz), low-beta (13-20Hz), and high-beta (20-30Hz). As the previous study selected an electrode channel from the left and right hemisphere and averaged over those electrode channels , we obtained the average pre-stimulus band power averaged over the C3 (left hemisphere) and C4 (right hemisphere) electrode channels. In addition, we obtained pre-stimulus band powers around the nine electrode channels (F3, Fz, F4, C3, Cz, C4, P3, Pz, and P4) used to evaluate motor imagery BCI performance. To evaluate statistical significance, pre-stimulus band powers before and after the stimulation session were compared using a paired Student’s t-test, and the Bonferroni correction was applied for the p-values obtained from the nine electrode channels.
2.5 Functional Connectivity during Motor Imagery
The functional connectivity between different brain regions in sensor (electrode) space can be assessed using phase synchronization. Through functional connectivity, the way the cortical regions communicate with each other and the way the information is transmitted between different regions during a cognitive task can be understood [29, 39]. One way to measure phase synchronization is by assessing phase-locking, which denotes a constant phase difference between two signals that remains constant for a certain period . In this study, a phase locking value (PLV) between the EEG channel was calculated with the following equation introduced in :
In which t stands for an extracted epoch (second) for each trial n up to N, N is the total number of trials, and the exponential term indicates the phase difference between two signals in the same trial, which denotes the difference between two phases extracted from the two signals [39, 40]. As a result, the PLV can measure the variability of the phase difference across trials; when the phase difference is small, the PLV is close to 1, and the PLVS is near 0 otherwise. In this study, the online motor imagery task phases before and after the stimulation session were used to calculate and compare the PLV. Specifically, the epochs were extracted from the EEG data from as long as [500-3500] ms to the stimulus onset, band-pass filtered with [8-30]Hz, and the same trial rejection and subject rejection criteria were applied as in BCI performance evaluation. With respect to scales of connectivity, this study employed a global PLV that can be obtained by the average connectivity over all electrodes around the central area (Figure 3), as investigated in . The nine electrode channels (F3, Fz, F4, C3, Cz, C4, P3, Pz, and P4) used to evaluate BCI performance and pre-stimulus band power were selected. In addition, we investigated the PLV with broader frequency bands including alpha (8-13Hz), low-beta (13-20Hz), and high-beta (20-30Hz) bands by calculating the PLV within these frequency bands . The PLV was calculated by MATLAB scripts using the FieldTrip toolbox .
In this study, the PLVs were compared before and after the stimulation session using a paired Student’s t-test for subjects with low and high BCI performance in the sham, vibrotactile stimulation, and tDCS groups to assess the stimulation effects on the average strength of the connection around the central areas during motor imagery.