Assessing the impact of vibrotactile kinaesthetic feedback on electroencephalographic signals in a center-out task

Objective. An important part of restoring motor control via a brain-computer interface is to close the sensorimotor feedback loop. As part of our investigations into vibrotactile kinaesthetic feedback of arm movements, we studied electroencephalographic signals in the δ, µ and β bands obtained during a center-out movement task with four conditions: movement with real-time kinaesthetic feedback, movement with static vibrations, movement without vibrotactile input, and no movement with sham feedback. Approach. Participants performed center-out movements with their palm on a flat table surface. One of three movement directions was cued visually before the movement. The palm position was tracked in order to provide real-time vibrotactile feedback. All analyses were performed offline. Main results. Movement-related cortical potentials exhibit minor discrepancies between movement conditions as well as between movement directions, in peak amplitude and shape. Classification of each movement condition and each direction against rest yields peak accuracies of 60%–65% using low-frequency amplitude features, and 90% using µ and β power features. Within-class accuracies of four-way classification between conditions based on low-frequency amplitude features are around chance level for the movement conditions with vibrotactile stimulation, slightly above chance level for the movement condition without stimulation, and considerably higher for the non-movement condition. Four-way classification between conditions based on µ and β power features yields within-class accuracies slightly above chance level for all movement conditions, and considerably higher for the non-movement condition. Within-class accuracies of three-way classification between directions are slightly above chance level for low-frequency amplitude features, and at chance level for power features. Significance. We found that the vibrotactile stimulation does not interfere with movement-related features in the δ, µ and β frequency ranges. Our feedback system may therefore feasibly be deployed in conjunction with a BCI based on movement-related cortical potentials or sensorimotor rhythms, without adversely affecting control.


Introduction
A substantial part of research in the field of braincomputer interfaces (BCI) concerns itself with ways to restore lost or degraded motor functionality for persons affected by traumatic injuries like spinal cord injury (SCI) or neurological disorders, in order to enhance the autonomy and life quality of these end users. Intact motor function operates akin 1 Author to whom any correspondence should be addressed.
to a fine-tuned control loop, where feed-forward processes are adapted in real-time, based on sensory inputs.
A large part of state-of-the-art research geared towards restoring motor control focuses on the feedforward aspects, usually closing the feedback loop with visual feedback (artificial or natural). The visual system certainly is very apt at processing dense information loads, however tying up the visual sense in observing the end effector may not be desirable in practical use cases. Moreover, there is evidence that the removal of kinaesthetic or haptic feedback negatively influences the execution of motor control [1] as well as cortical movement-related signals [2]. Apart from the performance aspect, somatosensory feedback contributes to the successful integration of a prosthesis as an extension of a user's body rather than a foreign object. According to a survey conducted by Bidiss et al [3], an absence of tactile feedback was frequently stated as a reason when prostheses are abandoned by users. Several further publications based on surveys with endusers [4][5][6] also emphasize the desire for tactile feedback.
Especially when an end effector is controlled via a BCI, feedback strategies have to be carefully designed and evaluated, so as not to interfere with control signals, either by introducing artifactual contamination, or by influencing the brain activity with the additional somatosensory input to an extent where control performance would be affected.
For a long time, the dominating control parameters for electroencephalography (EEG)-based BCIs for motor control have been derived from movementrelated modulations of sensorimotor rhythms [7][8][9][10][11]. This approach can provide a limited number of discrete control commands that reflect sustained activity, and generally require non-intuitive mappings between tasks and control functions. Recently, the focus is increasingly shifting towards low-frequency EEG, in the interest of developing more intuitive control strategies. One signal of interest is the movement-related cortical potential (MRCP), which manifests as a negative amplitude deflection in the delta frequency range when a movement is initiated [12][13][14]. The shape of an MRCP is believed to encode a number of characteristics of the respective performed, imagined or attempted movement. For instance, Slobounov et al [15] and Jochumsen et al [16] showed that the force of the movement is reflected in the MRCP. Furthermore, discriminability of different same-limb movements in the δ band has been demonstrated in several cases [17][18][19][20].
For sensorimotor rhythm-based control strategies, artificial somatosensory feedback of different modalities has been documented to have a beneficial effect [21,22]. Specifically, an enhanced power decrease in the mu frequency range has been found for motor imagery with proprioceptive feedback [21], and neuromuscular stimulation during motor imagery [22]. However, Hommelsen et al [23] raised concerns that the additional afferent input could increase false positives in sensorimotor rhythm-based BCIs since it elicits similar activation patterns to motor tasks in the mu band. With respect to low-frequency EEG, Mrachacz-Kersting et al [24] reported mixed results on whether tactile stimulation training in the form of deep pressure application improves the variability of movementrelated cortical potentials (MRCPs) in stroke patients.
In previous works employing artificial somatosensory feedback, stimulation has commonly been realised using mechanotactile [25,26], electrotactile [27,28], or vibrotactile modalities [29][30][31][32]. Vibrotactile feedback has been employed in several contexts, as force feedback [29], as coded patterns to provide feedback for discrete BCI applications [30,31], or to evoke sensations of movement on the skin [32]. Here, we focus on the latter, in order to provide kinaesthetic feedback of arm movements. The evocation of moving sensations relies on tactile illusions that arise as a consequence of inaccurate tactile perception, especially in loci of the human body where receptor density is low. More specifically, an apparent tactile motion, i.e. a perceived moving stimulus, is created when two stimuli are presented with a temporal offset, such that they are too close together to be perceived as separate from each other, yet too far apart to perceptually collapse into one merged stimulus [32][33][34]. Additionally, when two stimuli are active simultaneously and too close spatially for them to be identified as separate stimuli, they will be perceived as a single stimulus (a 'virtual stimulus' or 'phantom sensation' [35]) at a location in between the physical stimuli, where the location of the virtual stimulus depends on the relative intensities, and its perceived intensity on the absolute intensities. Several behavioural studies [32,36,37] indicate that it is most suitable to use a mapping based on a power relation between the desired location of the virtual stimulus and the intensities of the stimuli.
To our knowledge, research on the ramifications of vibrotactile feedback on EEG signals is still scarce, particularly in low frequencies. In this work, we assess the impact of continuous vibrotactile feedback on a variety of movement-related parameters in EEG, which can be utilized to control an end effector with a BCI, examining both neurophysiological features and classification based on both low-frequency amplitude features, and µ and β power features. We hypothesize that vibrotactile kinaesthetic feedback can be provided without adversely affecting these parameters, and might potentially have a facilitating effect.

Participants
The study was conducted with 15 (8 male, 7 female) able-bodied participants, aged 20 to 57 (average age 26). According to self-reports, all were right-handed. Participants were briefed on the experiment, and gave written informed consent for their participation. The study protocol was approved by the ethical committee of the Medical University of Graz. Due to heavy artifact contamination that could not be mitigated by the measures described below, the data sets of three participants were excluded from the analysis.  figure 1. In each trial, participants were informed on the condition, and the movement direction by a visual cue on a monitor located approximately 120 cm away from the participants. At the start of the trial, the cue, in the form of a coloured dot, appeared in the center of the screen. The cue was green in movement trials (conditions ME+VF, ME+Vstat, ME), and grey in non-movement trials (condition VF). In trials where vibrotactile stimulation was delivered (conditions ME+VF, ME+Vstat, VF), the cue was framed by a black border. One second after its appearance, the cue proceeded to move in one of three directions: right, left or up. Three example cues are illustrated in figure 2(A).
Participants were instructed to observe the moving cue, and remember the condition and direction. After the disappearance of the cue, the vibrotactile stimulation was turned on in the central position. One second later, a fixation cross appeared on the screen. Participants were told to fixate their gaze on the fixation cross, wait for a couple of seconds, and then initiate a movement in the direction indicated by the cue. The cue moving upwards cued a movement to the front-this case will be referred to as "up" from now on. The trial was over when a fixed distance from the central position was reached. Afterwards, a break followed. The next trial began when the participant moved their hand back to the starting position.
The paradigm was designed such that signals pertaining to the initiation of the movement (e.g. movement-related cortical potentials (MRCPs)) may be separated from responses evoked by the paradigm (visual and somatosensory evoked potentials due to the visual stimuli and the stimulation onset, respectively). Participants were instructed to perform each movement in a smooth and continuous fashion. Furthermore, they were informed of the following constraints concerning the movement, which would be enforced by the paradigm. Firstly, the movement should not be initiated earlier than 2.5 seconds after the appearance of the fixation cross, in order to avoid an overlap of the response evoked by the appearance of the fixation cross with the MRCP. Secondly, the movement should last at least 2 seconds. If any of the mentioned constraints was violated, the trial was aborted. Aborted trials were repeated at a random time point during the remainder of the run.
Following general instructions, participants performed a practice run, where they were led through all possible combinations of conditions and movement directions. In the first half of the practice run, they received additional written instructions on the screen while the paradigm was running to remind them of the meaning of the cues and the timings. In the second half of the practice run, the written instructions were removed. Each participant completed six runs of 36 trials, resulting in 54 trials per condition. Additionally, they performed five runs of 1 minute continuous rest.
The hand position was tracked by a Leap Motion controller, and fed back to the feedback control system in order to produce the real-time kinaesthetic feedback. The feedback was designed to mimic the movement of the hand position along the axis of the cued movement direction. At the beginning of each of the main runs, the central position was calibrated to the actual hand position.

Vibrotactile stimulation
Vibrotactile stimuli were delivered via C-2 tactors (Engineering Acoustics Inc. Casselberry, USA) that were driven by an ARM Cortex M4 micro-controller (STMicroelectronics, Geneva, Switzerland) and custom amplifiers. The driving signal was a sinusoid with a fixed frequency of 250 Hz. The tactor intensities were controlled in Python 2.7 via a serial interface. The tactors were mounted by attaching them to the inside of an elastic shirt with velcro strips. In order to accomodate participants' physiques, shirts were provided in several sizes. For the purposes of this study, 4 tactors were arranged as shown in Figure 1. Experimental paradigm. At the beginning of each trial (t = 0), the stimulation is off, and the visual cue is presented in the middle of the screen. Depending on the condition, the cue is a green (movement conditions) or grey (non-movement condition) dot. After 1 s, the cue slowly moves to the right, to the left, or up. After the disappearance of the cue, the stimulation turns on (conditions ME+VF, ME+Vstat, VF), and a fixation cross is displayed. The movement is initiated at a time chosen by the participant, and executed as one smooth and continuous movement to the direction dictated by the cue. The trial ends when a constant predefined distance from the starting position is reached. In the beginning of each experiment, a brief calibration protocol was executed to ensure that all tactors were perceived at an equal intensity level.
The real-time kinaesthetic feedback was realised by modulating the tactor intensities to evoke a virtual stimulus at the current measured hand position, with a sample rate of 10 Hz. The position was constrained to be on the axis spanned by the two tactors corresponding to the start and end point of the cued movement. Based on [32,36], and a small pilot study we conducted [37], we employed the following mapping between the desired location of virtual stimuli and the tactor intensities that is based on the following relation, with A 1 , A 2 the tactor intensities of the start and end tactor, respectively, and x v , A v the perceived location and intensity of the virtual stimulus: (1)

Signal aquisition
EEG and EOG was recorded from 61 and 3 actiCap electrodes, respectively, using two BrainAmp amplifiers (Brain Products GmbH, Gilching, Germany), at a sampling rate of 1 kHz. Electrodes were arranged according to the international 10/20 EEG system. The ground electrode was placed at the AFz location, with the reference electrode on the left mastoid. Frontal, temporal, and parieto-occipital channels were utilized for artifact removal and topographic visualization of features, but omitted for classification. Biosignals, behavioural signals and markers were streamed via the Lab Streaming Layer (Swartz Center for Computational Neuroscience, UCSD, USA, https://github.com/sccn/labstreaminglayer) framework, both for the online computation of the vibrotactile feedback and for offline recording.

Signal processing
For offline analyses, all signals were processed in Matlab 2016b (MathWorks, Inc. Nattick, USA).
We examined the variance of the EEG channels. Since no channels exceeded a variance threshold of 5 times the within-subject median, none were rejected. The EEG signals were preprocessed separately in order to compute event-related potentials, movement-related cortical potentials, timefrequency decompositions, as well as to prepare the signals for classification. A generalized overview of the preprocessing procedure is presented in figure 3.

Artifact removal
In order to mitigate artifact contamination, the EEG and EOG signals were filtered with a 4th order zerophase Butterworth band pass filter between 1 and Steps pertaining to artifact removal are shaded in grey. For each of the conducted analyses, i.e. slow potentials, time-frequency decomposition, and preparation for classification, the raw signals were preprocessed separately, with the appropriate band pass filter and downsampling rate. Signals were epoched relative to the onset of the fixation cross for event-related potentials, and relative to the movement onset for all analyses pertaining to movement-related parameters. Trials were rejected based on their peak amplitude and kurtosis. The remaining trials were subjected to an independent component analysis, in order to identify artifactual components.
60 Hz, and downsampled to 200 Hz. Trials were epoched from the onset of the fixation cross to the end of the trial. Trials where an amplitude threshold of +/-200 muV or a kurtosis of five times the standard deviation was exceeded, were rejected. Per subject, between 0 and 5 trials were rejected.
The remaining data was imported into the EEGLAB toolbox [38], where we performed an independent component analysis (ICA) using the Infomax algorithm [39,40]. Artifactual components were identified based on SASICA [41], complemented by manual inspection of components. Between 28 and 45 independent components were marked for rejection per subject. ICA weights were stored for later rejection of components.

Event-related potentials (ERPs)
In order to identify the responses elicited by the visual and somatosensory stimuli present in the paradigm, the raw EEG was filtered with a 4th order zero-phase Butterworth band pass filter between 0.1 and 5 Hz, and downsampled to 200 Hz. Trials were epoched with alignment to the onset of the fixation cross.

Movement-related cortical potentials (MRCPs)
For analysis of movement-related cortical potentials, the raw EEG signals were filtered with a 4th order zero-phase Butterworth band pass filter between 0.1 and 3 Hz, and downsampled to 200 Hz. Trials were epoched aligned to the movement onset. To determine the movement onset, the x and y coordinate of the palm position were extracted from the Leap Motion stream, and resampled to 200 Hz. This was accomplished by fitting a polynomial to the x (right/left) or y (up) coordinate recorded by the Leap Motion in a time window of ± 2 s around the time point where the corresponding coordinate crossed a subject-specific threshold after the fixation cross onset, and finding the intersection of the tangents of the baseline and the slope after the movement onset, respectively. Afterwards, a new time window was defined between [−0.6, 0.7] relative to the computed "rough" movement onset. In this window, the baseline and slope tangents were computed again, and their intersection was defined as the final movement onset.

Sensorimotor rhythms (SMRs)
To investigate power modulations in the µ and β frequency bands, the raw EEG signals were filtered with a 4th order zero-phase Butterworth band pass filter between 1 and 60 Hz, and downsampled to 200 Hz. Trials were epoched aligned to the movement onset. Subsequently, the signals were subjected to a wavelet decomposition using Morlet wavelets [42]. The temporal resolution was set to 3 seconds (full width half maximum) at the mother wavelet's center frequency of 1 Hz. From the Wavelet amplitude, we computed the power relative to resting EEG.

Classification
Classification was performed for several trial splits. Multi-class classifiers were trained to classify between conditions, and between directions. Additionally, separate classifiers were trained to classify each condition, and each direction against rest. Trials were epoched aligned to the movement onset. Rest trials were extracted from the continuous rest runs, and preprocessed in the same manner as the other trials. We performed classifications of the above comparisons separately based on amplitude features in the δ frequency band (henceforth: 'LF classification'), and based on power features extracted from the µ and β bands (henceforth: 'power classification'), using Linear Discriminant Analysis classifiers with shrinkage regularization (sLDA). To evaluate each comparison, a leave-one-out cross validation was conducted. For low-frequency classification, the raw signals were filtered with a 4th order zero-phase Butterworth band pass filter between 0.1 and 5 z, and downsampled to 10 Hz. Amplitude features were extracted with a window length of 1 s. For classification based on µ and β band power features, the raw EEG signals were separately filtered with 4th order zero-phase Butterworth band pass filters, between 7 and 14 Hz, and between 17 and 26 Hz, respectively. Both resulting sets of filtered signals were downsampled to 100 Hz. Common spatial patterns (CSPs) [43][44][45] were computed for each 2-class combination, and the logarithmic band power of the 4 CSPs carrying the most discriminative information was averaged with a moving average with a window length of 1 s. Figure 4 shows the evoked responses time-locked to the onset of the fixation cross (t = 0). The onset of the vibrotactile stimulation at -1.12 s elicits a somatosensory evoked potential (SEP), as evidenced most notably by the central P300 peak at −0.8 s in conditions ME+VF, ME+Vstat and VF ( figure 4(A)). After the fixation cross onset, a visual evoked potential (VEP) can be observed over occipital sensors in all conditions ( figure 4(B)). A smaller response appears subsequent to the disappearance of the visual cue (t = −2.12 s).

Behavioural measures
The distributions of movement onsets and movement durations are illustrated in figure 5. In the left plots, they are sorted by conditions, while on the right, they are sorted by movement directions.
Across all trials, 95% of movement onsets happened at least 1.4 s after the fixation cross onset. The median movement onset of all trials is 3.19 s after the onset of the fixation cross, with a discrepancy between trials with stimulation (median 2.96 s/2.92 s for conditions 1 and 2, resp.) and trials without stimulation (median 3.65 s). This discrepancy is significant in ten subjects between condition 1 and 3, and in 11 subjects between condition 2 and 3, according to Wilcoxon rank sum tests at the 5% significance level. Between movement directions across all movement conditions, movements to the left had an earlier onset (median 2.9 s) than to the right or up (median 3.3 s/3.3 s), though this is only significant in four subjects each. There is no significant difference in a majority of subjects in the movement durations between trials of different conditions. The durations of movements to the right are significantly different from up movements in a slim majority of seven subjects, and from left movements in half of the subjects, while there is no significant difference in the durations between left movements and up movements in a majority of subjects (eight subjects). Figure 6 depicts the grand average MRCPs for each condition and movement direction at Cz, and in topographic plots. The left panels depict the MRCPs for all conditions, while the right panel depicts the MRCPs for the three movement directions, considering only movement trials.

Movement-related cortical potentials
There is a prominent central negativity with a slight lateralization to the contralateral hemisphere in all movement conditions, while the response in condition VF hovers around zero. The negative deflection starts between 1 and 2 seconds prior to the movement onset, reaching its peak negativity at the time of the movement onset in condition ME+Vstat, and shortly thereafter in conditions ME+VF and ME. The grand average peak negativity for movements to the right is smaller than for the other two directions. Furthermore, the shape of the decending slope in the second prior to the movement onset deviates between the three directions. Figure 7 illustrates the grand average relative power for the four conditions and three movement directions in a time window of ± 3 s around the movement onset, between 1 and 40 Hz. The conditions are shown on the left, and the movement directions on the right. For the movement directions, only movement trials are considered. The spectrum and time-frequency maps exhibit a prominent negative peak in the µ band, and a smaller peak in the β band, in all movement conditions and directions. The specific bands marked in the plots were identified based on inspections of both the grand average and individual spectra. The power decrease  in both bands is strongest over the contralateral central/centroparietal region, and weaker on the ipsilateral side. There is a small discrepancy in the strength of the spectral peaks and the localization of the power decrease between movement conditions in the grand averages. The µ peak of condition ME+VF is marginally weaker compared to the other two movement conditions, whereas the µ power decrease spreads slightly more towards parietal regions in conditions ME+VF and ME+Vstat compared to condition ME. In the non-movement condition, there is a power increase in the µ band, which is localized over parietal regions, with a faint ipsilateral lateralization, and no effect in the β band. Between directions, the spectral profile is consistent, while the strength and spatial spread of the µ peak varies slightly.

Low frequency classification
Accuracy curves for LF classification of each condition against rest are depicted in figure 8. Peak accuracies of all movement conditions are between 60% and 65%, while for trials of condition VF, the mean accuracy does not exceed the significance threshold above chance level (dependent on the number of trials, [46]). Performing analogue classification of trials split by movement directions yields highly similar results to the movement conditions. For the split by directions, only movement trials were considered. Figure 9 shows the within-class accuracy curves for multi-class LF classification between conditions and directions, respectively. Condition VF is most discriminable from the other three conditions. Condition ME rises above the significance threshold to a smaller degree in a time window around the movement onset, while the accuracies for conditions ME+VF and ME+Vstat hover around the threshold. Between directions, accuracies for directions right and left slightly exceed the significance threshold.

Power classification
Accuracy curves for power classification of each condition against rest are depicted in figure 10. At 2 seconds before the movement onset, grand-average accuracies are approximately 75% in all cases. While the accuracy of condition VF is maintained at this level throughout the whole time window of classification, accuracies of the movement conditions rise to peak accuracies of approximately 90% within the first second after the movement onset, and stay near that level. Power classification of movement directions against rest mirror the results of the movement conditions. Figure 11 shows the within-class accuracy curves for multi-class LF classification between conditions and directions, respectively. The movement conditions are discriminable at a level slightly above the significance threshold between the movement onset and the end of the time window of classification, while the non-movement condition is more distinctly discriminable for the whole time window, with accuracies of up to 73% post movement onset. Accuracies for movement directions slightly exceed the significance threshold pre movement onset, but not post movement onset.

Discussion
In the presented study, we aimed to evaluate the influence of artificial kinaesthetic vibrotactile feedback in a center-out movement task by comparing four conditions-real-time kinaesthetic feedback, static vibration, no vibrotactile input, and no movement with sham feedback.
In low-frequency EEG signals around the movement onset, we found minor differences in movement-related features, i.e. grand average MRCPs and SMR modulations, as a function of the condition as well as the movement direction.
Classification of movement conditions and directions against rest yielded moderate accuracies that were well above chance level in all cases. For the nonmovement condition, accuracies were around chance level in the case of the LF classifier, and above chance level but considerably lower than for movement trials in the case of the power classifier. For classification between conditions, the non-movement condition was most discriminable, with within-class accuracies  of the three movement conditions slightly above or around the significance threshold above chance level. Similarly, within-class accuracies for classification between directions narrowly exceeded the significance threshold.

Vibrotactile feedback
Participants reported the vibrotactile stimulation to be pleasant and intuitive.
In the frequency ranges considered in this work, no artifacts resulting from the presence of the stimulation were found in the EEG signals. We deem our vibrotactile stimulation setup suitable for providing real-time feedback of unidirectional movements, and are planning to expand its capabilities to more generalized 2D patterns in the near future.

Design of the trial structure
The paradigm was carefully designed with the intention of minimizing any contaminating influences on low-frequency signals during the main time window of interest: from the start of the MRCP that originates before the movement onset until the termination of the center-out movement. In a visually guided paradigm, visual cues present an additional source of contamination in the frequency range of interest, e.g. in the form of evoked potentials. Similarly, we expected an SEP when turning on the vibrotactile stimulation. Therefore, we made an effort (1) to separate the cue from the movement execution, and keep the movement period free from visual stimuli apart from the fixation cross, and (2) to introduce both the fixation cross and the stimulation well before the movement onset. The offset between the fixation cross onset and the earliest allowed movement onset was chosen following several pilot experiments.
Timings and topographic representations of the low-frequency signals suggest that evoked responses to visual and somatosensory stimuli were successfully separated from the potentials associated with the initiation of the movement.

Movement vs. non-movement
In all movement conditions, a potential can be identified as a prominent negativity above central areas, time-locked to the movement onset. The negative slope begins around 1.5 seconds before the movement onset, followed by a negative peak shortly after the onset, and a return to baseline within approximately 2 seconds, matching the general properties of MRCPs documented in literature [12][13][14]. In the non-movement condition, an MRCP is clearly absent.
Furthermore, in movement conditions, a characteristic power decrease can be observed above contralateral, and to a smaller degree ipsilateral, central/centroparietal areas in the µ and β frequency bands [7]. This power modulation is present seconds before the movement onset, and most intense between 1 second pre movement onset until the end of the analysed time window. In the non-movement condition, a parietal power increase is found in a frequency range overlapping with the µ band. It appears at the time of  For movement trials, grand average classification accuracies of LF classification against rest reach peak accuracies of 60%-65% shortly after the movement onset for all splits according to conditions and directions. The power classification yields peak accuracies of approximately 90%. Peak accuracies are reached within 500 ms post movement onset in both cases. For the LF classification, accuracies peak for a duration of 0.5 to 1 second before returning almost to chance level by 2 seconds post movement onset, while for the power classification, the peak accuracies are sustained. This reflects the properties of the respective features, since an MRCP is a transient potential indicating the onset of a movement, while the movementrelated power decrease of SMRs is sustained throughout the duration of the same.
For non-movement trials, accuracies obtained by LF classification do not significantly exceed chance level at any time window included in the classification. For power classification, accuracies for nonmovement trials against rest are around 75% for the whole time window. Evidently, the power features carry information beyond indicating a movement state, as evidenced by the accuracies for condition VF against rest, as well as the accuracies pre movement onset in the movement conditions.
As could be expected from inspecting the features, non-movement trials are better discriminated from the other conditions in four-way classification between conditions, for both LF and power classification, though the LF classification yields lower accuracies than the power classification. Evidently, the movement-related power features are more consistently represented on a single-trial and single-subject level.

Difference between movement conditions
Variations in the shape of an MRCP have been documented to reflect properties of the movement, such as the type of movement [17][18][19][20], the level of force or effort involved in producing the movement [15,16], movement speed [16], or whether the movement is targeted towards a goal [47].
We would expect the negative slope of the MRCP in conditions ME+VF and ME+Vstat to be highly similar, since the two conditions are identical before the movement onset, when subjects would realise whether they were receiving feedback or meaningless vibrotactile stimulation. In the grand average, the slopes as well as the spatial profile of the negativity indeed seem more similar between these two conditions than compared to condition ME. In both ME+VF and ME+Vstat, the LF classification accuracy against rest rises above the significance level a little earlier than in ME, while in ME the peak accuracy is sustained for a longer time.
After the movement onset, the potential of ME+Vstat seems to diverge from ME+VF, before the latter reaches its peak negativity, and after approximately one second to assimilate to ME. This would support the hypothesis that after the initial recognition, the static stimulation becomes background noise to the subject, since it does not carry any meaningful information.
However, we refrain from attaching definitive claims to these differences, as they are not backed up by an objective measure of statistical significance.
Differences in the relative µ and β power of the movement conditions are less evident in the grand average power, though they are picked up by the power classifier at a level significantly exceeding chance from shortly before the movement onset until the end of the classification time window.
Concerning behavioural measures, movements are initiated significantly earlier in trials with stimulation (conditions ME+VF, ME+Vstat) than in trials without (condition ME), while movement durations are not significantly different between conditions.
While features based on sensorimotor-rhythms have previously been shown to benefit from vibrotactile or electrotactile input [21,22] during motor imagery, we are not able to confirm a beneficial effect in this data set. We are considering the possibility that artificial feedback bears less importance to the subject in the case of movement execution with intact proprioception. We are not aware of comparable results with respect to MRCP amplitudes or low-frequency EEG features in general. As discussed above, our results hint at a potential positive effect on MRCP amplitudes, but we are not able to substantiate these hints. Moreover, we have not found any increases in performance linked to either the vibrotactile feedback or the presence of the vibrotactile stimulation. We attribute the high levels of variance to the variability introduced by the mixing of conditions and directions as well as the small number of subjects, and are planning to perform experiments with more narrowly tailored protocols that allow for the collection of a larger number of trials per condition, as well as a larger subject pool, in order to shed light on these issues.

Difference between movement directions
Examining the MRCPs by movement directions, regardless of conditions, we conclude that they have a similar shape on the initial phase of the negativity, but diverge somewhat shortly before and around the movement onset. For left movements, the divergence takes place earlier than for the other directions. The peak amplitude is more negative for left and up movements. Similarly, the µ power decrease is weaker for right movements, compared to the other two directions. These differences may be attributed to discrepancies in cognitive effort, or in force production [15,16] between movements to the three directions. While we attempted to make the movements to the three directions as similar as possible in terms of effort and comfort, by carefully arranging subjects' seating and hand starting position and defining the maximal range of motion accordingly, the reported behavioural measures in conjunction with meta data concerning aborted trials and feedback from subjects indicate that there are slight inherent differences in difficulty.
Classification of individual directions against rest yields similar results to the movement conditions against rest, with the accuracy of left vs. rest rising slightly slower before the movement onset.
The minor differences in MRCP amplitude are reflected in accuracies slightly over the significance level above chance for classification between movement directions. Recent results obtained by Kobler et al [48] in a center-out task suggest that the directional discriminability in low-frequency EEG is encoded more strongly and more consistently across individuals at the time point where the movement direction is revealed, rather than the movement onset. In the present experiment, since the presentation of the cue is separated from the movement, a separate analysis at the beginning of the trial would be necessary. Unfortunately, the necessity of analysing this time window was not anticipated, and participants were given greater leeway with respect to artifacts.
Further more focused investigations into the potential impact on directional discriminability are warranted, since the makeup of this data set is not suitable to separate directions individually for each condition.

Conclusion
The presence of vibrotactile input, both in the form of real-time kinaesthetic feedback and mere static vibration, influence the shape of the MRCP as well as SMR modulations to a small degree, and are expressed in behavioural differences, with respect to the timing of the movement onset. Movement trials are classifiable against rest with accuracies well above chance level. The accuracies are comparable for all movement conditions, indicating that in a control scenario where similar features are utilized, the proposed vibrotactile feedback would not impede control performance.
Furthermore, we found minor differences related to the movement directions, in grand average MRCPs, classification, and behavioural measures, but these need to be investigated further in a more narrowly tailored paradigm.