Disentangling the percepts of illusory movement and sensory stimulation during tendon vibration in the EEG

Mechanical vibration of muscle tendons in specific frequencies - termed functional proprioceptive stimulation (FPS) - has the ability to induce the illusion of a movement which is congruent with a lengthening of the vibrated tendon and muscle. The majority of previous reports of the brain correlates of this illusion are based on functional neuroimaging. Contrary to the electroencephalogram (EEG) however, such technologies are not suitable for bedside or ambulant use. While a handful of studies have shown EEG changes during FPS, it remains underinvestigated whether these changes were due to the perceived illusion or the perceived vibration. Here, we aimed at disentangling the neural correlates of the illusory movement from those produced by the vibration sensation by comparing the neural responses to two vibration types, one that did and one that did not elicit an illusion. We recruited 40 naïve participants, 20 for the EEG experiment and 20 for a supporting behavioral study, who received functional tendon co-vibration on the biceps and triceps tendon at their left elbow, pseudo-randomly switching between the illusion and non-illusion trials. Time-frequency decomposition uncovered a strong and lasting event-related desynchronization (ERD) in the mu and beta band in both conditions, suggesting a strong somatosensory response to the vibration. Additionally, the analysis of the evoked potentials revealed a significant difference between the two experimental conditions from 310 to 990ms post stimulus onset. Training classifiers on the frequency-based and voltage-based correlates of illusion perception yielded above chance accuracies for 17 and 13 out of the 20 subjects respectively. Our findings show that FPS-induced illusions produce EEG correlates that are distinct from a vibration-based control and which can be classified reliably in a large number of participants. These results encourage pursuing EEG-based detection of kinesthetic illusions as a tool for clinical use, e.g., to uncover aspects of cognitive perception in unresponsive patients.


Introduction
Long term immobilization due to continuing intensive care after traumatic or non-traumatic brain lesions for example, can lead to additional cortical reorganization, diverting resources away from the underutilized motor system, which can impede successful motor rehabilitation ( Khan et al., 2017 ;Langer et al., 2012 ;Opie et al., 2016 ). One promising tool to enhance neurorehabilitation in these cases is functional proprioceptive stimulation (FPS) ( Proske and Gandevia, 2012 ;Roll et al., 2012 ;Taylor et al., 2017 ). The principle behind FPS is tactile stimulation of muscle tendons, usually with vibration frequencies between 30 -100 Hz Calvin-Figuière et al., 1999 ;Gilhodes et al., 1986 ;Ribot-Ciscar and Roll, 1998 ). This effect is best observed in the absence of visual or tactile feedback about the stillness of the body part (see Taylor et al. (2017) for a comprehensive review).
Brain scans with positron emission tomography (PET) or functional magnetic resonance imaging (fMRI) have previously identified the neural correlates of the movement illusion perception across a network of brain areas, including: the primary somatosensory (S1) and motor cortex (M1) -mostly contralateral to the stimulated limb, and -bilaterally -the supplementary motor areas, premotor cortices, bilateral cingulate motor areas, and bilateral parietal cortices ( Naito et al., 1999 ;Naito and Ehrsson, 2001 ;Romaiguère et al., 2003 ). The intricate interplay of all these areas is not fully understood yet, but there is good evidence to support the hypothesis that the perceived illusory movement results from the interaction between sensory areas processing afferent proprioceptive signals, primary and higher order motor areas integrating proprioceptive signals in motor planning ( Amemiya and Naito, 2016 ;Naito et al., 2005Naito et al., , 1999 and parietal multisensory areas involved in representing the body in space ( de Vignemont et al., 2005 ;Romaiguère et al., 2003 ).
Few studies have examined the electrophysiological correlates of motor illusion by investigating its neural generators through source localization -either in the electroencephalogram (EEG) or the magnetoencephalogram (MEG) -providing consistent evidence of the involvement of the sensorimotor network in the perception of illusory movement ( Balconi et al., 2018 ;Casini et al., 2006 ). Further reports, exploring the neural underpinning of FPS through electrophysiology, have focused on the event-related potential (ERP) and the event-related (de)synchronization (ERD/ERS) in response to the vibration.
Little is so far known about the ERPs during FPS. Münte et al. (1996) described a distinct late negative deflection starting around 400 ms after FPS start, spread across fronto-central electrodes, whose power was modulated by stimulation frequency across the whole vibration period. Similarly, Saradjian et al. (2013) corroborated this late ERP negativity pertaining to the tendon vibration while testing its changes with respect to motor preparation. A handful of studies have investigated the ERD/ERS -a measure of spectral power change compared to baseline -as the neural correlates of kinesthetic illusion induced by FPS and have found them to share key characteristics with ERD/ERS elicited during motor execution and motor imagery ( Keinrath et al., 2006 ;Neuper et al., 2006 ;Pfurtscheller and Lopes da Silva, 1999 ;Yao et al., 2014 ). In this context, a prominent marker of illusory movement is represented by a mu-band desynchronization during the evoked activity, as measured by the event-related desynchronization (ERD) ( Barsotti et al., 2018 ;Imai et al., 2017 ;Keinrath et al., 2006 ;Yao et al., 2014 ). The reliable generation of this pattern was further confirmed by the accurate single trial classification of left vs. right arm FPS in Yao et al. (2014) and motor imagery in conjunction with FPS vs. rest in Barsotti et al. (2018) . However, as none of the aforementioned EEG studies had vibrotactile stimulation in a no-illusion control condition it remains unclear whether the observed signatures are specific to the illusion percept or if they are mostly driven by the stimulus sensation.
Tendon co-vibration, the stimulation of both tendons of an antagonist muscle pair (e.g., biceps and triceps), mimics the afferent muscle spindle inputs of the whole muscle group during movement ( Calvin-Figuière et al., 1999 ;Gilhodes et al., 1986 ;Ribot-Ciscar and Roll, 1998 ). Gilhodes et al. (1986) showed that vibration of two antagonist muscle tendons with different frequencies chosen between 20 and 100 Hz will result in an illusory movement that feels equivalent to an elongation of the muscle with the faster stimulation frequency (e.g., arm extension in case the triceps vibration is faster than the biceps one). Illusion intensity (onset, velocity, perceived angle) scales with the difference in vibration frequency between the two stimulators. Most interestingly, tendon covibration with matching frequencies leads to no -or only insignificantmovement perception ( Calvin-Figuière et al., 1999 ;Chancel et al., 2016 ;Gilhodes et al., 1986 ).
Our study is, to the best of our knowledge, the first attempt to identify the neural EEG correlates of illusory movements induced through FPS by contrasting two different conditions of tendon co-vibration, which, depending on their frequencies, did or did not induce an illusion of movement. With this experimental design we aim to minimize the difference in sensory processing between conditions to uncover neural patterns pertaining only to the kinesthetic illusion perception based on the ERP and ERD/ERS.
To this goal, we recruited 40 subjects in total, 20 for an EEG-based FPS study and another 20 for a corresponding behavioral study. The behavioral experiment aimed at quantifying the perceptions to the two employed vibration schemes. In the EEG cohort, we explored illusion correlates with a focus on features that have proven to contain information about executed or imagined movement ( Barsotti et al., 2018 ;Casini et al., 2006 ;Imai et al., 2017 ;Keinrath et al., 2006 ;Münte et al., 1996 ;Saradjian et al., 2013 ;Yao et al., 2014 ). Specifically, from the EEG data, we extracted spectral (ERD/ERS) and time-series (ERP) features to search for specific markers discriminating illusion-inducing versus nonillusion stimulations. We then evaluated the significance of these results at single-subject level by training and testing a classifier discriminating the two conditions on ERD/ERS and ERP features.

Participants
For the EEG experiment we enrolled 20 able-bodied people in the study. The group consisted of 12 females and 8 males, all right-handed and on average 24.6 ± 3.2 years old (mean and standard deviation indicated here and elsewhere). For the behavioral experiment, we enrolled another age and sex matched 20 subjects, 12 female and 8 male, all right handed and on average 24.5 ± 4.1 years old. All subjects participated voluntarily and consented in writing to the experiment. The study was covered by the ethical protocol No. 142/09 from the Commission cantonale d'éthique de la recherche sur l'être humain (CER -VD) and in agreement with the Declaration of Helsinki ( World Medical Association, 2013 ).

Experimental setup
The subjects sat comfortably in a chair facing towards their right side so to not see the stimulated left arm, which could have hampered the illusion of movement created during the tendon vibration ( Guerraz et al., 2012 ;Seizova-Cajic and Azzi, 2011 ). Even though the illusions were designed to mimic elbow extensions, we ensured that even with full elbow flexion the forearm did not enter the visual field to prevent possible intersensory conflicts. While their right arm rested comfortably in the lap, the left arm was supported by a movable forearm rest which allowed two degrees of freedom in the horizontal plane ( Fig. 1 ). The reason for this were two-fold. First, the whole weight of the arm was taken over by the support structure and thus the participants could completely relax their arm muscles during the experiment, which has been shown to be important for strong and vivid illusions ( Cordo et al., 2005 ;McCloskey, 1973 ;Taylor et al., 2017 ). Second, it is known that the proprioceptive feedback of the arm touching an immobile object can prevent the kinesthetic illusion from forming ( Rabin and Gordon, 2004 ) which was the reason for choosing a movable arm rest.
Tendon vibration was achieved with electromechanical wireless vibrators set into a soft, elastic brace on the left elbow joint (Vibramoov, Techno Concept, Manosque, France). The left arm was chosen since it was demonstrated that illusions start faster and are more vivid in the non-dominant extremity ( Tidoni et al., 2015 ). One vibrator was sitting against the distal biceps tendon and the other against the distal triceps tendon on the same arm. Time information about the beginning of each stimulation was sent via a cable link to the computer. and an elastic brace with two inserted vibrators on the left elbow, forearm support, recording computer and Vibramoov stimulation device. EEG components and respective communication are marked in red while tendon vibration related components are highlighted in blue. Communication between the EEG and the computer, as well as between the stimulation device and its vibrators on the arm was wireless. Experimental triggers were sent via cable link from the Vibramoov to the computer. The setup was the same minus the EEG cap for the behavioral experiment.
For the EEG experiment, subjects additionally wore an EEG cap with built-in wireless amplifier (g.tec Nautilus, g.tec medical engineering, Graz, Austria) with 16 electrodes covering the sensorimotor cortex in the international 10-10 system at positions (Fz, FC3, FC2, FCz, FC2, FC4, C3, C1, Cz, C2, C4, CP3, CP1, CPz, CP2, CP4). The signals were recorded at 500Hz with a hardware-implemented bandpass filter between 0.1 and 100 Hz and sent to a computer in the same room. The reference electrode was placed on the right earlobe.

EEG study protocol
EEG was recorded continuously while delivering stimulation sequences consisting of two different vibration types. The first was designed to elicit an illusion of elbow extension and was produced by vibrating the distal biceps tendon at 90Hz and the distal triceps tendon at 50Hz. According to Gilhodes et al. (1986) , this 40 Hz difference in vibration frequencies should lead to the feeling of a 20-25°extension in the elbow joint. The second stimulation condition was designed to produce only a vibration sensation without any accompanying movement illusion and consisted of stimulating both tendons at 70Hz ( Calvin-Figuière et al., 1999 ;Gilhodes et al., 1986 ). Thus, the average frequency of stimulation delivered to the agonist-antagonist pair was the same between conditions, but one condition was designed to induce a clear illusion and the other no illusion at all (control).
Each stimulation lasted three seconds and consisted of one second of linear frequency ramp-up, one second of a stable frequency interval and one second of linear frequency ramp-down. The linear ramps started and ended 10 Hz below the target frequency for each stimulation type. The amplitude of the vibration was 2-3 mm. These parameters were based on Romaiguère et al. (2003) and all settings were kept constant throughout the whole recording session.
Each subject underwent three blocks of 72 vibrations (36 illusion, 36 control), arranged pseudo-randomly and different for each block. The same stimulus sequences were employed for each participant. Inter stimulus intervals (end of vibration until start of next one) varied between one and three seconds and were randomized within and between blocks in order to minimize stimulus onset anticipation. We ensured that all subjects were feeling illusory movements using open questions after every block ( "What have you been experiencing during stimulations? ").

Behavioral study protocol
To ensure the design parameters for the EEG experiment created in the intended illusion versus no-illusion effect, we conducted an additional behavioral experiment using the same stimulation conditions and parameters. Participants were firstly accustomed to the vibration and resulting sensations for 4 min and then were asked to use this experience to initialize a linear grading scale from 0 to 10, where zero meant the vibration did not induce any illusion of movement and 10 denoted the largest deflection experienced. The subjects were unaware that there were only two different stimulation conditions. Participants underwent two different experimental blocks, both between 6 and 7 min long, in a pseudo-randomly defined order to eliminate time-dependent factors. In one block, subjects kept their eyes open, and their head turned -like in the EEG experiment -and experienced 48 stimulations (24 illusion and 24 control, pseudo-randomly arranged). The participants were asked to perform a mouse click with their right hand whenever they perceived the start of a movement and called out a number between 0 and 10 after each vibration ceased, corresponding to the intensity of the experienced movement. Compared to the EEG experiment, the inter-trial interval was elongated to 5-6 s to allow for the subject feedback.
In the other block, participants were blindfolded and had their right arm sitting on an elbow rest on top of a wooden board with radial markings of 5°spacing (see supplementary Fig. S 1). They received 32 stimulations (16 illusion and 24 control, pseudo-randomly arranged) and were asked to mirror the perceived left-arm movement with their right arm after every stimulation. Their head was kept straight ahead as to not induce skewed angles based on the head position and the experimenter reset their right arm to the zero-position on the board after every trial. Subjects were instructed to adjust the mirror image with respect to an eventual offset in their perceived left-arm position at the beginning of a trial. The inter-trial interval was stretched again to 10 s to allow for enough time for the mirror movement and subsequent reset of the right arm.

Behavioral data analysis
We extracted three metrics of interest from the data collected in the behavioral experiment: (i) the subjective rating of perceived illusion strength -z-scored for better comparability between subjects, (ii) the movement onset timing and (iii) the perceived movement angle. We averaged individual data points in these three metrics for each participant and stimulation condition and then pooled the data for all participants. We performed one-sample t-tests on the differences between illusion and control trials for each metric. In addition, we also performed t-tests to identify if the observed angles (metric iii) were significantly different from the starting position.

EEG data analysis
We examined time-and frequency-dependent differences in the EEG response to tendon vibration during illusion and control conditions. All data processing was carried out in MATLAB (The MathWorks, Inc., Natick, MA, USA) using the FieldTrip toolbox ( Oostenveld et al., 2011 ) and custom scripts.
We had to exclude block 1 of subject 13 due to a corrupted file. In addition to the hardware-implemented bandpass filter between 0.1 and 100Hz, we applied an offline lowpass filter below 40 Hz (Butterworth, 4th order) to the continuous EEG data. EEG epochs were extracted from -0.5 before to 4 s after the vibration onset. Epochs containing physiological artifacts, such as eye movements or muscle activity, were identified using a ± 100 μV threshold and additional visual inspection and were excluded from subsequent analysis (8.3 ± 9.3%). The number of excluded trials did not differ significantly between the two conditions (t(12) = -0.59, p = .56). Next, we applied baseline correction by subtracting the average voltage in the 500 ms pre-stimulation interval. Finally, and only for the time-frequency decomposition, we applied a surface Laplacian spatial filter using spherical splines for neighboring electrodes in a 5 cm radius (Perrin et al. 1989).
For the ERP analysis we kept the right-ear reference of the recording. This was motivated by the small head coverage and by a very low variability of the EEG response across channels. In such condition, rereferencing to a common average reference or surface Laplacian would have removed the signal of interest.

Mu and beta band range time-frequency analysis
Movement related brain correlates are known to be encoded in oscillations that can be examined using time-frequency analysis ( Baker, 2007 ;Kiebel et al., 2005 ;Lopes da Silva, 2013 ). To this end, we used a wavelet decomposition (Morlet wavelet) to compute the power spectral density (PSD) for each electrode separately in 1 s windows sliding across a period of -0.5 s before trial onset to 1 s after trial end with a 50 ms increment. We then calculated the event related desynchronization/synchronization (ERD/ERS) using the following formula ( Pfurtscheller and Lopes da Silva, 1999 ): where PSD t is the PSD at time point t, PSD bsl is the mean PSD along the baseline and T denotes the set of time points in the trial. A negative value resulting from (1) is called event-related desynchronization (ERD), since the PSD at time point t is lower than during baseline ( = less synchronous activity). If the value is positive, it is referred to as event-related synchronization (ERS).
To contrast illusion and control trials, we also calculated the differences in the ERD and ERS for each channel. The grand average timefrequency maps were obtained by taking the average across all 13 subjects, either for the computed contrast or the ERD/ERS values per experimental condition. For each channel, we tested the observed ERD/ERS (power changes with respect to baseline) across subjects using a cluster permutation test ( = .05), correcting for multiple time points and frequency bands ( Maris and Oostenveld, 2007 ). We also performed equivalent cluster permutation tests on the contrast signal between illusion and control condition ERD/ERS.

Event-Related Potentials (ERP)
In addition to the time-frequency analysis, we computed the eventrelated potentials (ERP).
After preprocessing (see section 2.6 ) we calculated the mean ERPs by averaging peristimulus epochs, time-locked to each vibration onset, for each channel and condition separately. Statistical analysis was based on the mean across channels (average ERPs) for each experimental condition. At the group level, we tested differences of the channel-averaged ERPs across subjects using a cluster permutation test ( = .05) to identify significant time intervals while correcting for multiple comparisons across time points ( Maris and Oostenveld, 2007 ).

Classification
We implemented a single-trial classification algorithm to evaluate how well the differences between the two conditions, as based on the time-frequency decomposition or ERPs, can be detected in individual subjects. For the time-frequency classification, the input features were the mean ERD/ERS values in 0.5-second-long time windows, covering the whole stimulation period from vibration onset to the end of the vibration (0 to 3 s) for each channel and three averaged frequency bands: mu/alpha (8-12 Hz), low beta (13-20 Hz) and high beta (21-30 Hz), resulting in 16 × 6 × 3 = 288 features. For the ERP classification we took the mean across channels, smoothed the ERP time series of each trial with a 100ms moving average and then downsampled the data to 10 Hz. The input features for the ERP classification analysis were all data points inside the 0 to 3 second time window (30 features).
For each subject individually, we implemented a 10-fold nested cross-validation on all available data, where the optimal features for classification were chosen in an inner 9-fold cross validation. To reduce the dimensionality of the ERD/ERS feature space, we kept only features that passed a two-sample t-test at significance level = .05. For both classifiers we then employed a forward feature selection algorithm using a linear discriminant analysis (LDA) classifier. Subsequently, the selected features were used to train and evaluate another LDA classifier in the 10 outer folds ( Varoquaux et al., 2017 ;Xanthopoulos et al., 2013 ).
Since trial-by-trial variances were quite high, and to ensure that results were not a product of the stratified but random selection of folds, we employed following procedures to achieve more stable estimates: (i) averaging responses to each condition in the test set (sub-averages) to reduce physiological noise and (ii) repeating the whole classification procedure 10 times and averaging the results. For each subject, accuracy is computed as the percentage of correctly classified sub-averages, averaged firstly across the 10 outer folds of each repetition and secondly across the 10 repetitions (see Fig. 2 for a schematic representation of the classification procedure). We also reported mean classification accuracy plus standard deviation across all subjects.
To test if the single-subject results are significantly above chance level, we repeated each classification procedure 100 times with shuffled labels, fit a normal distribution to the results and tested the previously computed classification performance against this random-permutation classification distribution ( Pereira et al., 2009 ). To this aim, we employed a one-sided Wilcoxon signed-rank test with a chance level of = .05.

Open data
The EEG data collected for this study is freely available as Open-Neuro dataset https://openneuro.org/datasets/ds003343 . Likewise, the analysis code can be found in a Git repository at https://github.com/ DNC-EEG-platform/MotorIllusion .

Behavior
To quantize the effect elicited by the chosen vibration setup, we performed a behavioral test in a separate cohort of participants matched for age and gender to the cohort of the EEG experiment. When participants rated the perceived illusion intensity on a scale from 0 to 10, the difference in pooled z-scored (to account for subjective bias) averages was highly significant (t(19) = -6.3, p = 5 × 10 -6 , Fig. 3 left panel). Mean reported movement onset time was 1183 ± 566ms for the illusion, and 1516 ± 628 ms for the control condition, and the difference proved highly significant (t(19) = 5.3, p = 4 × 10 -5 , Fig. 3 center panel). When subjects mirrored the perceived movement with their right arm, they indicated a mean elbow extension of 21 ± 12.9°in the illusion condition and 5.2 ± 8.7°in the control condition. These numbers corroborate previously reported results and confirm our choice in the stimulation parameter design ( Calvin-Figuière et al., 1999 ;Gilhodes et al., 1986 ). The difference in perceived angles was highly significant (t(19) = -7.89, p = 2.1 × 10 -7 , Fig. 3 right panel). The angles reported in the illusion condition were highly significantly different from zero (t(19) = 5.3, p = 6.4 × 10 -7 ) while the ones in the control condition were still significantly different, but much less so (t(19) = 2.68, p = .015). No subject reported any contortion of body perception (e.g., lengthening, twisting or detachment of the arm).

Fig. 2.
Classification and cross-validation procedure. Columns are referred to by their number (.) from left to right. The classification procedure (1) was repeated 10 times for each subject (100 times for the random permutation test). Each classification procedure was a nested cross-validation consisting of 10 outer folds (2) and 9 inner folds (3). One fold was kept as test set in each cross-validation, here exemplary the Outer fold 10 and the Inner fold 9. Each test set consisted of multiple trials of the illusion (4) and control (5) condition. To minimize the influence of physiological noise in the evaluation of the classification, we averaged all trials of each condition obtaining a mean illusion and a mean control trial (4-5 bottom, with black outline). The classifier was then tested on these two averages. This averaging procedure for the trials in the test set was performed in the inner and the outer cross validation.
To compute the classification accuracy for each subject, we averaged first the accuracies of the 10 outer folds (2) to receive the classification accuracy of each classification procedure (1) and then we averaged again over the 10 classification procedures in (1).

Fig. 3.
Behavioral differences in illusion perception. Left: Boxplots of the subject average intensity rating, individually z scored. Center: Boxplots of the subject average movement onset time (zero is vibration start, 3000 is vibration end). Right: Boxplots of the subject average perceived movement angle (positive is extension, negative is flexion). Bars with three stars indicate statistical significance below p < .0001.

Time-frequency analysis
Grand average plots of the time-frequency decomposition of the whole trial duration showed strong mu (8-12 Hz) and lesser beta (13-30 Hz) rhythm desynchronization in both conditions ( Fig. 4 bottom row), especially in the right hemisphere -contralateral to the stimulated left arm. The observed ERDs during vibration differed significantly from baseline activity -specifically over the contralateral hemisphere -as the result of the cluster permutation tests ( Maris and Oostenveld, 2007 ) correcting for multiple time points and frequency bands on the grand averages for each channel, done separately for the illusion and the control condition. The significant clusters are marked with black contours in Fig. 4 (bottom row). The frontal electrode Fz displayed significant timefrequency clusters only in the illusion condition. Across electrodes, this significant desynchronization started around 0.5 s after vibration onset and lasted on average until 0.5 s after vibration end.
When contrasting conditions (ERD/ERS illusion minus ERD/ERS control), we did not obtain any significant time-frequency clusters ( Fig. 4 top). Looking at the differences in the ERD and ERS values between the two experimental conditions, we qualitatively observed a mirrored image between the right and left centroparietal electrodes in both mu and beta band ( Fig. 4 top). While not significant, the contrast image showed an increased mu band ERD contralateral and a decreased mu band ERD ipsilateral to the stimulated side. Inversely but also not sig-nificant, we found an increase in the beta band ERD ipsilaterally while seeing it diminished contralaterally.
Averaging the ERD/ERS contrast between the illusion and control condition across the whole stimulation duration (0 -3s) shows lateralization patterns in both the mu and beta band ( Fig. 5 ). Nonetheless, cluster permutation tests across electrode locations and time did not find statistically significant changes of the ERD/ERS between vibration types across the scalp topography. Examining the topographical distributions in half-second windows across the trial duration did not yield significant clusters neither, nor did a split of the beta band into a low (13-20 Hz) and high (21-30 Hz) subband (supplementary Fig. S 2)

Event-related potentials
Both conditions elicited a similar, asymmetrical spatial distribution (more negative over the contralateral hemisphere) across the voltage montage along the stimulation period ( Fig. 6 ). The grand average ERPs across the trial duration show clear deflections following vibration onand offset ( Fig. 7 ). Since a cluster-permutation test across channels and time yielded no significance, we additionally collapsed all 16 channels to obtain a more robust signal ( Fig. 7 , bottom). We observed a negative peak 350 ms after vibration onset, followed by a component in the positive direction after 420 ms (see arrows in Fig. 7 , bottom). Furthermore, we observed a strong positive peak in both conditions 400  ms after vibration ended. The first negative peak was stronger across the whole frontal part of the montage whereas from the second, more positive peak on, a lateralization formed with more negative values in the hemisphere contralateral to the stimulated side. After the end of the stimulation (3.3 -3.6 s), during the positive peak, the spatial distribution looked inverse to the one observed during the first negative peak (0.295 -0.375 s; see Fig. 6 ).
At the group level, the illusion condition displayed a stronger negativity from 310 to 990 ms after stimulation onset than the control condition ( Fig. 7 , bottom row in gray), which proved significant in a cluster permutation test at the = .05 level ( Pernet et al., 2015 ). The effect size of the difference in means within this time window was d = 0.33.

Classification
We applied our single-subject classification algorithm on timefrequency and ERP features.
By binning the frequency into three bands and computing the average ERD/ERS values in 0.5 second windows across the stimulation period, we obtained 288 features in the form 16 channels × 6 time windows × 3 frequency bands. This approach yielded a mean classifica- Fig. 6. Grand average spatial voltage distributions look very similar between the illusion condition (top row) and control condition (bottom row) across different periods of interest throughout the vibration trial (see Fig. 7 ). The color schemes across time points are centered around the mean ERP voltage during each window and can differ between the two conditions, notably for the third and fourth window, where the scalp voltages of the illusion condition were more negative overall. All color bars denote a span of 4μV from blue to red (2μV below and 2μV above the indicated voltage), where red is more positive. The ERPs in the time window before stimulation onset (-0.5 -0 s) are very close to zero due to the baselining procedure ( Section 2.6 ). The color-coded shaded areas indicate the standard error of the mean across subjects. The bottom plot shows the mean across all channel grand averages and topographical schematics of the grand average voltage differences between the illusion and control condition in several time windows. The gray box indicates the period of significant negativity of the illusion over the control condition, as evaluated by a cluster permutation test ( = .05). The black arrows denote negative peaks in both the illusion and control condition ERPs at 350 and 750 ms; white arrows indicate positively oriented peaks at 420 and 3400 ms after vibration onset (the latter equals 400 ms after vibration cessation). Note the traditionally inverted y-axis (up is negative). tion accuracy across subjects of 53 ± 7 %. 13 out of 20 subjects showed above-chance ( = .05) classification performance (see Table 1 in the supplementary material).
Using 100 ms time windows across the entire vibration period as ERP features (30 features), we obtained a mean accuracy of 56 ± 6%. Comparisons with the chance level revealed that in 17 out of the 20 subjects the classifier performed significantly better than chance at the = .05 level (see Table 1 in the supplementary material). Single subject classification results are depicted side-by-side for both feature spaces in Fig. 8 .
When comparing the two classification approaches, we saw that the ERP classification yielded more above-random results than the ERD/ERS-based classification (17 subjects to 13). The two approaches did not significantly differ in their overall classification accuracy (t(12) = 1.69, p = .11). All participants but one could be classified with above chance level accuracy by using at least one of the two approaches, 11 out of 20 achieved significance with both features. An attempted classification analysis based on the combined feature space (ERPs + ERD/ERS) did not improve classification performance. Fig. 8. Classification results. Classification accuracies for each subject for ERD/ERS features (blue) and ERP features (red). The whiskers denote the standard error of the mean across the 10 repetitions of the classification. Stars indicate results which were significantly above chance level as indicated by the random permutation tests ( = .05).

Discussion
In this study we investigated EEG-based correlates of kinesthetic illusions induced by functional proprioceptive stimulation and their applicability in a classification framework. We compared an FPS condition inducing an illusion of movement with a control condition which provided on average the same amount of vibrotactile stimulation but only subtle movement illusion. At group level, we found statistically significant differences between the event-related potentials in response to the two vibration types. Further, we demonstrated that differences based on ERPs and ERD/ERS lead to successful single-subject classification analysis of the perceived illusion versus control conditions in all participants but one.

Event-related desynchronization/synchronization
In both experimental conditions we observed the beginning of an ERD in the first second after vibration onset, which stayed on a stable level until the end of the trial and was stronger contralaterally to the vibrated extremity than ipsilaterally ( Fig. 4 bottom row, Keinrath et al. 2006 ;Yao et al. 2014 ;Barsotti et al. 2018 ). A very probable explanation for this is a strong activation of the contralateral sensorimotor cortex due to the vibratory stimuli. For obtaining further insight in differences specific to the illusion perception, we proceeded to cancel out the common sensory activation by contrasting the patterns obtained during illusion and control trials ( Fig. 4 top).
On the assumption that any canceled-out activity was purely sensory evoked, the remaining mu band ERD in the electrodes contralateral to the stimulated arm ( Fig. 5 ), together with an ipsilateral ERS, could suggest a stronger activation of the contralateral motor cortex with simultaneous inhibition of the ipsilateral motor cortex during the induced illusion ( Nam et al., 2011 ;Pfurtscheller and Lopes da Silva, 1999 ). This interpretation is consistent with the fMRI results by Romaiguère et al. (2003) , who, using a comparable experimental protocol, observed strong activations in S1 for all conditions -elicited by the tactile stimulation -while the contralateral M1 was only active during illusion conditions (see Casini et al. (2006) for similar findings using MEG). We also observed a widespread and temporally stable mu desynchronization from frontal to contralateral parietal electrodes ( Fig. 5 , supplementary Fig. S 2). Based on the theory that cortical mu/alpha rhythms are functionally inhibiting neuronal populations ( Klimesch et al., 2007 ;Jensen and Mazaheri, 2010 ), lower mu rhythm band power across frontal, central/contralateral and parietal areas could -under all caution necessary due to the non-significant results -hint towards a stronger fronto-parietal network activation during illusion trials. This would be consistent with previous reports that stated that the emergence of kinesthetic illusions is dependent on a fronto-parietal network, including the supplementary motor area, the premotor cortex and the cingulate motor cortex as well as the inferior parietal lobule ( Kenzie et al., 2018 ;Naito et al., 1999 ;Naito and Ehrsson, 2001 ;Romaiguère et al., 2003 ). We might speculate that lower mu-related inhibition in these regions reflects a high-level motor component associated to the movement illusion, as well as the involvement of multisensory parietal areas coding for the location of body parts in spacesee e.g. Bolognini and Maravita (2007) ; Azañón et al. (2015) . Continuing the interpretation of the observed beta band changes in the light of motor-related activity, the qualitative beta band ERS in the contralateral hemisphere is at first glance counterintuitive to the well-established connection of beta ERD during movement and motor imagery ( Neuper et al., 2006 ;Pfurtscheller and Lopes da Silva, 1999 ). However, an interesting approach is to understand the participant instruction for passiveness, i.e., to keep the arm relaxed throughout the vibration stimulus, as an active process during illusion, since subjects often reported an urge to move their limb in unison with their perception. This urge could be rooted in the activation of the SMA by the illusion ( Fried et al., 1991 ). Subsequently, one could speculate that the beta band ERS reflects the active attempt of the contralateral motor cortex to maintain the nonmoving status quo ( Engel and Fries, 2010 ).
Alternatively, the relative beta band power increase in the illusion versus the control condition might be caused by a heightened attentional status for sensory integration during the illusion perception ( Baker, 2007 ). Similar contrast patterns have been reported also for auditory illusions ( Leske et al., 2014 ;Theves et al., 2020 ), adding evidence for a non-motor related contribution to this pattern. Ultimately, since for the illusion-control contrast neither frequency-time pairs nor spatial locations showed significant differences, further investigations are needed to support or dispel such hypotheses.

Tendon vibration evokes many of the same EEG patterns even without illusion
Our results showed that tendon vibration, even without inducing illusory movement, induces neural patterns in the EEG that share key characteristics with the patterns present during the illusion condition. This includes the ERP-peaks after vibration on-and offset, as well as the long-lasting induced ERD over the contralateral motor cortex -EEG indices that have been previously reported for FPS, motor imagery and active or passive arm movements ( Barsotti et al., 2018 ;Imai et al., 2017 ;Keinrath et al., 2006 ;Yao et al., 2014 ). This suggests that sensory correlates like wide-spread mu-band suppression in the primary sensory cortex ( Kuhlman, 1978 ;Nikouline et al., 2000 ) during non-functional tendon vibration, i.e. without induced illusion, are likely creating a substantial proportion of the EEG correlates reported by previous EEG studies on FPS ( Keinrath et al., 2006 ;Yao et al., 2014 ). Thus, we conclude that it is indispensable to rule out the effects of vibrotactile stimulation to identify specific electrophysiological markers of illusory movement induced by FPS.

Event-related potentials
Our finding of a first negative peak around 350 ms and a subsequent late negativity starting from around 420 ms ( Fig. 7 ) corroborated earlier observations of EEG potentials evoked by FPS ( Münte et al., 1996 ;Saradjian et al., 2013 ). Münte et al. (1996) described a very similar, widespread fronto-central negativity forming around 400 ms after vibration onset, which lasted up to 250 ms after the end of their onesecond long trials. Saradjan et at. (2013) reported a first negative peak at 240 ms and the beginning of a large, long-lasting negativity around 400 ms after vibration onset. Also, the spatial location of our late ERP negativity, frontally towards the contralateral side of the stimulated arm (see Fig. 6 ), is in agreement with the results from Münte et al. (1996) .
As a novel result, we found the ERPs in the illusion condition significantly more negative from 310 ms until 990 ms after stimulation onset, compared with a tendon co-vibration control. It can be expected that somatosensory evoked potentials, which are correlates of early sensory processing (first 100 ms of the ERP), will not differ between conditions ( Cruccu et al., 2008 ;Desmedt et al., 1983 ). In contrast, the onset of the discovered significant negativity lies 310 ms after vibration start, and therefore well in the time window of other perception-related ERP components ( Chen et al., 2010 ;Galvez-Pol et al., 2020 ;Guillem et al., 1999 ). We speculate therefore, that the later ERP components, which display a significant difference between conditions, probably represent perceptual components relating to the illusory movement or movement onset. ( Codispoti et al., 2007 ;Münte et al., 1996 ;Parasuraman et al., 1982 ). This is supported by the data collected in the behavioral experiment (illusion onset reporting time), where -when we compensate for a typical reaction time of 200 to 300 ms in a speeded response task ( Jain et al., 2015 ) -the end of the differential ERP activity matches well with the decision time in our behavioral experiment. One could hypothesize that the larger ERP negativity indicates a perceptual processing and evidence accumulation period. Importantly, by contrasting the EEG response between the illusion and non-illusion (control) condition -which both included tendon vibration on the same arm -we can conclude that the magnitude of the first part of the observed late negativity (up to one second after stimulation onset) is not simply a product of the vibrotactile stimulation, but rather likely influenced by the illusion perception itself. This extends previously published FPS-EEG studies which did not include equivalent control conditions.

Classification
To the best of our knowledge, we were the first to classify EEG patterns during illusory movement against a control condition which involved the same amount of vibrotactile stimulation. Previous studies indeed classified illusion versus rest ( Barsotti et al., 2018 ) or left vs. right side movement illusion ( Yao et al., 2014 ), and actually reported higher classification accuracy. However, here, by contrasting two active stimulation conditions over the same arm region, we were able to classify features associated to the subjective component of the illusions in all but one participant. Further, we report a large inter-subject variability in the feature space, reflected in the span of the achieved accuracies from 44 to 74 % -similar to other EEG-based classification studies ( De Vos et al., 2014 ;Pfurtscheller et al., 2006 ;Pfurtscheller and Solis-Escalante, 2009 ;Yao et al., 2014 ).
Interestingly, the ERP features have yielded better overall classification and a smaller variance between different subjects than the ones obtained with the ERD -at least with the simple LDA classifier used. One possible reason for this might be the much smaller feature space in the ERP case (30 time points). Given the only around 200 vibration trials, the 288 ERD/ERS features facilitate overfitting when training the classifier -even with the feature space dimensionality reduction techniques employed -which consequently deteriorates test performance ( Mwangi et al., 2014 ;Verleysen and François, 2005 ). We suspect that substantially larger datasets would increase overall classification performance. However, this is in direct conflict with clinical feasibility -a cornerstone of our study design.

Limitations
The main limitation of our study is probably that illusion intensity and onset time were not collected from the participants in the EEG experiment. This is not so much a concern regarding the question if the subjects experienced kinesthetic illusions, since reliable illusion generation was confirmed in a separate behavioral experiment, but rather when the illusion started in each trial, and how strongly they were perceived. This information would have allowed us to dive much deeper into the analysis (time locking ERPs on illusion onset, correlating perceived intensity with ERD/ERS and ERP features) which might have helped in the interpretation of the results. However, the common approach of telling the participants to react at illusion onset, e.g. via a button press or by speaking out, would invoke additional motor-planning and subsequent execution patterns, like e.g. the readiness potential, which would probably obfuscate the illusion-related signals ( Caldara et al., 2004 ;Cunnington et al., 2002 ;Desmurget and Sirigu, 2009 ;Lee et al., 1999 ). The introduction of a Libet's-like paradigm ( Libet, 1985 ) -i.e., asking the participants to lately report the perceived time of the illusion onset with respect to an external event (e.g. a clock) -might be a viable option, if ensured that reporting will not elicit brain activations which overlap with the illusion-related ones. Granted, while diving deeper into this matter would certainly hold scientific value, we tried to design our study in a way that it can be extended to clinical use -specifically to people with disorders of consciousness ( Giacino et al., 2014 ) which require a fully passive protocol. In the end, not informing our analysis with the illusion start time or intensity underlines the fact that the observed phenomena are strong enough to be detected despite the uncertainty regarding the strength and onset.
Furthermore, one could argue that the differences in the ERPs between experimental conditions stem from the stimulation frequency difference between the two vibrations and therefore does not really represent the illusion perception. While we cannot entirely exclude this possibility, we have tried to mitigate this effect by considering the contrast between two conditions where the physical features (vibration depth, duration and average frequency) of the stimulation were matched on average. If the physical features of the stimulation were the main driver of the differential effect between the illusion and control conditions, we would have expected an earlier latency and longer lasting effect along the whole stimulation period. The observed effect, starting at 310 ms post-stimulus onset and disappearing around the one second mark, therefore supports the hypothesis that the ERPs are indeed modulated by the perceived illusion rather than a simple sensory effect.
The 16 electrode EEG montage provided a rather low spatial resolution for EEG signals and made source localization by inverse solution practically impossible ( Lantz et al., 2003 ;Michel et al., 2004 ). Nonetheless, the montage was covering the motor and premotor areas, which were the primary target regions in our study. Further, having demonstrated the classification of induced illusions with such a minimal setup, we open the door for applications that do not allow the time nor precision needed for high-density EEG setup, e.g., clinical and/or bedside use.
We did encounter not only high inter-subject variability, but also a high trial-to-trial variability for certain subjects. This compelled us to implement sub-averaging over the test set in the classification algorithm to raise the signal-to-noise ratio (see Fig. 2 ). Therefore, our measure of classification accuracy does not correspond to what is commonly referred to as single-trial classification accuracy, and applications would need to average across multiple trials in a real-time setting, as is standard in several brain-computer interface systems ( De Vos et al., 2014 ;Schalk and Mellinger, 2010 ). Our classification algorithms are thus applicable in all settings as long as fast decoding speed is not of primary concern.

Applications and further research
Many existing or envisioned applications of FPS revolve around motor learning, pain relief or rehabilitation after injury or illness ( Avanzino et al., 2014 ;Barsotti et al., 2018 ;Imai et al., 2017 ;Müller et al., 2002 ;Roll et al., 2012 ;Yao et al., 2014 ). Additionally, some studies focused on the consequences of spinal cord and brain lesions on the different stages of the sensory-motor network involved in generating illusory movements during FPS ( Desmurget and Sirigu, 2009 ;Müller et al., 2002 ). The ultimate goal of our study was instead to develop EEG-based markers of this illusion-generating network. A potential application would be to search for these biomarkers in patients with disorder of consciousness, where a critical clinical need is identifying patients with (at least partially) preserved conscious states in the absence of their ability to report them due to severe motor deficits ( Edlow et al., 2017 ;Pincherle et al., 2019 ). In this view, neural correlates of FPS could represent an implicit electrophysiological index of partially preserved consciousness (see Noel et al. (2020) for a similar approach). Finally, while fMRI has been the de-facto standard for examining cerebral areas involved in kinesthetic illusions, when it comes to regular use in diagnostic or therapeutic settings, it has large shortcomings in cost, availability and the possibility for bedside or out-of-clinic application. Our work showed that EEG has the potential to detect correlates of movement illusions with a lightweight and portable setup.

Conclusion
In this study, we demonstrated that the EEG contains significant and classifiable correlates of FPS-induced kinesthetic illusions as compared to a tendon co-vibration control condition. We validated the hypothesized strong perceptual difference between the conditions via a behavioral experiment and we could further corroborate findings of spectral power and ERP changes during vibration-induced illusions. Overall, our study extends current knowledge about the neural correlates of illusion processing in healthy participants in terms of latency and EEG topographic responses. In addition, we found that ERP correlates were the better features for classification and yielded levels of accuracy high enough to apply this setup in practical settings. This makes a strong case for pursuing EEG-based illusion detection as tool for prognostic or therapeutic use.

Declaration of Competing Interest
This study was part of a cooperation between the Acute Neurorehabilitation Unit (LRNA) of the Centre Hospitalier Universitaire Vaudois (CHUV) and Techno Concept. Techno Concept financed the first author during his employment at the LNRA. Furthermore, the Fondation CHUV and Loterie Romande provided grants for this work. Nevertheless, all involved authors were completely free in their choices concerning the experiment, analysis, writing and subsequent publication.
Data and Code Availability Statement: As also included in the materials section in the manuscript, we made all our data and code publicly available. The data collected for this study is freely available as OpenNeuro dataset https://openneuro.org/datasets/ds003343. Likewise, the analysis code can be found in a Git repository at https://github.com/DNC-EEG-platform/MotorIllusion.