Decoding three-dimensional reaching movements using electrocorticographic signals in humans

Initially, ECoG signals were visually inspected in both the spectral and temporal domain. Any channels that displayed non-physiologic activity or epileptic activity were excluded from all further analyses (Patient 1:7 electrodes, Patient 2:7 electrodes, Patient 3:10 electrodes, Patient 4:4 electrodes, Patient 5:19 electrodes). Remaining channels were included in all further analyses to avoid biasing the results by preselecting channels. Next, the spectral and temporal domain signals for each trial were examined. Trials containing non-physiologic spikes or interictal epileptic activity were rejected. The hand position data collected for each trial were also visually examined and trials in which the hand position left the sampling range of the Flock of Birds receiver and trials in which patients initiated a reach to an incorrect target were also rejected. After this screening process, the total number of trials analyzed was: 245, 76, 104, 221, and 202 for Patients 1–5 respectively. ECoG signals were then re-referenced to the common average. For 8 × 8 channel ECoG grids spanning multiple amplifiers, individual common averages were calculated for each recording amplifier. For 1 × 4, 1 × 6, and 1 × 8 channel strips, individual common averages were calculated for each strip of electrodes. Signals were band-pass filtered between 0.1 and 260 Hz using a 3rd order butterworth filter. Additionally, a 3rd order butterworth filter with a 5 Hz bandwidth was used to remove all noise harmonics below 260 Hz, including harmonics of the 60 Hz line frequency as well as harmonics of a 100 Hz artifact observed in the recording system used. Filters were run forwards and backwards to avoid inserting phase distortions into the signals.

Next, ECoG signals were segmented into 300 ms time windows with shifts of 50 ms between windows. Spectral power was estimated in 2 Hz bins with center frequencies from 3 to 253 Hz using the autoregressive maximum entropy method with a model order of 75 (Marple 1987). Although the 300 ms time window may limit the accuracy of power estimates for the lowest frequencies, it was chosen to balance the accuracy of spectral power estimates and temporal resolution. As ECoG power spectra are not normally distributed, power spectra were normalized using a log transform. Power spectra were then converted to z-score values, using the mean and standard deviation of spectral power from 200 ms after the beginning of the hold-A period until the end of the hold-A period. The z-score operation accounts for the 1/f fall-off in spectral power by ensuring that the mean and variance of spectral power is the same at each frequency. Additionally, positive and negative z-scores indicate respective increases and decreases in spectral power relative to the hold-A period. The hold-A period was chosen as the baseline period for the task as patients had not received any target information, were not moving, and were maintaining their hand in a similar position in all trials. Spectral power was then averaged into seven canonical frequency bands: theta (4–8 Hz), mu (8–12 Hz), beta 1 (12–24 Hz), beta 2 (24–34 Hz), gamma 1 (34–55 Hz), gamma 2 (65–95 Hz), and gamma 3 (130–175 Hz). The beta band was separated into two components because two distinct components have been described (Pfurtscheller et al 1997). These seven frequency bands were chosen to ensure inclusion of relevant frequency bands while avoiding all noise harmonics of the 60 Hz line noise and a 100 Hz harmonic observed in the recording system. Because the frequency bands used have varying bandwidths, their variances differed after averaging into frequency bands; therefore, the band-averaged spectral power was converted to a z-score. The two calculations were used to ensure that each frequency bin contributed equally to the band-averaged power estimates, and that the variance was similar across frequency bands, irrespective of the bandwidth.

In addition to the spectral features described above, the LMP, which is obtained by low-pass filtering time domain signals, has been shown to contain movement-related information (Schalk et al 2007, Pistohl et al 2008, Hotson et al 2014). To calculate the LMP, signals were segmented into 300 ms time windows with shifts of 50 ms. A 2nd order Savitzky–Golay filter with a −3 dB cutoff frequency at 3.6 Hz was run for each 300 ms window. Finally, the LMP signal was z-scored relative to the hold-A period. Therefore, the variance of the LMP was equalized to the variance of the spectral power in each of the canonical frequencies. All spectral frequency bands and the LMP amplitude were visually inspected to verify that they were normally distributed with equal variances.

Kinematic data was recorded as 3D positions with the positive x-axis oriented towards the patient in the anterior–posterior direction, the positive y-axis oriented laterally to the left, and the positive z-axis oriented inferiorly. The hand position data recorded in 3D space was differentiated to determine velocity in the three cardinal directions. Non-directional hand speed was calculated from the velocity at each time point. The onset of movement for each trial was determined as the time point at which hand speed crossed 20% of the maximum hand speed for the trial. To align kinematic information with ECoG spectral and temporal features, kinematic information was segmented into 300 ms windows with shifts of 50 ms between windows and was averaged within each window.

To examine the dimensionality of the kinematic dataset, we performed principle components analysis (PCA) on hand speed, velocity in the three cardinal directions, and position along the three cardinal axes after normalizing each of the 7 kinematic parameters between 0 and 1. The percent of information explained by each principle component was compared. Additionally, the cross-correlation matrix between each of the seven kinematic parameters was examined for any consistent relationships.

The ultimate goal was to develop a decoding model to predict time courses of speed, velocity, and position. However, the planning delay caused the task to contain both rest and movement periods. As previous BCI studies have found that using a hierarchical decoding model can improve decoding accuracy (Pfurtscheller et al 2010, Flamary and Rakotomamonjy 2012, Aggarwal et al 2013), we hypothesized that the relationship between neural activity and kinematics would differ between the movement and rest periods, leading to increases in kinematic decoding accuracy through the use of a hierarchical regression model. We implemented a two-step hierarchical decoding model to decode kinematics of reaching movements as shown in figure 2.

Zoom In Zoom Out Reset image size

The first component classified whether the patient was moving or not with a logistic regression model using elastic net regularization. A threshold of 10% of the maximum movement speed was applied to generate binary movement and rest labels. The features used to train the model consisted of the z-scores of each of the eight features (seven frequency bands and LMP). To account for potential differences in the optimal time lag between kinematics and neural activity for different channels and features, the Pearson correlation coefficient between training set hand speed and neural activity was calculated with the neural activity leading hand speed at time lags between 0 and 1 s in 50 ms steps. For each channel and feature, the training set data was used to select the lag producing the maximum absolute correlation coefficient with movement speed. To train the logistic regression model, we minimized the loss function shown in equation (1), where X is an n sample ×d feature input matrix of ECoG data, w is a d x 1 weight vector, y is an n x 1 vector of class labels (1 movement, −1 rest), and the ${\lambda }_{1}$ and ${\lambda }_{2}$ are the hyperparameter weights associated with the l1 and l2 norms respectively

$$\begin{eqnarray}&&L(w)=\displaystyle \frac{1}{n}\displaystyle \sum _{n}\mathrm{log}(1+{{\rm{e}}}^{X(n)*w*y})+{\lambda }_{1}{\Vert w\Vert }_{1}+{\lambda }_{2}{\Vert w\Vert }_{2}.\end{eqnarray} \tag{ 1 }$$

After finding the set of weights that minimize the loss function above, the output of the model was calculated as shown in equation (2)

$$\begin{eqnarray}&&p(y==1| X)=\displaystyle \frac{1}{{{\rm{e}}}^{-X*w}}.\end{eqnarray} \tag{ 2 }$$

The output of the logistic regression model can be considered as a probability that the current time window is from a movement period. Time windows with output probabilities above 0.5 were classified as movement periods and time windows with output probabilities below 0.5 were classified as rest periods. Within the training set, seven fold cross-validation was used to determine the optimal regularization hyperparameters. Finally, the entire training set was used to train the model and the testing set was used to evaluate the accuracy of the model.

The output from this classification step was used to switch between two regression models, one to predict each kinematic parameter while patients were at rest and a second model to predict each kinematic parameter during movement periods. We used a similar method to a prior study that demonstrated the ability to use a PLS regression model to decode continuous traces of kinematic data from ECoG in non-human primates (Chao et al 2010). The PLS model estimates a lower dimensional latent structure and uses this latent structure to fit a regression in order to avoid over-fitting (Wold et al 1984). For our model, the inputs consisted of the z-scored activity for each of the seven canonical frequency bands and the LMP at each channel. Individual feature vectors were produced for time lags between −500 and 1000 ms in 50 ms steps, with positive lags indicating neural activity leading kinematics and negative lags indicating neural activity lagging after kinematics. As the filters used during the preprocessing steps were run both forwards and backwards to avoid phase distortion, the neural activity in each time-window was not necessarily causal relative to the kinematic data, therefore we chose to use neural activity both leading and lagging the kinematic time courses to be predicted. Therefore, the results should be interpreted to represent the entire sensorimotor system and not necessarily causal activity related to motor planning alone. The outputs for the model consisted of movement speed, 3D movement velocity, and 3D hand position. The relationship between ECoG activity and kinematics is described by equation (3), with M(t) representing a kinematic parameter and w representing a prediction weight for a specific channel, feature, and time lag.

$$\begin{eqnarray}M(t) & = & \displaystyle \sum _{{\rm{Channels}}}\displaystyle \sum _{{\rm{Features}}}\displaystyle \sum _{\tau =-500\;{\rm{ms}}}^{1000\;{\rm{ms}}}w({\rm{ch}},{\rm{feat}},\tau )\\ & & \times Z{\rm{score}}({\rm{ch}},{\rm{feat}},t-\tau ).\end{eqnarray} \tag{ 3 }$$

Within the training set, we trained models with 1–24 latent features and used seven-fold cross validation to determine the optimal number of latent features that minimized the mean-squared prediction error. The final prediction model was generated from the full training set and the testing set was used to evaluate the accuracy.

To evaluate the accuracy of movement prediction, we generated 100 randomly sampled training and testing sets. The accuracy of the classification of movement from rest and of the PLS regression models used to predict kinematics were evaluated independently prior to evaluating the full hierarchical model combining the two models. Accuracy of the classification between movement and rest was calculated as the percent of the total test set windows that were correctly classified. Accuracy of the PLS regression model was evaluated by using the actual movement and rest labels to switch between the two PLS regression models and computing the Pearson correlation coefficient and root mean squared error (rMSE) between the actual and predicted values of each of the 7 kinematic parameters (speed, velocity: V_x, V_y, V_z, and position: X, Y, Z). Finally, the accuracy of the full model was quantified by calculating the Pearson correlation coefficient and rMSE between actual values and predicted values for each of the kinematic parameters.

We evaluated the statistical significance of the models using two surrogate methods to ensure that accurate predictions were generated from the specific patterns of ECoG activity and not from any systematic bias. First, we generated a surrogate kinematic dataset by randomly reordering the trials within the training set, randomly selecting a new trial onset from within each training set trial, and generating a new time course by wrapping data from the beginning of the trial to the end of the trial. This procedure ensured that the temporal relationship between the ECoG signals and kinematics was random, while maintaining the autocorrelation structure of the kinematics. The binary movement labels and continuous kinematic traces were reshuffled in the same fashion, so any differences between actual and surrogate predictions were not driven by the ability to decode movement and rest periods. For each training and testing set, we trained one model using the original kinematic data and a second model using the surrogate kinematics. Both models were tested using the original testing set. A second surrogate method that was similar to one used in other studies of kinematic prediction (Chao et al 2010, Hotson et al 2014) was also used to randomly vary ECoG features and channel assignments. For each training set, we reshuffled the channel and frequency assignments of the prediction weights 100 times and generated test accuracies using both the original and reshuffled weights. The statistical significance of both the movement classification and kinematic prediction models was evaluated using a Wilcoxon rank sum test to compare the median actual accuracy with the median surrogate accuracy. Bonferroni correction was used to correct for the total number of predictions (five patients × seven kinematic variables × two surrogate methods).

To determine the optimal method for generating movement trajectories, we compared the accuracy of the original position predictions with position predictions generated by concatenating the predicted velocity vectors. The overall accuracy was then quantified by determining the percent of trials in which the actual and predicted trajectories ended in the same quadrant.

Because the z-score operation equalized the variance of each feature prior to generating prediction models, the weights for each channel and feature could be compared to evaluate the importance of each feature type and electrode location. While the PLS model weights represent the optimal combination of ECoG signals to decode kinematics, they do not necessarily represent the strength of the relationship between movement kinematics and ECoG data (Haufe et al 2014). Therefore, the decoding model weights were converted to activation patterns describing the encoding of kinematics by ECoG signals (equation (4)) (Haufe et al 2014), where W represents the matrix of prediction weights, Σ_X represents the ECoG covariance matrix, Σ_S represents the kinematic data covariance, and A represents the activation patterns.

$$\begin{eqnarray}&&{\boldsymbol{A}}={{\rm{\Sigma }}}_{X}{\boldsymbol{W}}{{{\rm{\Sigma }}}_{S}}^{-1}.\end{eqnarray} \tag{ 4 }$$

To examine the importance of features in the classification of movement and rest conditions, logistic regression prediction weights were averaged across each of the 100 training sets. For each frequency band, the magnitude and location of average model weights was mapped onto an atlas brain using a weighted spherical Gaussian kernel (full-width at half-max: 9.2 mm) centered at each electrode location. In patients with overlapping electrode coverage, Gaussian kernels were combined using a weighted average based upon the distance from each electrode. For the PLS regression model predicting movement period kinematics, the relative importance the absolute value of the prediction weights and activation patterns for a given channel, feature, and lag were normalized by the sum of the absolute value of all prediction weights or activation patterns as shown in equation (5)

$$\begin{eqnarray}&&{\rm{Wnorm}}({\rm{ch}},{\rm{feat}},{\rm{lag}})\\ &&\quad =\displaystyle \frac{|w({\rm{ch}},{\rm{feat}},{\rm{lag}})|}{\displaystyle {\sum }_{{\rm{Features}}}\displaystyle {\sum }_{{\rm{Lags}}}\displaystyle {\sum }_{{\rm{Chans}}}|w({\rm{ch}},{\rm{feat}},{\rm{lag}})|}.\end{eqnarray} \tag{ 5 }$$

Normalized weights and activation patterns were calculated for each channel by summing the across features and time lags (equation (6)). Normalized channel weights and activation patterns were mapped onto a single atlas brain to evaluate the importance of cortical locations within the model

$$\begin{eqnarray}&&\bar{{\rm{Wnorm}}}({\rm{ch}})=\displaystyle \sum _{{\rm{Features}}}\displaystyle \sum _{\;{\rm{Lags}}}{\rm{Wnorm}}({\rm{ch}},{\rm{feat}},{\rm{lag}}).\end{eqnarray} \tag{ 6 }$$

The importance of the eight ECoG feature types used in the PLS model was evaluated in two ways. First the maximum PLS weight and activation pattern at each frequency was determined. Second, after identifying the 1% of model weights and activation patterns with the highest absolute values, the percentage of these weights in each frequency was calculated. The two methods identified features that had a large impact on PLS predictions either because of weights with very high values or a large number of moderately large weights.

To determine whether the proposed model can be applied to an asynchronous dataset, an additional model was trained for patient 5. After removing trials with artifacts, all remaining data including movement, rest, and ITI periods were used. Training and testing sets were generated using eight-fold cross validation with each testing set drawn from consecutive trials. The logistic regression and PLS regression models were calculated from the training sets as described in section 2.6.1 and evaluated on the test sets.

To evaluate the accuracy of the PLS regression model alone, actual movement and rest labels were used to switch between the PLS regressions. PLS model predictions of movement speed were significantly more correlated with actual hand speed than either of the surrogate models in all patients as shown in figure 4. Predictions of movement speed also had significantly lower rMSE values than either surrogate model. Because both the real and surrogate models were trained with speed amplitudes that were high during the movement model training data and low during the rest model training data, speed predictions are higher in the movement period for both actual and surrogate predictions. However, as can be observed in figure S2, the shape of speed predictions is significantly worse for both surrogate methods. Ideal test accuracies across patients are shown in figure 4 for velocity and position. The PLS model accuracy for both directional and non-directional kinematic variables are summarized in table 5. Predictions of velocity had significantly greater correlations and significantly lower rMSE values than both surrogate models for 12 of the 15 velocity components (three directions, five patients). Position predictions had significantly greater correlation coefficients and significantly lower rMSE values when compared to both surrogate models in 10 of 15 position components tested. Therefore, components of 3D kinematics were predicted by our hierarchical PLS model with accuracies better than chance in multiple patients. While the ability to predict movement and rest periods well may influence the size of the correlation values, the surrogate movement and rest models also differentiated movement and rest periods as can be seen in figure S2. Additionally, because movements were made in both directions, predicting movement and rest cannot solely account for the positive correlation values that were observed. Significant statistical differences between actual and surrogate model predictions were also found from the movement period alone (table S1). Therefore, the significant differences between actual and surrogate prediction accuracies were not caused by the difference between movement and rest but because the PLS regression models allowed us to predict kinematics with accuracies better than chance.

Zoom In Zoom Out Reset image size

Median correlations between actual and predicted kinematics are shown in table 6 for the full model predictions using both the logistic regression and PLS regression models. Accuracies for the full model prediction are lower than when actual movement labels are used to switch between PLS regression models (table 5). Additionally, velocity and position accuracies are highest along the x (anterior–posterior) axis and lowest along the y (lateral) axis. When using the predicted movement classes to switch between PLS regression models, predictions of speed were significantly more correlated to the actual movement speed and had significantly lower rMSE values than both surrogate methods in four of the five patients. Predictions of velocity components had significantly greater correlations and significantly lower rMSE than both surrogate models for 11 of 15 velocity components and position predictions had significantly greater correlations and significantly lower rMSE than both surrogate models for 9 of 15 position components. Exemplar time courses of actual and predicted kinematics using the full hierarchical PLS regression model are shown for a single patient in figure 5. As can be seen, time courses of predicted kinematics align well with the time courses of actual reaches.

Zoom In Zoom Out Reset image size

Because position predictions made by concatenating predicted velocity vectors were more accurate than the original PLS model position predictions, this method was used to compute trajectories for each trial. The accuracy of the full model was quantified by determining the percentage of trials in which the actual and predicted movement trajectories ended in the same quadrant. The final accuracy of target predictions was 49.0%, 44.5%, 45.1%, 40.0%, and 66.2% for patients 1–5 respectively, each of which is higher than the 12.5% accuracy expected by chance. The predicted trajectories were finally averaged for each of the eight target classes after normalizing the trial time course. Exemplar average movement trajectories are shown for a single patient in figure 6 demonstrating that movements were predicted accurately in each direction.

Zoom In Zoom Out Reset image size