Distance- and speed-informed kinematics decoding improves M/EEG based upper-limb movement decoder accuracy

Objective. One of the main goals in brain–computer interface (BCI) research is the replacement or restoration of lost function in individuals with paralysis. One line of research investigates the inference of movement kinematics from brain activity during different volitional states. A growing number of electroencephalography (EEG) and magnetoencephalography (MEG) studies suggest that information about directional (e.g. velocity) and nondirectional (e.g. speed) movement kinematics is accessible noninvasively. We sought to assess if the neural information associated with both types of kinematics can be combined to improve the decoding accuracy. Approach. In an offline analysis, we reanalyzed the data of two previous experiments containing the recordings of 34 healthy participants (15 EEG, 19 MEG). We decoded 2D movement trajectories from low-frequency M/EEG signals in executed and observed tracking movements, and compared the accuracy of an unscented Kalman filter (UKF) that explicitly modeled the nonlinear relation between directional and nondirectional kinematics to the accuracies of linear Kalman (KF) and Wiener filters which did not combine both types of kinematics. Main results. At the group level, posterior-parietal and parieto-occipital (executed and observed movements) and sensorimotor areas (executed movements) encoded kinematic information. Correlations between the recorded position and velocity trajectories and the UKF decoded ones were on average 0.49 during executed and 0.36 during observed movements. Compared to the other filters, the UKF could achieve the best trade-off between maximizing the signal to noise ratio and minimizing the amplitude mismatch between the recorded and decoded trajectories. Significance. We present direct evidence that directional and nondirectional kinematic information is simultaneously detectable in low-frequency M/EEG signals. Moreover, combining directional and nondirectional kinematic information significantly improves the decoding accuracy upon a linear KF.


Introduction
Inferring movement parameters from electrophysiological brain signals during different volitional states has been a central topic in brain-computer interface (BCI) research [1,2]. BCIs potentially allow individuals with tetraplegia to control computer cursors, robotic arms or neuroprostheses and thereby replace or restore lost or compromised 5 Author to whom any correspondence should be addressed. function [3][4][5][6]. In nonhuman primates, intracortical BCIs (iBCIs) typically decode kinematic and kinetic information of executed movements from neural spiking activity in primary motor cortex (M1), dorsal premotor cortex (PMd), primary somatosensory cortex (S1) and posterior parietal cortex (PP) [7][8][9][10]. Decoding kinematic parameters from attempted movements has allowed paralyzed individuals, equipped with an iBCI, to control end-effectors with a high accuracy [11][12][13][14]. Despite non-invasive functional neuroimaging techniques such as electroencephalography (EEG) and magnetoencephalography (MEG) have a lower spatial resolution and signal to noise ratio (SNR), they also carry information about movement kinematics and kinetics [15][16][17][18][19]. In this context, low-frequency activity in the range of the movement has been reported to encode kinematic information during different volitional states, such as executed [17,19,20], imagined [21][22][23] and observed movements [20,23]. Although numerous studies reported that kinematic information is encoded in M/EEG signals, it is not clear with which accuracy a non-invasive BCI can infer movement trajectories.
In studies that decode directional kinematics (e.g. position and velocity) from low-frequency activity, the accuracy of the decoded trajectories is typically quantified by computing correlations between the recorded and decoded trajectories [15,19,24]. Correlation is invariant to amplitude mismatches, rendering it insufficient to qualitatively evaluate the accuracy [25]. In fact, linear regression based decoders that minimize the mean squared error (MSE) tend to exhibit amplitude mismatches in discrete reaching [26,27] and continuous movement tasks [15,20,24,28]. The decoded trajectories are typically estimated too small, with the mismatch increasing as the correlation decreases [29,30]. Nonlinear methods reduce the amplitude mismatch [23], indicating that the low-frequency activity could carry kinematic information that is nonlinearly related to the directional kinematics.
There is evidence that the neural activity carries information about distance (length of the position vector) and speed (length of the velocity vector) [31]. Distance and speed are inherently nondirectional and nonlinearly related to position and velocity. They are also uncorrelated to position and velocity during commonly studied tasks (e.g. center out reaching, pursuit tracking). The nondirectional kinematic information is known to be encoded in neural spiking activity in M1 [31,32], in low-frequency local field potentials [33] and electrocorticographic (ECoG) activity [34]. Jerbi et al were the first to report speed related tuning effects in low-frequency MEG activity [16]. Using a continuous movement task, Hammer et al demonstrated that low-frequency ECoG activity simultaneously carries information about velocity and speed [34]. We recently corroborated their findings in low-frequency MEG activity and a 2D pursuit tracking task (PTT) [35].
Due to the nonlinear relationship between directional and nondirectional kinematics, a linear model can not combine the information encoded in the low-frequency M/EEG activity. We surmise that a decoding algorithm that models their nonlinear relation should achieve superior decoding accuracy and reduce the amplitude mismatch. Inspired by models used in iBCIs [32,36], we use an unscented Kalman filter (UKF) [37] to combine the directional and nondirectional information during a 2D PTT. In an offline analysis, we compared the UKF decoding accuracy with the one achieved by two linear models, a Kalman filter (KF) [38] and a Wiener filter (WF) [39].

Participants
In this study we reanalyzed experimental data obtained from two distinct groups of healthy human participants during executed and observed movements. The first group (experiment 1; 15 participants, nine female, six male, mean age of 23.8 years) was previously published in [20,40], and the second group (experiment 2; 23 recruited participants, 18 male, five female, mean age of 28.5 years) was previously published in [35]. As reported in [35], the data of four participants in the second group were excluded because of incomplete recordings or incorrect positioning in the MEG scanner. In both groups, all participants had normal or corrected to normal vision and self-reported to be right-handed.
The experimental procedure conformed to the Declaration of Helsinki and was approved by the ethics committee of the Medical University of Graz (group 1) and Osaka University Hospital (group 2). After the study participants were informed about the purpose and procedure of the study, they gave their written consent.

Experimental paradigm
The participants of experiment 1 underwent two experimental conditions (execution and observation) while their EEG activity was recorded. The experimental paradigm and setup have been described in detail in [20]. Figure 1(a) illustrates the sequence of events during a trial (top) and the virtual 2D environment displayed on a computer screen (bottom). In 180 trials (90 per condition, pseudo-randomly distributed), the participants completed a center-out task and a 16-s long PTT. In this study, we analyzed the PTT.
In the execution condition, we instructed the participants to smoothly track a moving target stimulus (large green ball) with a cursor (small gray ball) and their gaze. In the observation condition (blue target), the participants tracked the target only with their gaze; they lost control of the cursor and were instructed to maintain a resting position. In this condition, the paradigm replayed a previous, matching cursor trajectory.
The 2D target trajectories consisted of two independent and identically distributed (iid) components (horizontal and vertical). Each component was sampled from band-limited pink noise (0.3 to 0.6 Hz). The target trajectories were generated offline and were identical for all participants. The trajectory Figure 1. Overview of the experimental approach and methods. Data from two previous experiments was analyzed. In the first, EEG was recorded during right arm movements. In the second, MEG was recorded during right index finger movements. (a) Experimental paradigm. The sequence of events during a trial are displayed at the top. The paradigm distinguished two conditions (execution and observation). In the execution condition, the participants controlled a cursor in a virtual 2D environment on a screen. During a pursuit tracking task (PTT) they used the cursor to track a target stimulus for 16 s. In the observation condition, their hand/finger was resting; the paradigm replayed a matching cursor trajectory. (b) Mapping of right hand movements (experiment 1) and right index finger movements (experiment 2) to the 2D virtual environment. (c) Decoding algorithm evaluation approach. Three algorithms, namely, a partial least squares (PLS) regression based Wiener filter (WF), a Kalman filter (KF) and an unscented Kalman filter (UKF) were fit to decode cursor kinematics. (d) Experiment 1. Example of test set trajectories (recorded and decoded) during the PTT for a representative participant. Recorded target (yellow) and cursor (gray) position trajectories for the horizontal (top) and vertical (bottom) component. The decoded cursor trajectories (WF in red, KF in green, UKF in blue) are plotted on top of the recorded cursor trajectories (gray). (e) As in (d) for the second experiment. generation procedure ensured that the positions and velocities of the horizontal and vertical components were decorrelated, and in addition that the position and velocity of the same component were decorrelated at lag 0 [20].
In experiment 1, the participants controlled the cursor with 2D right arm movements on a table. Figure 1(b) illustrates how arm movements were mapped to cursor movements in the virtual environment on the screen. If the palm of the right hand moved forward (or rightward) by 5 cm from a preregistered resting position, the cursor moved upward (or rightward) to the edge of the virtual environment (50% of the workspace). The online transformation of hand movements into cursor movements introduced an average delay of 55 ms with a standard deviation (SD) of 1 ms.
The participants of experiment 2 underwent a similar experimental paradigm while their MEG activity was recorded. In this experiment, the participants were lying in an MEG scanner. The experimental conditions, paradigm and setup have been described in detail in [35]. The sequence of events during a trial were identical to experiment 1, except that we omitted the center-out task ( figure 1(a)). We recorded 160 trials (80 per condition, pseudo-randomly distributed) per participant.
In experiment 2, the participants performed 2D right index finger movements on a planar surface to the side of their body. If the tip of the index finger moved forward (or rightward) by 1.5 cm from a preregistered resting position, the cursor moved upward (or rightward) to the edge of the virtual environment (50% of the workspace). We used a custom motion capture system, to transform the index finger movements into cursor movements [35]. The system introduced an average delay of 190 ms.

Data acquisition
In experiment 1, we recorded EEG activity with 64 active electrodes (actiCAP, Brain Products GmbH, Germany). The electrodes were placed according to the 10-10 system; reference and ground were placed at the right mastoid and AFz. We recorded electrooculographic (EOG) activity with six additional active electrodes, placed next to the outer canthi and above and below both eyes. EEG and EOG signals were amplified and digitized at a rate of 1 kHz with a biosignal amplifier (BrainAmp, Brain Products GmbH, Germany).
In experiment 2, we recorded MEG activity with a 160 channel whole-head MEG system (MEGvision NEO, RICOH Ltd., Japan) housed in a magnetically shielded room. EOG activity was recorded with four electrodes placed at the outer canthi of both eyes, and above and below the left eye. EOG signals were recorded with a biosignal amplifier (Neurofax EEG 1200, Nihon Koden Corp., Japan). The MEG and EOG signals were recorded synchronously at a rate of 1 kHz. Inside the MEG scanner, the participant's head was resting on a cushion. We asked the participants to avoid head and shoulder movements during the course of the experiment. At the beginning of each run (16 trials), we additionally measured the head position with five marker coils, attached to the forehead.
Visual stimuli (60 Hz sampling rate) and motion capture signals (experiment 1: 100 Hz sampling rate; experiment 2: 60 Hz sampling rate) were recorded with the labstreaming layer protocol [41]. The paradigm was implemented in Python (version 2.7) and the simulation and neuroscience application platform [42].

Data preprocessing
Offline, we analyzed the data with a custom pipeline implemented in Matlab (version 2015b; Mathworks Inc. USA), EEGLAB toolbox (version 14.1.1) [43], Brainstorm (version 5 June 2018) [44] and R [45,46]. As a first step, we synchronized the M/EEG and EOG signals with the visual stimuli by aligning impulses captured by a photodiode with the associated markers. We resampled the synchronized signals at a rate of 200 Hz, before we separately preprocessed the cursor kinematics and the M/EEG signals. Supplementary figure 16 (available online at https://stacks.iop.org/JNE/17/056027/mmedia) summarizes the M/EEG preprocessing steps in both experiments.

EEG preprocessing (experiment 1)
The EEG and EOG signals were first high-pass (0.25 Hz cut-off frequency, Butterworth filter, eighth order, zero-phase) and notch (49 and 51 Hz cut-off frequencies, Butterworth filter, fourth order, zerophase) filtered to attenuate drifts and powerline noise. We visually inspected the signals to detect bad channels. A channel was considered as bad, if it was affected regularly by pops, drifts or muscle artifacts. On average, 2.1 channels were marked as bad. The signal of the bad channels was spherically interpolated using EEGLAB. Next, we applied the sparse generalized eye artifact subspace subtraction algorithm (SGEYESUB) to attenuate eye movement and blink related artifacts [47]. SGEYESUB was fit to calibration data recorded at the start and end of the experiment according to the paradigm presented in [48]. After the eye artifacts were attenuated, we excluded the EOG and fronto-temporal channels (FPz, AF7, AF8, AF8, FT9, F7, F8, FT10) from the further analysis and re-referenced the remaining 57 EEG channels to their common average.

MEG preprocessing (experiment 2)
We excluded 31 MEG channels at inferior temporal locations because they were contaminated by technical artifacts. To compensate for small head movements across runs, we spherically interpolated the MEG signals of the remaining 129 channels to their average location across runs. The mean distance of a channel to its average location across runs was 5 mm (0.4 mm SD). After concatenating the signals across runs, we also applied high-pass (0.25 Hz cut-off frequency, Butterworth filter, eighth order, zero-phase) and notch (59 and 61 Hz cut-off frequencies, Butterworth filter, fourth order, zero-phase) filters. To further attenuate technical and spatially stationary artifacts, we applied independent component analysis [49]. Using the extended Infomax algorithm, we decomposed the MEG signals (high-pass filtered at 0.4 Hz, instead of 0.25 Hz) into independent components (ICs) that explain 99.9% of the variance. We marked 8.6 (0.9 SD) out of 63.5 (0.1 SD) ICs (visual inspection). The marked ICs were then removed from the 0.25 Hz high-pass filtered signals. As in the EEG preprocessing, we applied SGEYESUB to attenuate eye movement and blink related artifacts and excluded the EOG and fronto-temporal channels (RT21, RT22, RT31, RT32, RT33, LT11, LT21, LT22, LT33) from further analysis.

Common M/EEG preprocessing steps
All subsequent processing steps were shared between experiment 1 and 2. The first common step was to attenuate transient drifts and pops in the continuous M/EEG signals with the high-variance electrode artifact removal algorithm [50].
Next, we detected bad epochs [51]. We first epoched the continuous, broadband data into 14-s epochs. Each epoch started 1 s after the beginning of the PTT. An epoch was considered as bad, if the broadband M/EEG activity of any channel exceeded a threshold (±200 µV or ± 5pT), had an abnormal probability, variance or kurtosis (more than 6 (probability), 5 (variance), 6 (kurtosis) SDs around the mean), the correlation between the EOG and target position signals were improbable (exceeding 4 SDs around the mean), or a tracking error happened (e.g. jerky or no cursor movement). We tested for abnormal probability, variance and kurtosis twice, in order to detect gross outliers in the first iteration and subtle outliers in the second iteration. All criteria combined, a total of 29.5 (9.8 SD) of 180.7 (2.8 SD) epochs were marked for rejection in experiment 1, and 27.3 (8.9 SD) of 152.2 (13.1 SD) in experiment 2.
To extract low-frequency activity, we applied a low-pass filter (2 Hz cut-off frequency, Butterworth filter, sixth order, zero-phase) to the continuous M/EEG signals and resampled the result at 10 Hz. In the last preprocessing step, we extracted the 14-s epochs and rejected the previously marked bad epochs.

Cursor kinematics preprocessing
We used Savitzky-Golay filters [52] (finite impulse response filter, third order polynomial, 21 filter taps, zero-phase) to compute smoothed cursor positions and velocities. The origin of the coordinate system was the center of the virtual 2D environment ( figure 1(b)). Similar to the M/EEG signals, we resampled the results at 10 Hz and extracted the 14-s epochs. We computed distance and speed by computing the Euclidean norm of the 2D cursor position and velocity vectors at each sample t. As a last step, we zscored the cursor kinematics (2D position, 2D velocity, distance and speed).

Encoding of cursor kinematics in low-frequency M/EEG activity
We fitted a general linear model (GLM) with the cursor kinematics as predictors, to identify how the cursor kinematics were expressed in the lowfrequency M/EEG activity at 11 lags l ∈ [−0.5s, 0.5s]. In detail, we fitted one GLM per experimental condition and lag l. By concatenating all observations (samples and trials) of one condition, the kinematics coulbed expressed as a 6 × n observations matrix K , d and s coding the horizontal position, horizontal velocity, vertical position, vertical velocity, distance and speed. For each lag l, we defined the GLM as with an n channels × n observations matrix X (l) containing the EEG activity at lag l, an n channels × 6 matrix A (l) reflecting the regression coefficients and an n channels × n observations noise matrix E (l) . The least squares estimate of A (l) is with an n channels × 6 cross-covariance matrix Σ X K (l) between the low-frequency M/EEG activity and the cursor kinematics, and a 6 × 6 covariance matrix Σ K K between the cursor kinematics. We used analytical shrinkage regularization [53] to estimate ΣKK.

M/EEG source imaging
As in [20,35], we applied M/EEG source imaging [54,55] to project the regression coefficientsÂ (l) to the cortical surface of the ICBM152 template boundary element (BEM) head model [56]. The BEM comprised three layers (cortex, skull, scalp) with relative conductivities (1, 0.008, 1) in experiment 1, and (1, 0, 0) in experiment 2. The cortex layer contained 5011 voxels. We co-registered the BEM and the digitized head points (including EEG and EOG electrode locations) using three anatomical landmarks (nasion, left and right preauricular points). To compensate deviations between participant and template anatomy, we finalized the co-registration by projecting floating EEG electrodes to the scalp layer. OpenMEEG [57,58] was used to compute the forward model for three orthogonal sources per voxel. sLORETA [59] was used to estimate the inverse solution. To reduce the effect of sensor noise in the inverse solution, resting data (similar preprocessing as the M/EEG data; recorded at the start and end of the experiments) were used to estimate a sensor noise covariance matrix (shrinkage regularization with 10% of its mean eigenvalue). After we projected the participant-specific regression coefficients to the cortical surface, we scaled the voxel level regression coefficients with the inverse global field power (GFP) [29]. We estimated the GFP by randomly selecting one sample from all trials, averaging across these samples, and computing the SD across voxels. This procedure was repeated 10 000 times. The final participant-specific GFP estimate was set to the median of the repetitions.

Source space statistics
We first tested for consistent (effect with the same sign across participants) expression of kinematic information at the group level. In detail, we tested the regression coefficients at the voxel level using one-sided permutation tests [60,61]. We obtained the observed voxel activation by first averaging over participants and then computing the norm across the three voxel components. To reduce the number of tests, we averaged the voxel activations associated with either position, velocity or amplitude (distance and speed) predictors. For example, the voxel activations of the horizontal and vertical positions were averaged. In doing so, we identified voxels that carried information about the vertical and/or horizontal position. Next, we used 10 000 permutations to obtain a random distribution. In each permutation, we randomly flipped the signs of a voxel's regression coefficients before averaging over participants. The p-value was then the fraction of random voxel activations that were larger or equal to the observed voxel activation. We computed tests for all voxels (5011), lags (11), type of kinematic signal (3; position, velocity and amplitude) and condition (2; execution and observation).
In a second test, we compared the two conditions. We specifically tested for voxels with stronger activation in either condition (i.e. consistently higher voxel norm in one condition compared to the other). We computed two-sided, permutation, paired t-tests [60,61] for all voxels (5011), lags (11) and type of kinematic signal (3).
The total number of tests for consistent activation and difference across participants was 496 089 in both experiments. The significance level of 0.05 was corrected for multiple comparisons in both experiments by controlling the false discovery rate (FDR) [62].

Cursor kinematics decoding
In a causal sliding window approach, we decoded the cursor kinematics from the low-frequency M/EEG activity. We compared three decoding algorithms, namely, a WF, a KF and an UKF. Figure 2 depicts the models underlying the three decoding algorithms. In all algorithms, we combined the M/EEG activity x t at the current time-point t and previous time-points (= lags) as an n channels · n lags × 1 feature vector x t = [x t ; x t−1 ; ]. The number of lags n lags used to create the feature vector defined the model order. For example, if only the previous lag was used, the model order was 1. We then projected the high-dimensional features x t to a latent subspace y t

Decoding algorithms
with an n dim × n features projection matrix L. As in [20,22], the SIMPLS algorithm [63] was used to estimate L so that the latent subspace y t contained the activity ofx t that maximally explained the variance of the zscored cursor kinematics k t In a nutshell, SIMPLS iteratively extracts the 1D subspace that explains most of the variance of the cross-covariance matrix between the predictor x t and dependent variables k t . This can be done via iterative singular value decomposition of the cross-covariance matrix, extraction of the dominant singular value and its associated projection vectors, Gram-Schmidt orthogonalization with the previously extracted subspaces and deflating the resulting subspace from the cross-covariance matrix. These steps are repeated until the desired stopping criterion is reached (e.g. a fixed number of latent components n dim ).
To reduce collinearity among the features and also the computation time of the KFs, we set n dim ≪ n features . We fixed n dim so that the latent subspace would explain 70% of the variance in the features . For a first order model, n dim was on average 32 (4 SD) in experiment 1 and 29 (2 SD) in experiment 2. The SIMPLS algorithm was also used to compute the coefficients Wof the WF to predict k t from y t (figure 2(a)).
The second algorithm was a KF [38]. We used it to infer the cursor kinematics (= state) from the latent subspace y t and the previous state. Since a KF is a linear state space model ( figure 2(b)), we only decoded z-scored positions and velocities k t The KF treats the noise corrupting successive observations y t as conditionally independent given the hidden state k t . Combining M/EEG activity at multiple lags x t into y t violates this assumption. However, this assumption is already violated for the M/EEG activity at a single lag x t because some noise sources (e.g. brain sources that do not encode kinematic information) are auto-correlated. In the next steps we outline the fitting procedure for a first order model. For higher order models (multiple previous states), we refer to [36]. The 4 × 4 state transition matrix can be estimated from the covariance matrix between k t and k t−1 and the covariance matrix of k t−1 . We used analytical shrinkage regularization [53] to estimate Σ kt−1, kt−1 from calibration data. The covariance matrix Σ qt qt of the state transition error can be estimated from the residuals The KF also relates the kinematic state k t with the latent subspace y t using a 4 × n dim matrix H. H can be determined by minimizing the mean squared error between y t and its prediction After fitting H to calibration data with n observations, we used the residuals r t defined as to estimate the measurement noise covariance matrix Σ rt rt . The third algorithm can be seen as an extension of the KF. We used the linear state transition matrix F from before (figure 2(c)). The key difference was, that we explicitly modelled the nonlinear relationship between position and distance or velocity and speed or both with the function If distance and speed are considered, the function essentially computes the extended kinematics k t , defined in (5). Next, we utilized the linear encoding of distance and/or speed information in the lowfrequency M/EEG signals to predict y t from ex (k t ) Using H and (11), the estimated latent subspace y t is As before, we used the residuals to estimate the measurement noise covariance matrix Σ rt rt . After the parameters of the state space models (figures 2(b),(c)) were fitted to calibration data, we used a square root KF [64] to infer the kinematics in test data for the linear model, and a square root UKF [65] for the nonlinear model.

Algorithm evaluation
We fitted and evaluated the algorithms with a tenfold cross-validation scheme. The model order was fixed to 3 (i.e. M/EEG activity and cursor kinematics of three previous lags). In ten iterations, we used nine folds to fit the parameters, and the held-out fold to test the performance. For the state space models, we estimated the initial state covariance matrix from the calibration data. We then applied the algorithms to each trial within the test set. The KF and UKF required an initial estimate of the state k 0 ; we used the recorded kinematics of the time-points before the first sample in the tested trials. All algorithms decoded the directional kinematics. We used equation (11) to infer the nondirectional kinematics from the directional ones in the test set. Three metrics were used to test the performance. We first computed Pearson correlation coefficients between the recorded test set kinematics and their decoded estimates. The second metric was the SNR as defined in [36]. For a recorded kinematic signal z t (e.g. horizontal position) and its decoded estimateẑ t the SNR is defined as the ratio of the variance of z t over the mean squared error (MSE): The third metric was the decoded signal to signal ratio (DSSR). We defined it as the variance of the decoded signal over the variance of the recorded signal: The DSSR captures the amplitude mismatch between the decoded and recorded signal. A DSSR of 0 dB indicates that the amplitude range of the recorded and decoded signal match. A negative DSSR indicates a bias of the decoder towards too small amplitudes.

Participant level significance
We assessed the significance of the correlation and SNR metrics at the participant level using a shuffling approach. To break the association between the M/EEG activity and the kinematics, we randomly permuted the kinematics across trials, before we split the data into two folds-one train fold with 90% of the data and one test fold with 10% of the data. The two folds were then used to fit and evaluate the algorithms. The shuffling procedure was repeated 1000 times. We set the p-value to the fraction of random results that were higher than the average cross-validation result. We controlled the FDR for 2 (conditions) * 3 (algorithms) * 4 (kinematics) * 2 (metrics) = 48 tests at a significance level of 0.05 [62].

Group level statistics
We compared the decoding accuracy of the algorithms in terms of correlation, SNR and DSSR for each experiment, condition and type of kinematic signal (directional, nondirectional). To identify significant differences, we used two-sided, permutation, paired t-tests (df = 14, 10 000 permutations, significance level = 0.05). The obtained p-values were FDR adjusted for 3 (metrics) * 3 (algorithms) * 2 (conditions) * 2 (type of kinematic signal) = 36 tests.

Tracking behavior
The tracking behavior during the PTT was inferred by computing cross-correlation curves between the target and cursor positions in the execution condition.

Neurophysiology
We used the fitted regression coefficients of the GLMs to identify cursor kinematics related effects in the low-frequency M/EEG activity. The low-frequency EEG signals encoded information about the cursor kinematics during the right arm PTT in experiment 1.
In figure 3, we show voxels whose activity was significantly tuned to specific cursor kinematics at the group level. Figure 3 displays only activity during negative lags; for example, a significant effect at the lag −0.1 s indicates that the EEG 100 ms ago predicted the current cursor kinematics. Supplementary figures 1 and 2 also include positive lags.
In the execution condition, we observed significant amplitude (distance/speed) related activity ( figure 3(a)) in contralateral motor and parietal areas as well as parieto-occipital, occipital and temporal areas of both hemispheres. There was a distinct peak in contralateral primary-and premotor areas at lags −0.1 to 0.0 s. Position ( figure 3(b)) and velocity (figure 3(c)) information were consistently expressed in parieto-occipital areas. Due to the large crosscorrelations between positions and velocities across lags (supplementary figures 14(b), (c)), we observed similar effects in both kinematic variables at different lags. In the case of velocity, the activity was strongest for the lags −0.1 to 0.0 s. Because the crosscorrelation between velocity and position peaked at a delay of 500 ms (supplementary figures 14(b), (c)), we observed a similar effect at the lags 0.4 and 0.5 for position (supplementary figure 1). For an effect to be significant in figure 3, the dipole orientations had to be consistent across participants. If the dipole orientation was allowed to vary across participants, M1, PMd, S1 and PP voxels encoded significant directional information in execution condition (supplementary figure 3).
In the observation condition, parieto-occipital areas also expressed information about position (figure 3(e)) and velocity (figure 3(f)). As before, the effect peaked for velocity at the lags −0.1 to 0.0 s. In contrast to the execution condition, we also observed significant activity at the lags −0.5 to −0.4 s. Regarding amplitude information ( figure 3(d)), we also observed significant activity in parieto-occipital and occipital areas. Compared to the execution condition, we neither observed sensorimotor activity at lags −0.1 to 0.0 s, nor activity in temporal areas.
In experiment 2, we also found that lowfrequency MEG signals encoded information about the cursor kinematics. Voxels with significant group level activity at negative lags are summarized in figure 4.
In the execution condition, we observed significant amplitude related activity (figure 4(a)) in contralateral sensorimotor areas (peak at lag −0.2 s), and to a weaker extent in occipital, parieto-occipital, parietal and ipsilateral premotor areas. Position (figure 4(b)) and velocity (figure 4(c)) information was consistently expressed in contralateral motor and parietal areas as well as parieto-occipital areas and to a weaker extent in occipital areas. As mentioned before, due to strong cross-correlations between positions and velocities across lags (supplementary figures 14(e), (f)), we observed similar effects in both kinematic variables at different lags. For velocity, the contralateral activity peaked at the lag −0.2 s (figure 4(c)) and consequently for position at lag 0.3 s (supplementary figure 4).
In the observation condition, the kinematics related activity (figures 4(d)-(f)) was generally smaller compared to the execution condition (figures 4(a)-(c)). Amplitude information (figure 4(d)) was expressed strongest in parietooccipital areas, followed by ipsi-and contralateral sensorimotor areas. Position (figure 4(e)) and velocity (figure 4(f)) related activity was strongest in PP areas.  Figures 1(d),(e) displays recorded and decoded test set trajectories for representative participants. Qualitatively, the UKF decoded trajectories (blue) matched the cursor trajectories (gray) best in terms of time course and amplitude range in the execution condition. The decoded trajectories of the other algorithms were either too small (red, WF) or too large (green, KF). In the observation condition, the congruence between the decoded and recorded trajectories was lower than in the execution condition.

Cursor kinematics decoding
The quantitative participant-level results are listed in supplementary tables 1-6. All three decoding algorithms directly decoded the directional kinematics. To quantify the goodness of fit for the nondirectional kinematics, we converted the decoded position and velocity trajectories into distance and speed trajectories and compared them with the recorded ones. The correlation and SNR between the recorded and decoded directional kinematics were significantly different from shuffled data for all participants in the execution condition. There were some combinations of algorithm, kinematic signal and condition that were not significant for specific participants. For two participants in experiment 1 the SNR of all algorithms and kinematic signals was not significantly different from chance in observation condition. The WF and KF correlations for distance and speed, inferred from the decoded positions and velocities, were not significantly different from chance for the majority of participants in experiment 1 (execution: 11 participants, observation: 13) and observation condition in experiment 2 (13 participants). Whereas the UKF correlations were not significant for the minority in experiment 1 (execution: 1 participant, observation: 5) and observation condition in experiment 2 (five participants).
The group level results are summarized in figure 5. The figure shows the performance metrics for each algorithm, condition and experiment. To reduce the number of plots, we grouped the results into directional (2D position, 2D velocity) and nondirectional (distance, speed) kinematics. The boxplots in figures 5(a),(d) summarize test set correlations between the recorded and decoded kinematics. The UKF achieved the highest test set correlations across experiments, conditions and kinematics. With regard to the directional kinematics, the UKF correlations in execution and observation condition in experiment 1 were 0.40 ± 0.06 SD and 0.31 ± 0.08 SD ( figure 5(a)). The WF correlations were significantly lower in execution (mean difference: 0.08 ± 0.02 SD) and observation (0.06 ± 0.02 SD) condition compared to the correlations reached by the other algorithms. The results were similar in experiment 2 ( figure 5(d)); the UKF correlations in execution (0.56 ± 0.05 SD) and observation (0.40 ± 0.08 SD) condition were significantly higher than the WF correlations. The difference in correlations between the UKF and KF did not reach significance in either experiment and condition. The difference in UKF correlations between experiments was significant (Wilcoxon rank sum tests, execution condition mean difference = 0.16, p-value = 2.1 * 10 -7 ; observation condition mean difference = 0.09, p-value = 3.4 * 10 -3 ). With regard to the nondirectional kinematics, we observed a large gap (mean difference: ≥ 0.1) between the UKF and the other algorithms in both experiments and conditions (figures 5(a),(d)). Recall that the UKF was the only algorithm that incorporated distance and speed information, encoded in the M/EEG signals, into the predictions of position and velocity. Figures 5(b), (e) summarize the performance of the algorithms in terms of SNR. The SNR of the UKF was always higher than the SNR of the KF. For the directional kinematics, we observed the highest SNR for the WF. It was on average 0.38 dB (0.16 dB observation) in experiment 1 ( figure 5(b)) and 1.08 dB (0.42 dB observation) in experiment 2 (figure 5(e)); it was significantly higher than the SNR of the KF in both experiments and conditions. The difference in SNR between the UKF and WF was significant in experiment 1 ( figure 5(b)), and in the observation condition of experiment 2 ( figure 5(e)), while the UKF and WF were on a par in the execution condition. Regarding the nondirectional kinematics, we observed a large gap (mean difference: ≥1 dB) between the UKF and the other algorithms. The differences were significant across experiments and conditions. While there was no significant difference in SNR between the KF and WF in experiment 1 ( figure 5(b)), the KF achieved a significantly higher SNR (mean difference: approx. 0.5 dB) in experiment 2 ( figure 5(e)).
The DSSR results are summarized in figures 5(c), (f). The differences in DSSR between the algorithms were significant in call cases (kinematic groups, conditions, experiments). Compared to the KF and UKF, the WF DSSR was lower by at least 5 dB. This mismatch in amplitude is clearly visible in figures 1(d), (f). The UKF DSSR was closest to the target value of 0 dB in experiment 1, while in experiment 2 the KF DSSR was closest to 0 dB. The position, velocity, distance and speed specific group level results are displayed in figure 6 (correlation and SNR) and supplementary figure 11 (DSSR). In experiment 1, the algorithm-specific correlation, SNR and DSSR were generally similar across the directional kinematics (figures 6(a), (b)), supplementary figures 11(a)). In experiment 2, the correlation and SNR were higher for the vertical component (vertical position, vertical velocity) than for the horizontal component (figures 6(c), (d)). We additionally observed that the UKF and KF correlations were higher for the positions than for the velocities ( figure 6(c)). The results for distance and speed were consistently lower than the ones for position and velocity ( figure 6). Compared to distance, the correlations for speed were similar (figures 6(a), (c)) while the SNR was higher (figures 6(b), (d)).
Next, we evaluated the effect of including only distance or speed or neither in the UKF model. In figure 7, we show the results in terms of correlation and SNR. If we did not consider speed or distance in the UKF model, the associated test set correlation and SNR declined. For example, if we excluded distance in equation (11), the correlation and SNR for the test set distance trajectories dropped and also to a smaller extent for the position trajectories. If neither distance nor speed was included, the UKF model was equivalent to the linear KF model and achieved a similar accuracy as the KF. Taken together, using distance and speed in (11) resulted in the highest decoder accuracy for the UKF.

Discussion
We showed that the low-frequency EEG and MEG activity encoded rich kinematic information during the PTT. At the group level, posterior-parietal, parieto-occipital (both conditions) and contralateral premotor and primary sensorimotor (execution condition) areas encoded information about directional (position and velocity) and nondirectional (distance and speed) kinematics. Based on the low-frequency M/EEG activity of the current and preceding timepoints, the three decoding algorithms could reconstruct the cursor kinematics at the current time-point with moderate congruence. Across experiments and directional kinematics, the correlations of the UKF were 0.49 (0.10 SD) in the execution and 0.36 (0.09 SD) in the observation condition. They were on average 0.08 (0.02 SD) higher than the correlations of the WF . A linear state space model (KF) was sufficient to significantly improve the correlations upon the WF (∆ correlation = 0.06 ± 0.01 SD) at the cost of significantly lower SNR (∆ SNR = −1.5 dB ± 0.6 dB SD). Using information about nondirectional kinematics, encoded within the low-frequency M/EEG activity, the UKF could significantly improve the SNR compared to the KF (∆ SNR = 0.6 dB ± 0.1 dB SD) and at the same time reduce the amplitude mismatch between recorded and decoded trajectories to a minimum (DSSR = −0.3 dB ± 0.5 dB SD).
Despite we observed significantly lower correlations for the WF compared to the state space models (KF, UKF), its SNR was significantly higher for the directional kinematics. The optimization criterion of the WF was to minimize the MSE which is equivalent to maximizing the SNR defined in (14). Whereas the state space models did not directly maximize the SNR; they learned the distribution of the kinematics in the train set (e.g., the amplitude range), and thereby also optimized the DSSR to be close to 0 dB. Intuitively, correlation, SNR and DSSR are related. Since the WF maximized the SNR of the directional kinematics, it achieved the highest SNR at the cost of poor DSSR, while the state space algorithms achieved a DSSR close to 0 dB at the cost of lower SNR for the directional kinematics ( figure 5). If the correlations were high, as in experiment 2 (execution condition), the UKF could maintain the same SNR as the WF and at the same time a DSSR close to 0 dB (figures 5(d)-(f)). Knowing that the KF and UKF essentially differed in the sense that the UKF utilized the nondirectional Figure 5. Group level cursor kinematics decoding results for each experiment, condition and algorithm. The kinematics were summarized by grouping them into directional (horizontal position, horizontal velocity, vertical position, vertical velocity) and nondirectional (distance, speed). (a) Experiment 1 (EEG, 15 participants). Boxplots showing the correlation between the recorded and decoded kinematics for execution (left) and observation (right) conditions. The algorithms are color-coded. Each dot summarizes the average CV test set correlation of a participant. Significant differences between the algorithms are highlighted (two-sided, permutation, paired t-tests, df = 14, 10 000 permutations, sig. level = 0.05, FDR adjusted p-values for n metrics * n algorithms * n conditions * ngroups = 36 tests). (b) As in (a) for the signal to noise ratio. (c) As in (a) for the decoded signal to signal ratio. All results were significant (p < 0.001). (d)-(f) As in (a)-(c) for experiment 2 (MEG, 19 participants). Significant differences in (d) and (e) are highlighted (two-sided, permutation, paired t-tests, df = 18, 10 000 permutations, sig. level = 0.05, FDR adjusted p-values for n metrics * n algorithms * n conditions * ngroups = 36 tests). All results in f were significant (p < 0.001).
kinematic information, and that we observed a significantly higher SNR for the UKF across all kinematics, it follows that utilizing the nondirectional information improved the decoder accuracy. In summary, the UKF achieved the highest qualitative (figures 1(d), (e)) and quantitative in terms of correlation and DSSR (figures 5(a), (c), (d), (f)) decoding accuracy among the three algorithms. In terms of SNR, the accuracy of the UKF was on a par with the WF in one case (experiment 2, execution condition) and lower in the other three cases (figures 5(b), (d)).
The UKF was the only algorithm that used distance and speed information to improve the position and velocity estimates. Using either speed or distance in equation (11), had a negligible effect on the decoder accuracy of the directional kinematics and a negative effect on the nondirectional kinematics (figure 7). If neither distance nor speed were included, the accuracy of the UKF declined to the level of the KF. Yeom et al explored a different way to incorporate speed information into the velocity estimates [66]. In an MEG study, they first decoded velocity and speed during center-out movements, then normalized the decoded velocity and multiplied the result with the decoded speed estimate. Compared to linear decoders, their two step approach could slightly, yet significantly improve the accuracy. We also varied the model order to identify how the number of previous lags affected the decoder accuracy (supplementary figure 13). Compared to the third order model (figure 5), the correlations and SNR were marginally smaller for a first order model (∆ correlation = −0.01 ± 0.01 SD; ∆ SNR = −0.1 dB ± 0.1 dB SD) and indifferent for a fifth order model (∆ correlation = 0.0 ± 0.01 SD; ∆ SNR = 0.0 dB ± 0.1 dB SD). As reported in [20], the WF accuracy increased for higher orders (supplementary figure 13). Despite the increase, the correlations (supplementary figures 13(a), (d)) for a fifth order WF were still lower compared to a fifth order UKF (∆ correlation = 0.05 ± 0.02 SD for the directional kinematics). Taken together, the UKF still outperformed the other algorithms in terms of correlation as the model order increased; higher order models did not improve the accuracy upon a third order UKF.
Since the spectral filter characteristics have a large impact on the performance of lowfrequency decoders [66], we tested three low-pass filter cut-off frequencies 0. For the nondirectional kinematics, the correlations increased, since frequencies above 0.8 Hz had a considerable contribution to the power spectral density of distance and speed (supplementary figure 15). Increasing the cut-off frequency had a positive effect on the UKF SNR (0.5 dB directional, 1 dB nondirectional) and DSSR.
Previous works in the context of low-frequency M/EEG kinematics decoding have mainly studied center-out tasks [19,27,[67][68][69][70]. Since correlation and SNR are sensitive to various parameters (e.g., preprocessing, statistical properties of the signals being compared), a direct comparison across studies is not straightforward [25]. Bradberry et al were the first to show that velocity trajectories can be decoded from low-frequency MEG [67] and EEG [19] activity during center-out reaching tasks. They decoded 2D velocity trajectories from MEG during a visuomotor adaptation experiment in [67]. Their preexposure condition was closest to our execution condition. They obtained average (five participants, two dimensions) correlations of 0.4. In the EEG study [19], they studied 3D self-paced center-out movements and decoded velocities with average (five participants, three dimensions) correlations of 0.3. Kim et al decoded continuous, repetitive 3D arm movement trajectories from low-frequency EEG activity with multiple linear regression (MLR) and kernel ridge regression decoders [23]. After correcting for eye artifacts, they obtained grand average (ten participants, three dimensions) correlations of 0.5 (0.38 MLR) in executed movements and 0.4 (0.3 MLR) in imagined (and observed) movements. Their correlations were in a similar range, compared to our findings in the PTT in experiment 2 (figures 5(d)-(f)). In an even simpler task (rhythmic 2D circular arm movements), Georgopoulos et al reconstructed arm movement trajectories from low-frequency MEG activity with high congruence (correlation = 0.85) [15]. Recently we reported slightly lower correlations (correlation = 0.68) for low-frequency EEG activity [29].
In the invasive domain, executed arm movement kinematics have been decoded from ECoG activity of humans [24,28,[71][72][73]. Schalk et al reported that the low-frequency ECoG activity contributed the most to the decoding of trajectories in a simple, repetitive 2D PTT [24]. Using low-frequency and spectral features, they predicted positions and velocities with correlations of 0.49 (0.10 SD). Using similar features, Pistohl et al predicted cursor trajectories with a linear KF (correlation = 0.33 for six participants; correlation = 0.43 for three participants with ECoG grids covering hand/arm areas) in a more complex 2D continuous movement task [71]. Nakanishi et al studied a 3D continuous reaching and grasping task in three participants [28]. They decoded wrist and elbow positions and elbow and shoulder joint angles. The average correlation for the 3D wrist position was 0.51 and the average normalized root mean squared error (nRMSE) was 0.29. In a follow up study with three additional participants, they obtained average correlations of 0.62 and a root mean squared error (RMSE) of 55 mm [74]. That is, the M/EEG results in terms of correlation reported here are within the range of those using implanted ECoG grids during continuous, goal-oriented movement tasks, whereas the low nRMSE and RMSE confirm the higher SNR of ECoG compared to M/EEG [75] Expression of executed arm movement kinematics in neural spiking activity has been studied in nonhuman primates (NHPs) [8,31,36,76,77]. An early work by Paninski et al studied neural spiking activity during a comparable 2D PTT in three NHPs [77]. They decoded hand positions from the spike rates of few M1 neurons (5 to 19 neurons) with a linear regression model and reported an SNR of 1.3 dB (0.8 dB SD). Our execution condition results for the UKF (−0.8 dB experiment 1; 0.9 dB experiment 2) and WF (0.4 dB experiment 1; 1.1 dB experiment 2) were lower. Later, Li et al used neural spiking activity from M1, S1, PMd and PP neurons of two NHPs to decode cursor kinematics during a simple 2D PTT and a center-out task [36]. They applied a similar UKF that used distance and speed information to decode positions and velocities. In the PTT, they reported that a first order UKF decoded positions (correlation = 0.85, SNR = 5.1 dB) and velocities (correlation = 0.52, SNR = 1.5 dB) with high congruence. Our execution condition correlations for velocity were in a similar range, while the correlations for position and the SNRs were lower ( figure 6). The discrepancy in correlations could be explained by the task difficulty. Their 2D trajectories were repetitive, while ours were pseudo-random. They also reported results for a WF. In their study, the UKF improved the correlations upon the WF for the positions (∆ correlation = 0.09) and velocities (∆ correlation = 0.06). We also observed significantly higher correlations for the UKF compared to the WF (∆ correlation = 0.06).

Recently Makin et al compared the UKF proposed by
Li et al to a recurrent exponential family harmonium (rEFH) filter during a consecutive 2D reaching task to multiple targets [78]. The rEFH filter achieved an SNR of 6 dB compared to 4 dB for the UKF. In a study with two NHPs, Stavitsky et al compared the decoder accuracy between population and spiking activity with a linear KF [79]. In a center-out reaching task, they reported similar correlations for the lowfrequency population (0.78) and spiking (0.77) activity. Compared to M/EEG, the fine spatial resolution of intracortical electrode arrays enables the reconstruction of movement trajectories with a high accuracy in terms of correlation and SNR.
In both experiments, we studied the PTT in two experimental conditions. Compared to the execution condition, we observed a drop in performance in the observation condition ( figure 5). For the UKF the average difference in correlation and SNR was 0.13 (0.08 SD) and 1.1 dB (0.6 dB SD). Kim et al also reported lower correlations (∆ correlation = 0.1), if the movements were observed and imagined [23]. Korik et al used spectral features to decode repetitive 3D reaching movements from EEG activity [68]. They reported a stronger drop in correlations (∆ correlation = 0.2), if the reaching movements were imagined rather than executed. Using low-frequency EEG activity, Ofner et al reported correlations of 0.32 (0.22 SD) for rhythmic movement imaginations [22]. They were substantially lower than the correlations (0.70 ± 0.12 SD) reported earlier in executed 3D movements [80]. The drop in decoding performance could be explained by the brain encoding less kinematic information. A body of work has investigated how cortical areas, involved in movement control, modulate to different volitional states, including executed, imagined, attempted or observed movements [81][82][83][84][85]. In supplementary figures 3 and 6 we compared the group level kinematics related voxel activations between conditions. In both experiments, voxels in primary sensorimotor cortex (SM1), supplementary motor cortex, PMd and PP encoded significantly more kinematic information in execution condition. Executed movements have been repeatedly shown to activate motor and sensory areas stronger than observed movements [81,82,84].
The activity in SM1, PMd, PP and PO simultaneously encodes information about various movement parameters in executed and observed movements [8, 20, 33-35, 73, 86-88]. Our group level encoding results (figures 3 and 4) indicate that directional and nondirectional kinematics were expressed differently across cortical areas. We observed distinct peaks in SM1 for the nondirectional kinematics in both experiments (figures 3(a) and 4(a)), supplementary figures 1,4,7,9). The SM1 activity peaked between lags −0.1 s and 0.0 s in experiment 1 and lags −0.3 s to −0.2 s in experiment 2. Considering the delays between hand/finger and cursor movements (55 ms in experiment 1; 190 ms in experiment 2), the SM1 peak was phase-locked to the executed hand/finger movements. Tuning effects of SM1 to nondirectional kinematics have been reported in ECoG [34], MEG [16,35,89] and EEG [18]. We observed that the directional kinematics (position and velocity) were consistently encoded in PO (figures 3(b), (c), (e), (f)) and PP (figures 4(b),(c),(e),(f)) in both conditions. Voxels in contralateral SM1 and PMd also consistently encoded directional information in the execution condition of experiment 2 (figures 4(b),(c)) but not in experiment 1 (figures 3(b),(c)). In figures 3 and 4 we report consistent group level effects. That is, for an effect to be significant, the dipole orientations have to be similar across participants. If only the strength of the effect is considered (= dipole orientations are allowed to vary across participants), SM1, PMd and PP voxels encoded significantly more directional information in execution condition in experiment 1 (supplementary figure 3). This means that the directional information was consistently expressed in PO and PP areas at the group level, and participant specific in PMd and SM1. The less consistent effect in SM1 during experiment 1 could to some degree also be attributed to a lower spatial localization accuracy of EEG source imaging for template head models compared to MEG source imaging [90].
In spite of similar task dynamics (supplementary figure 14) and visual stimuli ( figure 1(a)) in both experiments, the movement tasks differed in execution condition; we studied arm movements in experiment 1 and index finger movements in experiment 2 ( figure 1(b)). The two movement types activate fronto-parietal networks for reaching and grasping [91]. Among the fronto-parietal networks, reaching movements mainly activate the dorsal reach system and grasping movements mainly the lateral grasp system [87,[92][93][94]. As expected, the kinematic information was encoded strongest in medial areas along the dorsal stream in experiment 1 (figures 3(a)-(c), supplementary figures 1, 3 and 7) and more lateral areas in experiment 2 (figures 4(a)-(c), supplementary figures 4, 6 and 9). Since the fingers and hand have a stronger representation in SM1 than the elbow and shoulder [95], more kinematic information could have been expressed in experiment 2. This could explain the higher decoding accuracy (∆ correlation(UKF) = 0.16 ; ∆ SNR(UKF) = 1.7 dB) in execution condition in experiment 2 (figure 5). In the observation condition the hand/finger was resting. Nonetheless, the decoder accuracy was higher in experiment 2 (∆ correlation(UKF) = 0.09 ; ∆ SNR(UKF) = 1.0 dB). This difference could be attributed to the recording modality (EEG in experiment 1; MEG in experiment 2) or the spatial coverage (64 channels in experiment 1; 129 channels in experiment 2).
Despite the promising results, our study suffers from a number of limitations that need to be addressed in future research. The most critical one is that we analyzed the data offline with non-causal filters. Online, the spectral filters would introduce large delays in the order of a second, rendering lowlatency control infeasible. If the filter characteristics were relaxed, the decoding accuracy of the WF would drop substantially (supplementary figure 12). Fortunately, the results indicate that the accuracy of the KF and UKF would be marginally affected. In an exploratory closed loop BCI study, we recently demonstrated that a KF decoded arm movements online with low-latency so that healthy individuals could track a target with a robotic arm during a PTT [30]. In [30] we used the KF to smooth the estimates of a WF. As a result we observed a similar amplitude mismatch as for the WF in this study (figures 5(c), (f)); the amplitudes of the decoded position trajectories were three to four times smaller than the ones of the actual movements, rendering feedback training impracticable. We believe that at the start of feedback training, the BCI users prefer larger errors over not being able to produce movements in the same range as during manual control. A second limitation, concerns the transfer to paralyzed individuals. In this study, we decoded executed and observed movements. As the participants were continuously moving in the execution condition, the decoders had access to a mixture of efferent (movement intention) and reafferent (feedback) activity. The reafferent activity is either completely or partially missing in paralyzed individuals. Nonetheless, the observation condition decoding accuracy should transfer to paralyzed individuals. Moreover, we think that if a movement attempt strategy is used rather than mere observation, a similar performance to execution condition should be feasible [84,85]. Using intracortical BCIs, paralyzed individuals achieved similar decoding accuracy in attempted movements as NHPs in executed movements [6,96,97]. Even if the offline results can be transferred to online experiments, the moderate decoding accuracy obtained here is not sufficient to accurately control an end-effector. Using motor imagery of distinct movements and band power features, Wolpaw and McFarland demonstrated that 2D cursor control could be significantly improved with feedback training [4]. After training for several weeks, the study participants could control a cursor in a 2D center-out task with average correlations of 0.63. Edelman et al used a similar control strategy and reported that their participants could improve the control skill significantly faster in a 2D PTT compared to a center-out task [98]. In another online study, Bradberry and colleagues reported that they could decode intended movement trajectories from low-frequency activity with high congruence [99]. Within a single-session experiment the participants could control an end-effector in a center-out movement task. However, they did not report whether the task performance was significant, leading to an inconclusive debate [100,101]. Although there is evidence that individuals can learn to voluntarily modulate low-frequency activity with feedback training [102], it remains unclear how fast the kinematics decoding accuracy can be improved.
There were also some potential confounding factors that might have affected the presented results. The first one concerns eye movement artifacts. Allowing eye movements enabled us to study similar visuomotor (and oculomotor) tasks as invasive studies. However, they also give rise to prominent EOG artifacts in low-frequency M/EEG signals [103,104]. Due to the nature of the PTT, they co-varied considerably with the directional kinematics. We used a state-ofthe-art correction algorithm to attenuate these artifacts and maintain brain activity [47]. Although we cannot rule out that residual eye artifacts contributed to the decoding results, the encoding results indicate that the brain activity dominated. If residual eye artifacts had dominated, we would have expected strong activity in prefrontal and anterior temporal cortex (cortical areas closest to the eyes). Across conditions and experiments, the kinematics related activity at prefrontal and anterior temporal areas was much lower compared to cortical areas along the dorsal stream (figures 3 and 4; supplementary figures 1, 2, 4, 5, 7-10). Movement artifacts could have also confounded the execution condition decoding results. Significant activity at the inferior temporal lobe might indicate the presence of an arm movement-related artifact that was phase-locked to the nondirectional kinematics in experiment 1 ( figure 3(a)). Alternatively, this activity could have originated in deeper brain structures (e.g. thalamus, cerebellum) which we did not consider in the source imaging model. As in [24,36], we decoded cursor trajectories and ignored potential effects related to the target stimulus. This raises the question to which extent target related information was utilized by the decoders. During the PTT, the target and cursor trajectories were highly correlated (figures 1(d),(e)), rendering a disentanglement of cursor and target related effects impracticable. Studying multiple tasks with different dynamics could help to disentangle cursor (arm/finger) and target related effects [105].

Conclusion
Previous work suggested that low-frequency EEG and MEG signals encode information about directional (position and velocity) and nondirectional (distance and speed) arm movement kinematics. We found direct evidence that directional and nondirectional kinematic information is simultaneously detectable in low-frequency M/EEG signals. Moreover, movement trajectories could be reconstructed with significantly higher accuracy in terms of correlation and DSSR using an UKF that explicitly models the nonlinear relation between directional and nondirectional kinematics than by using linear Kalman (KF) and WFs which do not combine both types of kinematics. In terms of SNR the UKF outperformed the linear KF in all cases, while it was outperformed by the WF in three out of four cases. Further research is necessary to improve the decoder accuracy, particularly the strengths of the UKF in terms of correlation and DSSR and the WF in terms of SNR should be combined. Apart from improving the algorithms, it is essential to assess how the results transfer to attempted movements in paralyzed individuals and to which extent the accuracy can be improved with feedback training.