Using machine learning to reveal the population vector from EEG signals

Objective. Since the discovery of the population vector that directly relates neural spiking activity with arm movement direction, it has become feasible to control robotic arms and neuroprostheses using invasively recorded brain signals. For non-invasive approaches, a direct relation between human brain signals and arm movement direction is yet to be established. Approach. Here, we investigated electroencephalographic (EEG) signals in temporal and spectral domains in a continuous, circular arm movement task. Using machine learning methods that respect the linear mixture of brain activity within EEG signals, we show that directional information is represented in the temporal domain in amplitude modulations of the same frequency as the arm movement, and in the spectral domain in power modulations of the 20–24 Hz frequency band. Main results. In the temporal domain, the directional information was mainly expressed in primary sensorimotor cortex (SM1) and posterior parietal cortex (PPC) contralateral to the moving arm, while in the spectral domain SM1 and PPC of both hemispheres predicted arm movement direction. The different cortical representations suggest distinct neural representations in both domains. Significance. This direct relation between neural activity and arm movement direction in both domains demonstrates the potential of machine learning to reveal neuroscientific insights about the dynamics of human arm movements.


Introduction
The ability to generate volitional movements is essential for independent everyday interactions. The control of an upper limb neuroprosthesis or a robotic arm using brain activity in people with spinal cord injury has enticed many researchers for more than four decades (Georgopoulos et al 1982, Taylor et al 2002, Wolpaw et al 2002, del R. Millan et al 2004, Rupp et al 2015, Müller-Putz et al 2017. Such control has been achieved by decoding kinematic signals recorded invasively from brain areas involved in motor control (Hochberg et al 2006, Collinger et al 2013, Ajiboye et al 2017. The control of an upper limb neuroprosthesis or robotic arm has also been achieved noninvasively, for example through detecting distinct 4 Author to whom any correspondence should be addressed. mental imaginations within electroencephalographic (EEG) brain activity (Pfurtscheller et al 2000, Müller-Putz et al 2005, Rohm et al 2013, Meng et al 2016. However, it has so far remained unclear how arm directional information is represented in EEG and how this representation relates to the behavioral response that describes the kinematics of the movement.
In the 80's, Georgopoulos et al (Georgopoulos et al 1982, 1983, Georgopoulos 1988) measured the activity of single neurons in rhesus monkeys during center-out reaching tasks. From the summation of single-neuron activity, they found a vector that pointed into the direction of the arm movement, and named it the population vector. Ten years later, Caminiti et al (1990) showed also in non-human primates that the population vector would still point in the direction of the movement when performing discrete movements in similar directions but using different sets of muscles. Researchers have attempted to decode directional information also from the human brain signals recorded invasively during center-out tasks (Mehring et al 2003, Ball et al 2009) and during continuous (Schalk et al 2007, Pistohl et al 2008 or attempted movements (Truccolo et al 2008, Chadwick et al 2011. Impressive results showing the control of robotic arms of up to ten degrees of freedom, as well as upper-limb neuroprostheses using invasively acquired directional information from the human brain have been published in recent years (Collinger et al 2013, Wodlinger et al 2014, Ajiboye et al 2017. In parallel to the developments in the invasive field, also non-invasive EEG movement covariates have been used either to control neuroprostheses which enable hand grasping in individuals with spinal cord injury (SCI) (Pfurtscheller et al 2000, Vučković and Sepulveda 2012, Rohm et al 2013, Rupp et al 2015, or to control the movements of a robotic arm (Baxter et al 2013, Meng et al 2016). Studies applying traditional non-invasive methods to decode covariates of arm movement typically target whole arm movements (Hammon et al 2008, Waldert et al 2008, Bradberry et al 2010, Ofner and Müller-Putz 2012, Yeom et al 2013. Moreover, the EEG frequencies that encode such movement covariates are under ongoing debate. Most studies focused on low-frequency (below 6 Hz) linear decoding (Bradberry et al 2010, Ofner and Müller-Putz 2012, Yeom et al 2013, Kobler et al 2018. Others argue that the mu and the beta band do not hold information about the trajectory of the arm (Waldert et al 2008, Ball et al 2009, Nakanishi et al 2013, Yeom et al 2013. However, it is possible that some of these results were misinterpreted due to the properties of linear regression methods (Antelis et al 2013). Korik et al (2018) state that the spatiotemporal power pattern of various frequency bands (including mu and beta) hold significant information about movement parameters. Seeber et al (2016) also showed that there are modulations in the mu and beta frequency bands when opening and closing the hand rhythmically. Furthermore, Lv et al (Lv et al 2010) state that both low-frequency (0.1-4 Hz), mu and beta frequency bands (12-30 Hz) encode information about arm movement velocities.
Although the concept of Georgopoulos' population vector cannot be applied directly to non-invasive recordings, we hypothesize that a population vector can be inferred non-invasively from trajectoryrelated directional information encoded in the EEG signals. Using machine learning methods that comply with the linear mixture of sources model (Nunez and Srinivasan 2006), we identified brain sources that co-modulate with arm movement direction. More specifically, we found brain sources whose activity in the low-frequency time domain (~0.5 Hz) and brain sources whose activity in the beta band-power domain (20-24Hz) resemble the movement direction of a rhythmic circular movement of the right arm.

Experimental design 2.1.1. Participants and experimental setup
Twelve healthy right-handed participants, aged 20-29 years (average age of 24 years, eight male participants) with normal or corrected to normal vision, participated in this study. The study conformed to the Declaration of Helsinki and was approved by the local ethics committee (approval number 29-058 ex 16/17). After the participants were informed about the purpose and the procedure, they gave their written informed consent to participate in the study. Figure 1(a) depicts the experimental setup. The participants sat in a comfortable chair in an electrically shielded room. A Kinect sensor (Microsoft Inc., USA) positioned at a distance of 1.2 m relative to the participants, was used to capture right arm movements. We mounted a black cardboard above the right shoulder to prevent visual feedback about the right arm movements. For the first participant, we did not use any occlusion apparatus; the next two participants were wearing cardboard goggles to occlude visual feedback. The experimental paradigm was displayed on a monitor located 1.5 m in front of the participants.

Experimental paradigm
The experiment consisted of twelve runs with a break after each run (figure 1(b)). A run contained seven trials. Figure 1(c) depicts the structure of a trial. Each trial started with the appearance of a blue fixation cross, which marked the beginning of a preparation period. The participants were advised to focus their gaze on it to minimize eye movements during the trials. The preparation period lasted between 3 and 4 s (pseudo-randomly distributed). Next, a change in fixation cross color indicated the task. The paradigm separated two tasks, namely movement and rest, with a respective ratio of 4:1. In the movement task (green fixation cross), participants were asked to perform a continuous circular right arm movement, while in the rest task (red fixation cross), they were asked to stay in a resting position, with their body relaxed and both hands on their lap. The disappearance of the fixation cross marked a break between trials. During the 8 to 10 s break (pseudo-randomly distributed), the participants were instructed to return their right arm to the lap and relax.
Regarding the movement task, we instructed the participants to perform a continuous, circular, clockwise movement on the vertical plane as if they would draw a circle to the right of their body (figure 1(a)). We asked the participants to perform the circular movement with a diameter of the size of their torso, with uniform speed and with a duration of 2 s per (a) Illustration of the experimental task and setup from three different perspectives. A black cardboard occluded movements of the right arm. Circular arm movements were performed on a 2D plane with a diameter of the size of the torso. A Kinect sensor (Microsoft Inc., USA) was used to capture right arm movements. (b) Outline of the experimental procedure. The experiment consisted of 12 runs of 7 trials each. (c) Experimental paradigm. One trial would start with a blue fixation cross in the center of the screen (preparation). Next, the cross would either turn green (movement task) or red (rest task) for 20 s, followed by a break of 8 to 10 s. cycle; i.e. at rate of 0.5 Hz. The experiment started after the participants practiced the movement and felt comfortable performing the task. We recorded ten rather than twelve runs for the first six participants.

Data acquisition
EEG signals were recorded with 61 active Ag/AgCl electrodes arranged according to the international 10-10 system (actiCAP and BrainAmp, Brain Products GmbH, Germany) and sampled at 1 kHz (anti-aliasing filter cuf-off frequency 450 Hz). Reference and ground electrodes were placed at the right mastoid and AFz, respectively. EOG activity was recorded with three additional active Ag/AgCl electrodes that were placed above the nasion and below the outer canthi of the eyes ( figure 1(a)). Electrode to scalp interface impedances were kept below 20 kΩ. The electrode positions on the scalp were acquired using an ultrasound based 3D positioning system (ELPOS, Zebris Medical GmbH, Germany). The right arm movements were recorded with a motion capture system (Kinect, Microsoft Inc., USA) at a rate of 30 Hz. It captured the positions of body joints (hand, wrist, elbow, shoulder, …) during the experiment. We used the lab streaming layer (LSL) protocol (https://github.com/sccn/labstreaminglayer) to record and store all signals.

Data pre-processing
The data processing was implemented in MAT-LAB (version R2015b, Mathworks Inc., USA) and EEGLAB toolbox (version 14.0.0b) (Delorme and Makeig 2004). Before we pre-processed the motion capture and EEG signals, we synchronized them offline.

Motion capture signals pre-processing
The 3D motion capture signals of the right hand position were band-pass filtered (infinite impulse response (IIR) Butterworth filter; fourth-order; zero phase; cut-off frequencies 0.3 and 0.8 Hz). We then used the Euclidean norm of the filtered 3D hand position signals to detect trials during which the sensor lost track of the hand or if a participant's circular movements were inconsistent compared to the other trials. In detail, we computed the mean of the Euclidean norm for every trial. Assuming a normal distribution across trials, we estimated the distribution parameters and detected an improbable trial if its mean was outside 3 standard-deviations around the mean across trials. We also tested the variance of the Euclidean norm in a similar fashion. On average, 3.2 trials were marked for rejection. Supplementary table 4 lists the detailed results per participant.
Next, we applied principal component analysis (PCA) to estimate the 2D plane on which a participant mainly performed the circular movement. We extracted the first two principal components (PCs). To maintain consistent orientation of the new 2D coordinate system across participants, we rotated the PC space so that the axes were maximally aligned to the horizontal and vertical axes of the original 3D space. As a result, the first axis was maximally pointing in the RIGHT direction and the second axis maximally in the UP direction (Supplementary figure 1 which is available online at (stacks.iop.org/JNE/17/026002/mmedia)). Consequentially, we labeled the projected signals as RIGHT and UP components. All subsequent steps were performed in the RIGHT and UP component space.
We detected peaks in each of the four directions, namely UP (local maximum in the UP component), RIGHT (local maximum in the RIGHT component), DOWN (local minimum in the UP component), LEFT (local minimum in the RIGHT component). For each direction, we rejected improbable peaks with a Mahalanobis distance larger than six to the distribution center. We also rejected peaks within a trial which were too close or too distant in time; i.e. two consecutive peaks appeared earlier than half of the nominal duration (= 0.5 s for a quarter of a cycle) or later than double the nominal duration. We only accepted cycles which were complete, i.e. a sequence of UP, RIGHT, DOWN, LEFT and UP peaks. Samples within complete cycles were labelled according to the quarter in which they were located.

EEG pre-processing
The EEG signals were first resampled at 200 Hz (Matlab function resample, finite impulse response anti-aliasing filter, 100 th order, zero phase, cut-off frequency 96 Hz). We then applied a high-pass filter (IIR Butterworth filter; fourth-order; zero phase; cutoff frequency 0.3 Hz) and a notch filter (IIR butterworth; second-order; zero phase; cut-off frequencies 49 and 51 Hz) to remove electrode drifts and line noise. We visually inspected the recorded signals for bad channels. A channel was considered as bad, if it was regularly contaminated by (1) pops and drifts, or (2) muscular (high-variance and high-frequency) activity. On average, 2.3 channels were marked as bad (supplementary table 4). The bad channels were spherically interpolated using the eeg_interp function (EEGLAB toolbox). Thereafter, the data were rereferenced to the common average of all channels. We then applied robust principal component analysis (Robust PCA) (Candès et al 2011) to attenuate occasional electrode pops and low-frequency drifts (Kobler et al 2019).
To identify outlier trials, we extracted 20 s epochs of EEG data, starting 1 s before the task cue in each trial. A trial was marked for rejection if the EEG signal of any channel exceeded a threshold of ± 250 µV or had an abnormal probability or kurtosis (more than five standard deviations beyond the mean). We applied the joint probability and kurtosis rejection criteria twice to detect large outliers in the first iteration, and subtle outliers in the second iteration.
We applied independent component analysis (ICA) to the epoched EEG data to identify artifactual sources. In detail, we used the extended infomax algorithm to extract independent components (ICs) that explained 99.9% of the variance. We then used SASICA (Chaumon et al 2015) to classify ICs into artifactual and non-artifactual. We adjusted the classification result by visually inspecting each IC's spectrum, topographic distribution and activation signal. On average, 14.3 of 53.7 ICs were marked as artifactual (supplementary table 4). The artifactual ICs were removed from the continuous EEG signals. Finally, the cleaned EEG channels were rereferenced to their common average. The entire pre-processing pipeline is outlined in Supplementary figure 8.

Low-frequency time domain analysis
We applied a low-pass filter (IIR Butterworth filter; second-order; zero phase; cut-off frequency 0.8 Hz) to the continuous EEG data and downsampled the dataset to 20 Hz before we extracted the 20 s epochs. Then, the previously marked outlier trials were rejected. We applied time-warping to the quarters of each complete cycle so that their duration after time-warping was 0.5 s. As a last data pre-processing step, we computed the power of the EEG channels for each quarter and rejected a quarter if the power of any channel was outside an interval of five standard-deviations around the mean. The extracted data was then used to decode directions.

Decoding directions from low-frequency EEG
As outlined in figure 2(c) (blue inset), we applied partial least squares (PLS) regression (Wold et al 2001) to decode directions from the pre-processed EEG x[t] (Ofner et al 2015, Kobler et al 2018 We assessed within participant significance by randomly rotating the 2D hand position signal of each trial. This approach maintained the autocorrelation structure of the direction signals within trials but  (c)), projected to the cortical surface (EEG source imaging). In source space, we averaged the norm of each voxel over participants. (g) as in (f) for the patterns associated to the SPoC models.
broke the structure across trials. We repeated the random rotation approach 1000 times per participant to estimate the chance level distribution of correlation coefficients. Its 95 percentile was used as significance level.
We also computed decoder patterns a PLS associated to the extracted weights  and scaled them with the decoded signal's standard-deviation to convert the patterns to µV at the scalp level (figure 2(c), blue inset).

EEG source imaging
We applied EEG source imaging (Michel et al 2004) to summarize the time domain results at group level. Participant specific head models were generated in the open source software Brainstorm (Tadel et al 2011) by co-registering the ICBM152 (Fonov et al 2011) template boundary element head model (BEM) with the participant specific electrode positions. The BEM consisted of three layers (cortex, skull, scalp) with relative conductivities (1, 0.008, 1). The cortex was modelled with 5011 voxels. BEM and electrode positions were co-registered by three anatomical landmarks (nasion, left and right preauricular points). Due to deviations between participant and template anatomy, we completed the co-registration by projecting floating electrodes to the scalp. Open-MEEG (Kybic et al 2005, Gramfort et al 2010 was used to compute the forward model; that is, to describe the propagation of the electric fields from cortex to scalp. sLORETA (Pascual-Marqui 2002) was applied to compute the corresponding inverse model for unconstrained sources (three components per voxel). We used the preprocessed EEG during rest trials to estimate the noise covariance between the EEG channels. The resulting noise covariance matrix was regularized by adding an identity matrix which was scaled by 10% of the noise covariance matrix's largest eigenvalue.
Before we mapped a participant's decoder patterns a PLS to source space, we normalized them with the global field power (GFP) to reduce participantdependent scaling effects (Kobler et al 2018). In detail, the source space pattern a (s) associated to the pattern a was with K being a n voxels · n voxelcomponents × n channels projection kernel (estimated with sLORETA) and g being the participant specific GFP. The scalar g was set to the median of the GFPs associated to the 1000 randomization models. For each randomization model r, the GFP g (r) was the standard-deviation across channels of the pattern a (r) .

Source space statistics
Group level significance of the decoder patterns was assessed in source space for the norm of each voxel; i.e. the Euclidean norm of its 3 components. Each voxel's pattern activity norm was compared to a baseline with a two-tailed, paired, nonparametric, permutation t-test (Nichols and Holmes 2002) (1000 permutations). We estimated the baseline through averaging the activity norms of the 1000 randomization models. Regarding multiple comparisons, we controlled the false discovery rate (FDR) at a significance level of 0.05 (Benjamini and Hochberg 1995).

Beta band-power domain analysis
We applied a band-pass filter (IIR Butterworth filter; sixth-order; zero phase; cut-off frequencies 20 and 24 Hz) to the continuous EEG data before we extracted the 20 s epochs. Then, the previously marked outlier trials were rejected. The remaining data was partitioned into overlapping segments. Each segment corresponded to the duration of a quarter (90 degrees). The overlap between consecutive segments (quarters) was 11.25 degrees. Subsequently, we applied time-warping to each segment so that its duration after time-warping corresponded exactly to 0.5 s. As a last step, we computed the power for each channel and each segment and rejected a segment, if the power of any channel was outside an interval of five standard-deviations around the mean.

Decoding directions from EEG beta band-power
We used the lambda variant of the Source Power Comodulation (SPoC) algorithm  to decode directions from the 20-24 Hz band-power. The goal of SPoC is to estimate a latent signal by linearly combining the recorded EEG signals x[t] such that the covariance between the power of the latent signal and a specified target signal z[t] is maximal (yellow inset in figure 2(c)). This approach is compliant with the assumption that EEG signals (typically below 1 kHz) on the scalp are a linear combination of neural sources (Nunez and Srinivasan 2006). Hence, the patterns associated with the SPoC weights are meaningful with regard to neurophysiology .
The UP and RIGHT components of the hand movement served as target signals z UP [t] and z RIGHT [t]. For each target signal, we fitted one SPoC model. Estimating the power of a signal requires a temporal window. We used the 0.5 s segments, which corresponded to quarters of the circular movement. For each 0.5 s segment e, SPoC required a scalar target value z[e]; we set the value to the average of the target signal within the segment (yellow inset in figure 2(c)). SPoC also required the target values to follow a standard normal distribution (zero-mean, unit-variance). Before fitting the model, we z-scored them. To reduce the sensitivity of SPoC to noise, we reduced the dimensionality of the predictor variables x[t] (EEG channels) to N, with N being the number of dimensions that explain 99% of the variance of x[t]. Then SPoC was applied to fit the model weights w SPoC .
The SPoC decoded signal z SPoC [t] was computed after the weights w SPoC were applied to the bandpass filtered and epoched EEG signals, the power was estimated with a 0.5 s sliding window and the z-scoring was reverted. To achieve a similar SNR improvement as in the low-frequency domain, we applied the same high-pass (0.3 Hz cut-off frequency) and low-pass (0.8 Hz cut-off frequency) filters to the band-power signal. Pearson correlation coefficients between the actual z[t] and SPoC decoded signal z SPoC [t] were used as evaluation metric.
As in the time domain analysis, we assessed within participant significance of the correlations by randomly rotating the 2D hand position signal of each trial. The significance threshold was set to the 95 percentile of 1000 repetitions.
We also computed patterns a SPoC associated to the extracted weights, and scaled them with the decoded signal's standard-deviation to convert the patterns to µV at the scalp level (yellow inset in figure 2(c)). We then applied the same EEG source imaging approach as in the time domain analysis to project the SPoC patterns from the scalp to the cortical surface.

Time-frequency analysis
After epoching the pre-processed data, the previously marked outlier trials were rejected. We applied the SPoC weights to the EEG signals to get the broadband activity of the latent UP and RIGHT components. The latent components were transformed into a time-frequency (TF) representation, using complex Morlet wavelets (Morlet et al 1982) for frequencies ranging from 2 to 70 Hz in 1 Hz steps. We traded a coarser temporal resolution of higher frequencies for a finer spectral resolution by linearly decreasing the full width half maximum of the Gaussian kernel from 1.5 s at 2 Hz to 0.375 s at 70 Hz, since we investigated slow power modulations in the frequency range of the movement (0.5 Hz). After wavelet convolution, we computed the instantaneous power and resampled the result at 20 Hz. For each trial, we time-warped the instantaneous power signals of the first seven complete cycles to their nominal duration. Trials whose first seven cycles were incomplete were excluded. The time-warped instantaneous power was then averaged over trials.
The trial averaged instantaneous power signals were transformed to an event related de-/synchronization (ERDS) (Pfurtscheller and Lopes da Silva, 1999) representations. In detail, the ERDS representation was computed for each channel in decibel (dB) by dividing the trial averaged power by the average (trials and samples) power during rest trials, applying the base ten logarithm and multiplying the result by ten. A low-pass filter (5 s triangularwindow, zero-phase) was used to separate sustained and modulated ERDS activity (Seeber et al 2016). The ERDS activity (combined, sustained and modulated) was then averaged over participants to estimate the group level response.

Behavioral results
After extracting the participant specific 2D movement plane and time-warping the quarters, we computed the participant specific and grand average cycles ( figure 2(b)). The 2D movement plane explained on average 98% of the variance of the 3D movement, with the UP component explaining 62% and the RIGHT component 36% (supplementary table 1). The difference in explained variance was significant (two-sided, paired Wilcoxon signed rank test; p-value = 0.000488). That is, the participants performed an ellipse, with the semi-major axis oriented along the UP direction, rather than a circle.

Low-frequency time domain analysis
The time domain analysis was confined to the frequency range of the movement. In figure 2(d), we show the cycle-averaged PLS UP and RIGHT components at the participant and grand average level. The instantaneous directions within the grand average PLS estimated cycle matched the grand average recorded cycle, while the size of the PLS estimated cycle was smaller (figure 2(b), figures 3(a) and (e)).
We report Pearson correlation coefficients (CCs) to summarize the linear relationship between the recorded and decoded UP (figures 3(b)) and RIGHT (figure 3(f)) components. All CCs were above the participant-specific significance level, which was estimated with the shuffling approach. The distribution over participants is summarized by their probability density. The CCs of both components ranged from approx. 0.5 to 0.85. Regarding the RIGHT component (figure 3(f)), we observed a skewed distribution with its mode at 0.8. The CCs are also listed in supplementary table 2.
It is important to know not only how well the decoded signals predicted the UP/RIGHT component, but also which sources contributed to the result. Using EEG source imaging (Michel et al 2004), we projected the participant-specific patterns to the cortical surface (supplementary figure 4). At the cortical surface, we computed the grand average patterns for the UP (figures 2(f) and 3(c)) and RIGHT (figures 2(f) and 3(g)) components. The patterns indicate where the estimated components were encoded. For the circular movement, they resemble the cycle-averaged low-frequency activity at the maximum along the UP (t = 0 s) and RIGHT (t = 0.5 s) directions (supplementary figure 5b).
To identify voxels with significant pattern activity, we computed two-sided, permutation paired t-tests between the patterns and the estimated baseline. Voxels with significant pattern activity for the UP component are displayed in figure 3(d). They cover mainly contralateral pre-motor, primary sensorimotor and posterior parietal areas. Regarding the RIGHT component, we observed significant pattern activity in contralateral pre-and primary motor areas, in the occipital lobe and lateral areas of the temporal lobes ( figure 3(h)).

Beta band-power domain analysis
The cycle-averaged estimates of the SPoC UP and RIGHT components for each participant and the grand average are shown in figure 2(e). The SPoC estimates clearly followed the circular movement, while their sizes were considerably smaller than the recorded ones ( figure 2(b)).
The CCs between the recorded UP and RIGHT components and their SPoC based estimates are shown in figures 4(a) and (c). The CCs were significant for 10 out of 12 participants for the UP component, and eight out of 12 participants for the RIGHT component.  To identify the cortical sources that contributed to the SPoC components, we computed participantspecific patterns and projected them to the cortical surface (supplementary figure 6). At the cortical surface, we computed the grand average patterns (figures 2(g), 4(b) and (d)). The grand average SPoC UP component pattern ( figure 4(b)) exhibited the strongest activation in parietal and primary sensorimotor areas, with a slight lateralization towards the ipsilateral hemisphere, while the SPoC RIGHT component pattern ( figure 4(d)) was the strongest in the ipsilateral parietal cortex.
We investigated the significance of the SPoC patterns at the level of regions of interest (ROIs). Based on the findings of previous studies (Schaal et al 2004, Culham and Valyear 2006, Seeber et al 2016, Marty et al 2018, we defined eight ROIs that covered the intraparietal sulcus (IPS), superior parietal lobule (SPL), premotor (PM) and primary sensorimotor (SM1) areas of both hemispheres ( figure 5(a)). We tested for differences in pattern norms between the SPoC UP (figures 5(b)) and RIGHT (figure 5(c)) components against baseline. We observed significant pattern norms for the SPoC UP component in the IPS, SPL and SM1 ROIs of both hemispheres (table 1). Regarding the SPoC RIGHT component, we observed a significant effect for the right IPS, while the effect at the right SPL did not reach significance (table 1).

Time-frequency analysis
To explore whether the power modulations of the SPoC components were only confined to the beta band (20-24 Hz), we applied the SPoC weights to the broad-band EEG signals. The resulting signals were transformed into the time-frequency domain. In the time-frequency domain, we analyzed the changes in power relative to rest trials by computing event-related de-/synchronization (ERDS) maps (Graimann et al 2002). Figure 6 displays the grand average ERDS maps for the SPoC UP and RIGHT components. We decomposed the ERDS activity into sustained (figures 6(b) and (f)) and modulated (figures 6(c) and (g)) effects. Regarding the SPoC UP component, a cyclic co-modulation in the beta (20-24 Hz) band was apparent (figures 6(a) and (c)). It was also specific to the beta band (figure 6(c)). Figure 6(d) shows the beta band modulation in detail. The grand average modulation effect was tightly coupled with the recorded UP component of the cyclic arm movements (t ≥ 2.0 s) and exhibited less variance across participants than the sustained effect. Regarding the sustained effect ( figure 6(b)), we observed an average ERD of 1.8 dB in the beta (20-24 Hz) band during the movement phase, and additionally an average ERD of 2.3 dB in the mu (10-12 Hz) band. Regarding the SPoC RIGHT component (figures 6(e)-(h)), we also observed a tightly coupled modulation effect specific to the beta band (figures 6(g) and (h)). The sustained ERD in the beta (1.0 dB) and mu (1.0 dB) bands (figures 6(f) and (h)) was weaker compared to the SPoC UP component.
We observed the strongest sustained beta-band ERD (2.4 dB) at the EEG channel closest to the contralateral rolandic region (channel C3; supplementary figures 7(a) and (b)). The sustained ERD at the channel closest to the ipsilateral rolandic region was 2.0 dB on average (channel C4; supplementary figure 7(a)). We did not observe co-modulations with the arm movement in the ERDS maps of any EEG channel.

Discussion
We found that the human EEG signals in the lowfrequency and in the beta domain carry directional information about the continuously moving arm. We applied machine learning techniques, which comply with the linear mixing of sources model, to reconstruct the arm movement from EEG activity in either domain. At each time instance, the hand position and its reconstruction-the non-invasively decoded population vector-were characterized by two directional components. Using EEG source imaging, we recovered the brain sources associated to the two directional components. In the low-frequency domain, the brain sources were located in motor and parietal areas contralateral to the moving arm.
In the beta domain, the brain sources were located in primary sensorimotor and parietal areas of both hemispheres.
The cortical representations of continuous movement kinematics have been investigated in different frequency ranges using various modalities to acquire neural information (Moran and Schwartz 1999, Georgopoulos et al 2005, Schalk et al 2007, Kim et al 2015, Kobler et al 2018. Here, we investigated the non-invasive EEG activity as participants performed continuous, circular arm movements. We hypothesized that arm movement direction is encoded in both, low and high frequency domains. We formulated this hypothesis following the study of Lv et al (2010), in which they state that both lowfrequency (0.1-4 Hz) and beta band (12-30 Hz) components encode information about arm movement velocities. Our approach was to infer from each of the two frequency domains not only the latent EEG components-the non-invasively decoded population vector-but also their cortical representation.
In the low-frequency domain (0.3-0.8 Hz), we studied the cortical representation in terms of activity phase locked to the 2D hand movement. We found a contralateral activation in primary sensorimotor and posterior parietal areas (figures 3(d) and (h), supplementary figures 6(a) and (b)). Similar locations have been reported previously in a related study (Kim et al 2015). They showed that the infinity-like trajectory of 3D arm movements can be decoded from EEG activity. Specifically, they found that a low-frequency component encodes information about the hand movement velocity, as well as elbow velocity. Other studies (Schalk et al 2007, Pistohl et al 2008 demonstrated that circular hand movement trajectories can be approximately inferred from ECoG signals at electrodes covering the contralateral motor cortex. They obtained an average correlation for 2D position of approximately 0.5. In the current study, we obtained correlations for direction in the range between 0.5 and 0.85 (figures 3(b) and (f), supplementary table 2). It is important to note that the continuous movement tasks in the previous studies were slightly different than ours, so a direct comparison is not straightforward.
In the higher frequency domain, we focused on the beta band (20-24 Hz), which has been previously reported to encode non-directional (Hammer et al 2016, Marty et al 2018 and directional kinematic information (Schalk et al 2007, Lv et al 2010, Korik et al 2018. Specifically, we studied power modulations of latent EEG components. Using the SPoC algorithm  we successfully estimated beta-band latent components whose power co-modulated with arm movement direction (figures 2(e) and (6)). At single cycle level, the average correlations were 0.30 for the UP and 0.25 for the RIGHT component (figures 4(a) and (c)). They were above the significance level for 10 out of 12 participants for the UP and 8 out of 12 participants for the RIGHT component.
Using EEG source imaging, we found that the brain sources, whose beta band-power co-modulated with the movement, covered primary sensorimotor and posterior parietal areas. Similar findings in the beta band have been reported by Seeber et al (2016). They found, in an EEG source imaging study that in addition to sustained ERD, the mu (10-12 Hz) and beta (18-24 Hz) band-power in the contra-and ipsilateral primary sensorimotor areas were co-modulated with the phase of a rhythmic hand opening and closing task. We also observed a sustained ERD effect of the SPoC UP and RIGHT components in the mu (10-12 Hz) band (figures 6(b) and (f)). However, we did not observe co-modulations with movement direction in the mu band (figures 6(c) and (g)). Since SPoC was fit to the beta band, this discrepancy could be explained twofold. First, if the co-modulating sources in the beta and mu band had a different spatial configuration, SPoC would identify the beta band source. Second, if other non-task related sources were active  in the mu band, SPoC would not effectively attenuate them. In either case, the result would be a weaker co-modulation in the mu band compared to the beta band.
In a recent MEG source imaging study (Marty et al 2018) that investigated repetitive hand pinching movements, Marty et al found that the coherence between the amplitude of beta oscillations (15-25 Hz) and the 1D acceleration signal coming from forefinger movements emerged in bilateral primary sensorimotor (SM1) cortex. In contrast to Seeber et al and Marty et al who investigated 1D signals, we identified two sources coupled to the 2D arm movement direction. Other approaches to study beta band oscillatory components have been reported in several studies (Chao et al 2010, Bundy et al 2016, Korik et al 2018. However, in these studies the band-power was estimated at the channel level before latent components were identified. Dähne et al demonstrated that this approach disregards the linear mixing of sources and, therefore, hinders the interpretability of the model patterns . EEG/MEG source imaging techniques applied in (Seeber et al 2016, Marty et al 2018 respect the linear mixing of sources model. Compared to SPoC, they require additional anatomical information in terms of a forward model. Interaction with the environment requires flexible processing of body-related and spatial information derived through different sensory channels. For instance, movement direction can be conveyed both by vision and proprioception. It has been shown that multisensory information about the environment and the body converges in the posterior parietal cortex (PPC) (Gallivan and Culham 2015). In the current study we occluded the visual feedback coming from the circular arm movement and focused on proprioception. In general, processing in PPC has been related to the transformation from sensory to motor-related references frames (Culham and Valyear 2006, Scott 2012, Vesia and Crawford 2012, Gallivan and Culham 2015. In a recent review, Battaglia Mayer (Battaglia-Mayer 2019) discusses that information related to the same process (e.g. directional information in reaching movement) is encoded in multiple, simultaneous neural representations. It is likely that the activation in PPC and SM1 both in the low-frequency domain (figure 2(d)) and beta domain (figure 2(e)) reflect the integration and transformation of directional information. In the low frequencies, the directional information was mainly expressed in SM1 and PPC contralateral to the moving arm, while in the beta range SM1 and PPC of both hemispheres predicted direction. We interpret that the low-frequency and beta domain activity might reflect direction related activity of distinct neural representations.
Based on the neural representation of the two components both in the low-frequency domain and the beta domain, we observed that the UP component had stronger activations relative to the RIGHT component. From the behavioral analysis, we found that the participants performed an ellipse rather than a circle with semi-major axis along the UP direction. On a closer look on the participant-specific 2D movement planes (supplementary figure 1), we observed that the 2D plane orientations varied more along the RIGHT direction than the UP direction. Hence, we believe that there are two potential explanations for a stronger representation of the UP component in the group level decoder patterns. First, the magnitude of the vector on the vertical direction was larger due to the shape of the movement (ellipse, rather than a circle). Second, the movement planes were more consistent along the vertical direction than along the horizontal direction (supplementary figure 1). That is, the orientation of the horizontal direction in the 2D movement plane varied more in the original 3D space, which could have led to different cortical activation across participants (Caminiti et al 1990) and, therefore, to a weaker pattern at the group level.
In terms of correlation between the recorded and decoded UP and RIGHT components, we observed higher correlations at the single cycle level in the low-frequency domain (figures 3(b) and (f)) compared to the beta domain (figures 4(a) and (c)). This observation is consistent with findings of previous EEG and ECoG studies who investigated both domains (Schalk et al 2007, Waldert et al 2008, Bundy et al 2016, Hammer et al 2016. Despite similar processing pipelines (figure 2(c)), the obtained correlations could potentially depend on the different methods (PLS regression, SPoC), limiting a direct comparison of the two domains. Moreover, the correlation metric is scale invariant; it does not capture amplitude mismatches (Antelis et al 2013). We observed mismatches in the average cycles for both frequency domains, with a larger mismatch in the beta domain (figures 2(b), (d) and (e)). Both algorithms (PLS, SPoC) were applied to estimate the UP/RIGHT components at the single cycle level. Averaging over repetitions typically improves the SNR. The observed amplitude mismatches in the average cycles can, therefore, be attributed to the lower SNR at the single cycle level.
Although we observed for all the participants in our study a rotational shift of the brain patterns in meaningful hand-related areas over the course of the movement (figures 2(f) and (g); supplementary figures 6(c) and (d)), we also observed patterns at temporal and occipital sites in the low-frequency domain (figures 3(d) and (h)). Since there were no visual stimuli (except the fixation cross) present during the movement condition and the arm movement was occluded, we did not expect activation of occipital cortex. The combination of rejecting independent components with SASICA and fewer electrodes at inferior locations might have led to residual movement artifacts being projected to temporal and occipital sites. Other potential limitations consist in the variability of the orientation of the 2D plane among participants (supplementary figure 1) and in the difference in visual feedback for the first participant, who could observe his/her circular arm movements.

Conclusion
Recording non-invasively from a large neural population allowed us to extract arm movement direction related activity on a short timescale. Using machine learning and source imaging techniques, we revealed the population vector with a high congruence to the kinematic behavior from EEG activity in temporal and spectral domains. The ability to interpret a fast change in neural activity with a high temporal resolution is not only useful but essential for a reliable and continuous control of robotic limbs or neuroprostheses. The non-invasively decoded population vector may further provide a versatile concept for a quantitative understanding of dynamic coding of motion within large scale brain networks.