|Home | About | Journals | Submit | Contact Us | Français|
Electrocorticography (ECoG) signals have emerged as a potential control signal for brain-computer interface (BCI) applications due to balancing signal quality and implant invasiveness. While there have been numerous demonstrations in which ECoG signals were used to decode motor movements and to develop BCI systems, the extent of information that can be decoded has been uncertain. Therefore, we sought to determine if ECoG signals could be used to decode kinematics (speed, velocity, and position) of arm movements in 3D space.
To investigate this, we designed a 3D center-out reaching task that was performed by 5 epileptic patients undergoing temporary placement of ECoG arrays. We used the ECoG signals within a hierarchical partial-least squares regression model to perform offline prediction of hand speed, velocity, and position.
The hierarchical partial-least squares regression model enabled us to predict hand speed, velocity, and position during 3D reaching movements from held-out test sets with accuracies above chance in each patient with mean correlation coefficients between 0.31 and 0.80 for speed, 0.27 and 0.54 for velocity, and 0.22 and 0.57 for position. While beta band power changes were the most significant features within the model used to classify movement and rest, the local motor potential and high gamma band power changes, were the most important features in the prediction of kinematic parameters.
We believe that this study represents the first demonstration that truly three-dimensional movements can be predicted from ECoG recordings in human patients. Furthermore, this prediction underscores the potential to develop BCI systems with multiple degrees of freedom in human patients using ECoG.
Recently, brain-computer interfaces (BCIs), systems designed to translate neural activity into a machine command to control an external assistive device, have emerged with the potential to restore function to patients suffering from conditions such as spinal cord injury, ALS, and stroke. While the idea of using neural activity related to motor intentions to control a prosthetic system has existed for several decades (Brindley and Craggs 1972, Craggs 1975), the past several years has seen significant developments in the ability to design and implement BCI systems using a variety of recording modalities. At one end of the spectrum, penetrating multi-electrode arrays have high spatial resolution and have been used to control closed-loop BCI systems in non-human primates (Wessberg, Stambaugh et al. 2000, Serruya, Hatsopoulos et al. 2002, Taylor, Tillery et al. 2002, Carmena, Lebedev et al. 2003, Velliste, Perel et al. 2008) and human patients (Kennedy and Bakay 1998, Hochberg, Serruya et al. 2006, Kim, Simeral et al. 2008, Collinger, Wodlinger et al. 2013, Wodlinger, Downey et al. 2015). However, these systems require invasive surgical procedures with risks of damage to cortex near the implant site (Bjornsson, Oh et al. 2006) and recording signal quality decreases over time due to the brain’s immune response (Williams, Hippensteel et al. 2007). At the opposite end of the spectrum, both normal controls and patient populations have controlled non-invasive BCI systems (Wolpaw, McFarland et al. 1991, Pfurtscheller, Muller et al. 2003, Pfurtscheller, Neuper et al. 2003, Wolpaw and McFarland 2004, McFarland, Sarnacki et al. 2010, Foldes and Taylor 2013). However, the efficacy of these systems has been limited by the decreased spatial resolution (Cooper, Winter et al. 1965, Pfurtscheller and Cooper 1975) and susceptibility to artifacts (Cooper, Winter et al. 1965, Wolpaw, McFarland et al. 2003). Because of these challenges, electrocorticography (ECoG) recordings made from the surface of the brain have emerged as a potential recording modality for BCI system development due to the balance between signal quality and invasiveness. Furthermore, ECoG signals have been used to develop BCI systems in non-human primates (Rouse and Moran 2009, Rouse, Williams et al. 2013), motor-intact human subjects (Leuthardt, Schalk et al. 2004, Wilson, Felton et al. 2006, Felton, Wilson et al. 2007, Schalk, Miller et al. 2008), and in a quadriplegic human patient (Wang, Collinger et al. 2013), and are stable over several months (Chao, Nagasaka et al. 2010).
Many studies demonstrate the potential for BCI systems by using machine-learning algorithms for open-loop prediction of behaviors such as motor movements. While these open-loop studies do not model the effect of learning and adaptation during closed-loop BCI control, they enable investigation of the feasibility of BCI applications using various recording methods and decoding algorithms with limited data sets. Toward this end, a large body of previous work has shown that ECoG signals can be used to decode motor movements. In animals, ECoG signals have been used to decode kinematic trajectories of two-dimensional (2D) and three-dimensional (3D) limb movements (Chao, Nagasaka et al. 2010, Slutzky, Jordan et al. 2011, Flint, Lindberg et al. 2012, Chen, Shin et al. 2013, Marathe and Taylor 2013). In humans, ECoG signals have also been used to decode trajectories of 2D movements (Schalk, Kubanek et al. 2007, Pistohl, Ball et al. 2008, Sanchez, Gunduz et al. 2008). While several studies have used human ECoG recordings to decode kinematics of movements that were not constrained to two dimensions, the movements decoded were biased because speed and movement direction were correlated in at least one dimension (Hotson, Fifer et al. 2012, Nakanishi, Yanagisawa et al. 2013, Hotson, Fifer et al. 2014). Therefore, the extent of kinematic information that can be decoded from human ECoG recordings is uncertain.
This study utilized a reaching task requiring movements containing 3 independent degrees of freedom to determine the extent of kinematic information that can be decoded using ECoG signals from human patients. Previous studies using magnetoencephalography and ECoG have found movement-related local motor potential (LMP), mu, beta, and high gamma band activity with the LMP amplitude and high gamma band most related to directional information (Miller, Leuthardt et al. 2007, Waldert, Preissl et al. 2008, Anderson, Blakely et al. 2012, Yeom, Kim et al. 2013). Therefore, we hypothesized that mu and beta band power changes could be utilized to decode movement and rest, and that the high gamma band and LMP would allow for predicting movement kinematics. Through the use of a hierarchical partial-least squares regression model, this study provides the first demonstration that ECoG signals recorded from human patients can be used to perform offline decoding of 3D kinematics in the form of the speed, velocity, and position of arm reaching movements.
The participants in this study were patients with intractable epilepsy who underwent temporary placement of subdural ECoG arrays for localization of their epileptic foci and mapping of eloquent cortices for pre-surgical planning. ECoG electrodes were implanted for a period of approximately one week (5 days – 14 days), during which time the recordings utilized in this study were collected. The Washington University School of Medicine Institutional Review Board approved the study protocol, and all patients provided written informed consent prior to participating. Table 1 contains demographic and clinical information for all 5 patients. Four patients were right-handed, and all electrode arrays were located in the hemisphere contralateral to the dominant hand. Notably, Patient 2 had weakness of the arm contralateral to the site of electrode implantation due to mass effects from the array. This weakness did not noticeably alter movement kinematics, but limited the number of trials that the patient was able to perform.
Clinical ECoG arrays (PMT Corporation, Chanhassen, MN or Ad-Tech, Racine, WI) were utilized for this study. Electrodes were platinum-iridium disks surrounded in silastic sheets with a diameter of 4 mm (2.3 mm diameter exposed) and an inter-electrode distance of 1 cm. Electrode arrays were configured in 8x8 grids, 1x4 strips, 1x6 strips, or 1x8 strips as shown in Figure 1A. In addition to the cortically facing recording electrodes, a 1x4 or 1x6 strip of electrodes was implanted facing the skull for use as ground and reference signals. Recordings were made using g.USBamp biosignal amplifiers (g.tec, Graz, Austria), which utilized 24-bit resolution analog-to-digital converters, an internal sampling rate of 38.4 kHz, and an internal anti-aliasing filter at 5 kHz. The BCI2000 software package (Schalk, McFarland et al. 2004) was used to record ECoG signals with a sampling rate of 1200 Hz and no additional external filtering.
Electrode locations were solely dependent upon the clinical needs of each patient. Electrode coordinates in atlas space were estimated from lateral radiographs collected after electrode implantation. The getLOC package (Miller, Makeig et al. 2007) was used to approximate electrode coordinates with an accuracy of approximately 1 cm. Figure 1B displays the electrode locations for each of the 5 patients.
To examine the relationship between ECoG signals and reaching movements, a 3D center-out reaching task was used. Hand positions were collected using a Flock-of-Birds six degree-of-freedom motion capture system (Ascension Technology, Shelburne, VT). A single sensor was fixed to the index and middle fingers of the moving arm to track hand position. 3D hand positions were sampled at 37.5 Hz. Kinematic data was recorded and synchronized with ECoG signals using a custom-programmed Flock-of-Birds filter that was integrated into the BCI2000 system.
The center-out reaching task consisted of cued reaches to 8 targets positioned at the corners of a physical cube with 50 cm long sides that was set in front of the patient. All reaches began from a target at the center of the cube and progressed to one of the 8 corners of the cube. LED lights were placed at the center target and each of the external targets to provide stimulus cues and reward feedback synchronized to the ECoG and kinematic recordings through a custom-built microcontroller circuit that was interfaced to BCI2000 via a USB interface and custom-programmed application module. During task performance, patients were seated in their hospital bed in a semi-recumbent position with the center target placed at the patient’s midline approximately 40 cm away from their chest. Figure 1C shows the physical apparatus used for the reaching task.
At the beginning of each session, the task was calibrated to determine the location of the target positions and to account for any limitations in patient-specific range-of-motion. Figure 1D shows the time-course for a correct trial. Each trial began with a visual cue for patients to move their hand to the center hold position at which time a hold-A period began, lasting for 500 ms for Patient 1 and 1000 ms for Patients 2–5. During the hold-A period no other stimuli related to the current target was provided. Next, a 2 second plan period began, during which time one of the external targets was illuminated and patients were instructed to plan but not initiate a reaching movement to the target indicated. At the conclusion of the plan period, the indicated external target changed colors, cueing the patient to initiate a reaching movement to the target. Upon reaching the external target, the LED at the specific target turned to green, indicating that the patient had correctly reached the target and cueing the beginning of a 500 ms hold-B period. At the conclusion of the hold-B period, the center target and each of the 8 external targets were illuminated in green to indicate a successful trial completion. If patients reached to an incorrect target or did not reach the correct target within the 4 second time period allowed for the movement, the trial was aborted and all LED lights were illuminated in red to indicate an unsuccessful trial. In Patients 4 and 5, a trial was also aborted if the patient moved before the end of the hold-A, plan, or hold-B periods. The 8 targets were presented in a random order and patients completed multiple runs with 2–4 trials to each target for a total of 16–32 total trials per run. The total number of trials collected and duration of each run were adjusted based upon each patient’s stamina and comfort and are shown in Table 1.
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 8x8 channel ECoG grids spanning multiple amplifiers, individual common averages were calculated for each recording amplifier. For 1x4, 1x6, and 1x8 channel strips, individual common averages were calculated for each strip of electrodes. Signals were band-pass filtered between 0.1 Hz 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 Hz 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 7 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, Stancak et al. 1997). These 7 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 local motor potential, which is obtained by low-pass filtering time domain signals, has been shown to contain movement-related information (Schalk, Kubanek et al. 2007, Pistohl, Ball et al. 2008, Hotson, Fifer 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 7 kinematic parameters was examined for any consistent relationships.
To examine the ability to decode reaching movement kinematics using ECoG signals, the datasets were each divided into training and testing sets. We generated training sets by a simple random sampling of 7/8th of the trials, the remaining trials were held out as a test set. Through the random sampling, each trial was included in multiple test sets. The full training and testing sets were constructed by concatenating each trial from 2 seconds before the onset of movement until the end of the trial.
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, Allison et al. 2010, Flamary and Rakotomamonjy 2012, Aggarwal, Mollazadeh 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.
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 8 features (7 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 seconds and 1 second 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 x 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 λ1 and λ2 are the hyperparameter weights associated with the l1 and l2 norms respectively.
After finding the set of weights that minimize the loss function above, the output of the model was calculated as shown in Equation 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, 7 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 partial-least squares (PLS) regression model to decode continuous traces of kinematic data from ECoG in non-human primates (Chao, Nagasaka 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, Ruhe et al. 1984). For our model, the inputs consisted of the z-scored activity for each of the 7 canonical frequency bands and the LMP at each channel. Individual feature vectors were produced for time lags between −500 ms 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.
Within the training set, we trained models with 1–24 latent features and used 7-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: Vx, Vy, Vz, 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, Nagasaka et al. 2010, Hotson, Fifer 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 (5 patients x 7 kinematic variables x 2 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, Meinecke et al. 2014). Therefore, the decoding model weights were converted to activation patterns describing the encoding of kinematics by ECoG signals (Equation 4) (Haufe, Meinecke 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.
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.
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.
The importance of the 8 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 8-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.
All patients were able to consistently and accurately perform reaching movements to the target locations with no differences in movement times to between targets. To evaluate the true dimensionality of the reaching movements performed, the percent of variance explained by the principle components of the seven kinematic parameters used: speed, velocity (Vx, Vy, Vz), and position (X, Y, Z) is shown in Table 2. The first principle component explained at most 35% of the variance and each of the first four principle components explained at least 11% of the variance. Additionally, the average cross-correlation between kinematic parameters (Table 3) shows that, as expected for a center-out reaching task, while each directional component of velocity is positively correlated with its respective position, speed is not correlated with any direction and there are no cross-directional correlations.
The first level of the decoding model was a logistic regression to classify each time point as movement or rest. While the ability to classify movement and rest periods using ECoG is not novel in itself, it is important to demonstrate as the first component of the hierarchical decoding model used. Exemplar binary predictions of movement and rest are shown for a test set from a single patient in Figure 3. The predicted movement states match the actual labels closely, showing a good ability to predict whether a patient is moving or not using ECoG. Accuracies of predictions of movement and rest are shown in Table 4. Final prediction accuracies were compared to chance prediction levels produced with two surrogate datasets derived by shuffling the temporal relationship between movement class labels and ECoG data, or by shuffling the channel and feature assignments of the model weights. Supplementary Figure S1 shows exemplar surrogate predictions for a single patient. The true prediction accuracy is significantly higher than accuracy produced by either of the surrogate datasets in all patients after Bonferroni correction for the total number of comparisons. Because the movement and rest classes were unbalanced, prediction accuracies produced by the surrogate methods were similar to the overall proportion of movement and rest classes in the dataset. This can be seen in Figure S1 as the surrogate models classify the majority of windows as rest.
The second level of our machine learning model was composed of 2 PLS regressions, one to describe the relationship between ECoG features and speed, velocity, and position during non-movement periods, and a second model to characterize the relationship between ECoG signals and kinematics during movement periods. To examine the ability of this hierarchical PLS regression model to predict the time courses of kinematics during novel time periods, we trained each of the PLS regression models on a training set and used a separate test set to evaluate the prediction accuracy. We first evaluated the ideal performance of the PLS regression in isolation from the movement classification and later evaluated the combined model accuracy.
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 (3 directions, 5 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.
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 4 of the 5 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.
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 8 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.
We sought to evaluate the importance of individual ECoG feature types and cortical locations. Average model weights for the logistic regression model used to classify movement and rest periods across patients are shown for selected features in Figure 7. Beta band and LMP features, in particular, demonstrate the largest amplitude prediction weights with negative weights centered in primary sensorimotor areas, indicating that decreases in sensorimotor cortex beta band power and LMP amplitude are associated with movement periods. The lags for the largest amplitude weights were between 0 ms and 150 ms with neural activity leading the prediction.
To evaluate the importance of individual features and cortical regions for the prediction of continuous kinematic parameters, we examined the movement period PLS regression models. In addition to the decoding model weights, activation patterns were calculated to examine kinematic encoding at each location and are plotted in Figure 8. While the largest weights were located over sensorimotor regions, large amplitude weights are also observed from diffuse regions of cortex. In contrast, the activation patterns are much better localized to sensorimotor regions with a clear peak in the superior motor cortex. While each dimension (X, Y, and Z) was predicted well, at the current electrode size, no cortical representation of directionality was observed. Across temporal and spectral features, the maximum normalized weight amplitude and percent of the largest 1% of prediction weights demonstrated similar trends (Figure 9). In particular, LMP and high gamma features in the 65–95 Hz and 130–175 Hz ranges represented the most important prediction weight features. For the activation patterns, strong encoding of kinematics was observed for the LMP, mu band, beta 1 band, and high gamma bands. Additionally, while both causal and non-causal temporal lags were considered, the peak activation patterns were observed at lags between 250 ms and 400 ms with neural activity leading the predicted kinematics.
To simulate on-line BCI performance, a model was trained for patient 5 without time-locking to individual trials. Patient 5 was chosen because this patient had the most accurate time-locked predictions (Table 6). Figure 10 shows an exemplar period of this pseudo real-time, continuous prediction spanning several consecutive trials. A video demonstrating trajectory predictions made by concatenating the predicted velocity vectors is included in the on-line supplemental materials.
This study demonstrates that kinematic parameters, in particular speed, velocity, and position, of 3D reaching movements can be decoded from ECoG signals in human patients. The ability to decode kinematics of reaching movements with several degrees of freedom builds upon previous demonstrations that ECoG signals can be used to decode movement trajectories (Schalk, Kubanek et al. 2007, Pistohl, Ball et al. 2008, Sanchez, Gunduz et al. 2008, Ganguly, Secundo et al. 2009, Chao, Nagasaka et al. 2010, Shimoda, Nagasaka et al. 2012, Marathe and Taylor 2013) and to control closed-loop BCIs (Leuthardt, Schalk et al. 2004, Wilson, Felton et al. 2006, Felton, Wilson et al. 2007, Schalk, Miller et al. 2008, Rouse and Moran 2009, Rouse, Williams et al. 2013). Although a few previous studies have decoded movement trajectories from movements that were not constrained to two dimensions (Hotson, Fifer et al. 2012, Chen, Shin et al. 2013, Hotson, Fifer et al. 2014), the movements used were correlated with speed in at least one dimension, reducing the dimensionality of the movements. In the task used in this study, the first kinematic principle component explained at most 35% of the variance in any patient and the first 4 principle components each explained at least 11% of the variance in all patients, indicating that the movements performed and decoded included multiple independent dimensions. Therefore, we believe that this study provides the first true demonstration of the ability to use ECoG recordings from human patients to decode 3D kinematics. While 3D kinematic decoding has also been demonstrated from ECoG recordings in non-human primates (Chao, Nagasaka et al. 2010, Shimoda, Nagasaka et al. 2012), the demonstration of kinematic decoding in humans substantially adds to these findings. Non-human primates typically require extensive periods of behavioral training prior to the recording of neural activity, leading to more consistent motor movements than would be expected in real-world behavior. Additionally, as long-term practice of a motor task leads to an increase in the size of the cortical area activated during task performance (Karni, Meyer et al. 1995), decoding of behavioral intentions during a trained motor task may not generalize to more typical novel behaviors.
In addition to spectral power changes, the LMP signal was striking for the strength of its PLS model weights and activation patterns. As the LMP is derived from low-pass filtering ECoG signals, it is likely that there is a correspondence to movement-related cortical potentials such as the bereitschaft potential that are characterized by a negative shift in signal amplitude with an increasing negativity immediately before the movement onset (Kornhuber and Deecke 1965, Deecke, Scheid et al. 1969, Shibasaki and Hallett 2006). Similarly, the LMP signal was negatively correlated with movement speed in sensorimotor regions. A number of previous studies have shown that LMP amplitude can be used to decode position and velocity during arm movements (Schalk, Kubanek et al. 2007, Pistohl, Ball et al. 2008, Ganguly, Secundo et al. 2009, Hotson, Fifer et al. 2014), therefore the large amplitude of prediction weights and activation patterns is unsurprising.
In designing our hierarchical PLS regression model, we chose to use the logistic regression and PLS regression methods in part because they are both linear. Therefore, as the ECoG features were transformed to z-score values prior to training the prediction model, we were able to easily interpret and compare the importance of each cortical location and feature. As the prediction accuracy was improved when real movement and rest labels were used, future work could improve the model accuracy either by improving the classification of movement and rest, or through improving the kinematic decoding. Potential methods to improve decoding accuracy include: calculating spatial filters through unsupervised methods such as independent component analysis (Makeig, Bell et al. 1996, Kachenoura, Albera et al. 2008), supervised learning methods such as common spatial patterns (Koles, Lazar et al. 1990, Koles 1991, Ramoser, Muller-Gerking et al. 2000, Blankertz, Tomioka et al. 2008, Ince, Gupta et al. 2010, Marathe and Taylor 2013), using a model such as a kalman filter that incorporates temporal filtering into the prediction model (Wu, Black et al. 2003, Marathe and Taylor 2013), or by using an alternative linear or non-linear machine learning method (Wessberg, Stambaugh et al. 2000). It should be noted that these alternative methods are not guaranteed to increase prediction accuracies as they may be poorly suited to fit the relationship between ECoG activity and movement kinematics. Additionally, although several variants of the common spatial patterns algorithm have been designed to reduce overfitting (Blankertz, Tomioka et al. 2008, Foldes and Taylor 2011, Lotte and Guan 2011, Falzon, Camilleri et al. 2012, Samek, Vidaurre et al. 2012), these algorithms may still be susceptible to overfitting due to the large feature space and limited amount of training data available. Along with speed, velocity, and position, there are a number of other movement parameters that describe movements such as joint angles or muscle activations. One factor in choosing to begin by decoding speed, velocity, and position is their prevalence in prior studies of single-unit motor control (Georgopoulos, Kalaska et al. 1982, Moran and Schwartz 1999, Wang, Chan et al. 2007) and movement decoding with ECoG (Schalk, Kubanek et al. 2007, Pistohl, Ball et al. 2008, Ganguly, Secundo et al. 2009, Chao, Nagasaka et al. 2010, Marathe and Taylor 2013). While we have shown that speed, position, and velocity can be decoded from ECoG recordings in human patients, algorithms designed to decode alternative components of reaching movements may lead to increases in the accuracies of decoded trajectories.
Along with alternative decoding algorithms, the ability to decode movement kinematics would likely be improved by optimizing the ECoG implants. As the ECoG arrays were implanted for clinical purposes, the location, size, and extent of electrode coverage were all based upon clinical necessity. Recently, it has been found that ECoG electrodes with smaller sizes and spacing can be used to examine neural activity with increased spatial resolution (Leuthardt, Freudenberg et al. 2009, Wang, Degenhart et al. 2009, Kellis, Miller et al. 2010). Additionally, as the arrays used here missed the most superior portions of the motor cortex, optimizing implant locations for decoding neural activity through the use of pre-operative functional imaging techniques (Wang, Collinger et al. 2013) would likely improve decoding accuracy. Taken together, we believe the level of accuracy in this study represents a lower bound that can be improved upon by future studies optimizing electrode arrays and implant locations for BCI systems in patient populations. While these clinical electrode implants required large craniotomies and were highly invasive, future generations of BCI-specific ECoG arrays are likely to substantially reduce the invasiveness of implant procedures (Schalk and Leuthardt 2011).
This study shows that arm movement kinematics have broad cortical representations that can be recorded with ECoG electrode arrays. There are several potential alternative considerations that should be noted. First, due to the margin of error in localizing electrodes onto cortex and intersubject variability in neural anatomy, there is some variance inherent in the localization of decoder weights onto an atlas brain. Next, all of the recordings were made in chronic epilepsy patients. Although care was taken to ensure that all seizures occurred at least 2 hours before or after recordings were made, and that trials with interictal activity were removed prior to analysis, it is difficult to determine if the results of this study were affected by the patient population used. As 4 of the 5 patients had epileptic foci located in the temporal lobe away from sensorimotor areas, we believe that the results were not significantly affected by focal epileptic activity. Additionally, as the reaches used involved the entire arm, it is uncertain if these results will generalize to isolated movements of more distal body parts. Because of this, the results of this study must be interpreted within the context of the activity of the entire motor system and not of individual musculature. Regardless of the role of stabilizing movements, the fact that we can decode information with multiple degrees of freedom still demonstrates the potential for BCI system development. Our decoding model used neural activity both leading and lagging kinematic movements, which would not be possible in a causal BCI system, however, the ability to use ECoG signals to decode movement kinematics underscores the degree of information that can be decoded from ECoG activity, representing a significant step towards the development of ECoG BCI applications. Finally, the broad distribution of decoder weights may function either to amplify movement-related neural activity or to remove non-movement related covariant activity. The activation patterns, representing the features with the strongest relationship with kinematics, were more localized to sensorimotor regions with a clear peak in the more superior areas of motor cortex as expected.
While we have shown that ECoG signals from human patients can be used to decode kinematics of 3D reaches, the optimal method to evaluate the accuracy of a BCI system is ultimately through closed-loop experiments. Open-loop decoding of motor activity allows us to assess the feasibility of BCI algorithms, but there are a number of important differences between open-loop and closed-loop methods. In particular, closed-loop BCI control provides feedback to the user, allowing them to learn to adapt their neural activity to improve the accuracy of the BCI system (Taylor, Tillery et al. 2002, Jarosiewicz, Chase et al. 2008, Rouse and Moran 2009, Rouse, Williams et al. 2013). Additionally, co-adaptively modifying the decoding model in conjunction with the subject’s ability to learn has led to further improvements in on-line BCI control (Williams 2013), indicating that the open-loop accuracies demonstrated in this work could be further improved upon during closed-loop BCI control. Because decoding accuracy can improve though closed-loop adaptation, BCI systems should be designed to enhance the ability of the BCI user to learn and adapt. For example, although open-loop decoding accuracies can be maximized by using algorithms with slow temporal update times, closed-loop decoding accuracies are maximized by algorithms with fast temporal updates allowing for better on-line adaptation (Cunningham, Nuyujukian et al. 2011). Therefore, it will be important for future studies to “close the loop” in order to ensure that algorithms that increase the accuracy of off-line kinematic decoding will also have corresponding increases in BCI control performance. A further question is how valuable the local motor potential (LMP) will be for closed-loop control. In line with previous studies (Schalk, Kubanek et al. 2007, Pistohl, Ball et al. 2008, Hotson, Fifer et al. 2014), we found that the LMP amplitude had a significant contribution to the decoding of kinematics. Although the LMP signal is clearly important for movement decoding, very few studies have used time domain signals for BCI control (Kennedy, Andreasen et al. 2004, Kennedy, Kirby et al. 2004). Because the LMP is derived from low-pass filtering ECoG signals, it changes on a relatively slow time scale. This slow time scale may make it difficult for patients to adapt their LMP signals during BCI control. Future studies using LMP signals within closed-loop BCI systems will be important to understand the advantages and disadvantages of the LMP signal for BCI control.
Taken together, the results of this study provide an exciting demonstration of the ability to use ECoG recordings to decode 3D kinematics of reaching movements in human patients. This underscores the potential to utilize ECoG to develop BCI systems with multiple degrees of freedom for translation into multiple patient populations.
This publication was supported by the James S. McDonnell Foundation and the Washington University Institute of Clinical and Translational Sciences grants UL1 TR000448 and TL1 TR000449 from the National Center for Advancing Translational Sciences.
ECL has stock ownership in the company Neurolutions. DTB has received consulting fees from the company Neurolutions.