Hidden Markov model and support vector machine based decoding of finger movements using electrocorticography

ECoG data were recorded from four patients who volunteered to participate (aged 18–35 years, right handed male). All patients received subdural electrode implants for pre-surgical planning of epilepsy treatment at UC San Francisco, CA, USA. The electrode grid was solely placed based on clinical criteria and covered cortical pre-motor, motor, somatosensory, and temporal areas. Details on the exact grid placement for each subject can be obtained from supplementary figures SI1–SI4 available from stacks.iop.org/JNE/10/056020/mmedia. The study was conducted in concordance with the local institutional review board approved protocol as well as the Declaration of Helsinki. All patients gave their informed consent before the recordings started. The recordings did not interfere with the treatment and entailed minimal additional risk for the participant. At the time of recording all patients were off anti-epileptic medications. They were withdrawn when they were admitted for ECoG monitoring.

The recordings were obtained from a 64-channel grid of 8 × 8 platinum–iridium electrodes with 1 cm centre-to-centre spacing. The diameter of the electrodes was 4 mm, of which 2.3 mm were exposed (except for subject 4 who had a 16 × 16 electrode grid, 4 mm centre-to-centre spacing). The ECoG was recorded with a hardware sampling frequency of 3051.7 Hz and was then down-sampled to 1017 Hz for storage and further processing. During the recordings patients performed a serial reaction time task in which a number appeared on a screen indicating with which finger to press a key on the keyboard (1 indicated thumb, 2: index, 3: middle, 5: little finger; the ring finger was not used). The next trial was automatically initiated with a short randomized delay of (335±36) ms after the previous key press (mean value for S1, session 1; others in similar range). The patients were instructed to respond as rapidly and accurately as possible, with all fingers remaining on a fixed key during the whole 6–8 min run. For each trial the requested as well as actually used finger was recorded. The rate of correct button presses was at or close to 100% for all subjects. The classification trials were labeled according to the finger that was actually moved. Each patient participated in 2–4 of these sessions (table 1). Subjects performed slight finger movements requiring minimal force. Topographic finger representations extend approximately 55 mm along the motor cortex (Penfield and Boldrey 1937). The timing of stimulus presentations and button presses was recorded in auxiliary analogue channels synchronized with the brain data.

Due to clinical constraints, the time available for ECoG data acquisition is restricted, providing a challenging dataset for any classifier. The numbers of trials for classification are listed in table 1. For pre-processing the time series was first highpass filtered using a cut-off frequency at 0.5 Hz as well as notch filtered around the power line frequency (60 Hz) using frequency domain filtering. The electric potentials on the grid were re-referenced to common average reference (CAR).

Next, the time series were visually inspected for the remaining intervals and artefacts caused by the measurement hardware, line noise, loose contacts and sections with signal variations not explained by normal physiology were excluded from the analysis. Channels with epileptic activity were removed since epileptic spiking has massive high-frequency spectral activity and can distort informative features in this range. We also removed any epochs with spread of the epileptic activity to normal brain sites to avoid high-frequency artefacts. The procedure was chosen to ensure comparability across subjects and sessions. Trial rejection was carried out without any knowledge of the corresponding class labels and in advance of any investigations on classification accuracies. Finally, the trials were aligned with respect to the detected button press being the response for each movement request.

The methods described in the following assemble an overall test framework as illustrated in figure 1. The classifier, being the main focus here, is embedded as the core part among other supporting modules including pre-processing, data grouping according to the testing scheme, channel selection and feature extraction. The framework includes two separate paths for training and testing data to ensure that no prior knowledge about the test data is used in the training process.

Zoom In Zoom Out Reset image size

Two different feature types have been extracted for this study. First, time courses of low-frequency components were obtained by spectral filtering in the Fourier domain. These time courses were limited to a specific interval of interest around the button press containing the activation most predictive for the four classes. The chosen frequency range of the lowpass filter covers a band of low-frequency components that contain phase-locked event-related potentials (ERPs) and prominent local motor potentials (LMP) (Kubanek et al 2009). These are usually extracted by averaging out uncorrelated noise over multiple trials (Makeig 1993). Here we focus on single-trial analysis which may not reflect LMP activity (Picton et al 1995). As a consequence, we refer to this feature type as low-frequency time domain (LFTD) features. Boundaries of time interval (length and location relative to the button press) and lowpass filter cut-off were selected by grid search for each session and are stated in the results. The data sampling was reduced to the Nyquist frequency.

Second, event-related spectral perturbation (ERSP) was measured using a sliding Hann-window approach (window length n_w = 260 ms). For each window the square root of the power spectrum was computed by fast Fourier transform (FFT). The resulting coefficients were then averaged in a frequency band of interest. Taking the square root to obtain the amplitude spectrum assigned a higher weight to high-frequency components compared to the original power spectrum. The resulting feature sequence can be used as a measure for time-locked power fluctuation in the high gamma band, neglecting explicit phase information (Graimann et al 2004, Crone et al 1998, Makeig 1993).

The time interval relative to the button press, and the window overlap and exact low and high cut-off frequency of the high gamma band most discriminative for the four classes were selected by grid search. The resulting parameters for each session are stated in the results section.

A combined feature space including both high gamma (HG) as well as LFTD features was also evaluated. The ratio p of how many LFTD feature sequences are used relative to the number of HG sequences was obtained by grid search. Subset selection beyond this ratio is described in the next subsection.

2.4.1. General remarks on subset selection

2.4.2. DB index

The index corresponds to a cluster separation index measuring the overlap of two clusters in an arbitrarily high dimensional space, where ${\boldsymbol r} \in \mathbb{R}^T$ is one of N sample vectors belonging to cluster C_i:

$$\begin{equation} {\boldsymbol \mu }_{\bf i} = \frac{1}{N}\sum\limits_{{\bf r} \in C_i } {\bf r} \end{equation} \tag{ 1 }$$

$$\begin{equation} d_i = \frac{1}{N}\sum\limits_{{\bf r} \in C_i } {\| {{\bf r} - {\boldsymbol \mu }_{\bf i} } \|} \end{equation} \tag{ 2 }$$

$$\begin{equation} R_{ij} = \frac{{{\rm d}_i + {\rm d}_i }}{{\| {{\boldsymbol \mu }_{\bf i} - {\boldsymbol \mu }_{\bf j} } \|}}\quad i,j \in N^ + ,\,i \ne j \end{equation} \tag{ 3 }$$

Using the cluster centroids ${\boldsymbol \mu}_i \in \mathbb{R}^T$ and spread d_i, matrix ${\boldsymbol R} = \{ {R_{ij} } \}_{i = 1 \ldots N_c ,j = 1 \ldots N_c}$ contains the rate between the within-class spread and the between-class spread for all class combinations. The smaller the index, the less two clusters overlap in that feature space. The index is computed for all features available.

2.4.3. Feature selection procedure

2.5.1. Hidden Markov models

Arising from a state machine structure, hidden Markov models (HMMs) comprise two stochastic processes. First, as illustrated in figure 2, a set of Q underlying states is used to describe a sequence O_t of T O-dimensional feature vectors. These states q_t ∈ {s_i}_{i = 1...Q} are assumed to be hidden, while each is associated with an O-dimensional multivariate Gaussian M-mixture distribution to generate the observed features—the second stochastic process. Both, the state incident q_t as well as the observed output variable o_t at time t are subject to the Markov property-–a conditional dependence only on the last step in time t − 1.

Zoom In Zoom Out Reset image size

Among other parameters the state sequence structure and its probable variations are described by prior probabilities P(q₁ = s_i) and a transition matrix

$$\begin{equation*} {\bf A} = \{ {a_{ij} } \}_{i,j = 1 \ldots Q} = \{ {P(q_{t + 1} = s_j | {q_t = s_i )} } \}_{i,j = 1 \ldots Q} . \end{equation*}$$

A more detailed account of the methodological background is given by Rabiner (1989).

For each of the four classes one HMM had to be defined, which was trained using the iterative Baum–Welch algorithm with eight expectation maximization (EM) steps (Dempster et al 1977). The classification on the test set was carried out following a maximum likelihood approach as indicated by the curly bracket in figure 2.

In order to increase decoding performance, the above standard HMM structure was constrained and adapted to the classification problem. This aims at optimizing the model hypothesis to give a better fit and reduces the number of free parameters to deal with limited training data.

2.5.2. Model constraints

2.5.3. Model structure

2.5.4. Initialization

2.5.5. Implementation

2.5.6. Support vector machines

As a gold-standard reference we used a one-vs-one SVM for multi-class classification. In our study, training samples come in paired assignments (x_i, y_i)_{i = 1...N}, where x_i is a 1 × (T · O) feature vector derived from the ECoG data for one trial with its class assignment y_i = {–1,1}. Equation (4) describes the classification rule of the trained classifier for new samples x. A detailed derivation of this decision rule and the rationale behind SVM-classification is given e.g. by Hastie et al (2009):

$$\begin{equation} F({\bf x}) = {\mathop{\rm sgn}} \left( {\sum\nolimits_i {\alpha _i y_i K({\bf x}_{\bf i} {\bf ,x}) + \beta _0 } } \right) \end{equation} \tag{ 4 }$$

The α_i ∈ [0, C] are weights for the training vectors and are estimated during the training process. Training vectors with weights α_i ≠ 0 contribute to the solution of the classification problem and are called support vectors. The offset of the separating hyperplane from the origin is proportional to β₀. The constant C generally controls the trade-off between classification error on the training set and smoothness of the decision boundary. The numbers of training samples per session (table 1) were in a similar range compared to the resulting dimension of the SVM feature space T · O( ∼ 130-400). Influenced by this fact, the decoding accuracy of the SVM revealed almost no sensitivity with respect to regularization. Since grid search yielded a stable plateau of high performance for C > 1, the regularization constant was fixed to C = 1000. The linear kernel, which has been used throughout the study, is denoted by K(x_i, x). We further used the LIBSVM (www.csie.ntu.edu.tw/∼cjlin/libsvm/) package for Matlab developed by Chih-Chung Chang and Chih-Jen Lin (Chang and Lin 2011).

In order to evaluate the findings presented in the following sections, a five-fold cross-validation (CV) scheme was used (Lemm et al 2011). The assignment of a trial to a specific fold set was performed using uniform random permutations.

Each CV procedure was repeated 30 times to average out fluctuations in the resulting recognition rates arising from random partitioning of the CV sets. The number of repetitions was chosen as a compromise to obtain stable estimates of the generalization error, but also to avoid massive computational load. The training data was balanced between all classes to treat no class preferentially and to avoid the need for class prior probabilities. Feature subset selection and estimation of the parameters of each classifier was performed only on data sets allocated for training in order to guarantee a fair estimate of the generalization error. Due to prohibitive computational effort, the feature spaces have been parameterized only once for each subject and session and are hence based on all samples of each subject-specific data set. In this context preliminary experiments indicated minimal differences to the decoding performance (see table 2) when fixing the feature definitions based on the training sets only. These experiments used nested CV with an additional validation set for the feature space parameters.

In order to identify optimal parameters for feature extraction a grid search for each subject and session was performed. For this optimization a preliminary feature selection was applied aiming at a fixed subset size of O = 8, which was found to be within the range of optimal subset sizes (see figure 6 and table 3). Results for all sessions are illustrated by the bar plots in figure 3. A sensitivity of the decoding rate has been found for temporal parameters such as location and width for the time windows (ROIs) containing the most discriminative information about LFTD or HG features (figures 3(a) and (c)).

Zoom In Zoom Out Reset image size

The highest impact on accuracy was found for the HG ROI (mean optimal ranges: S1: [−187.5 337.5] ms; S3: [−400 475] ms; S4: [−280 420] ms).

Further sensitivity is given for spectral parameters namely the high cut-off frequency for the LFTD lowpass filter and the low and high cut-off for the HG bandpass (figures 3(b) and (d)). For the former, optimal frequencies ranged from 11–30 Hz and the optimal frequency bands for the latter all lay within 60–300 Hz, depending on the subject. There was hardly any impact on the decoding rate when changing the sliding window size and window overlap for the HG features. These were hence set to fixed values. According to these parameter settings, the mean numbers of samples T were 26.5 (subject 1), 21.3 (subject 2), 31 (subject 3) and 22 (subject 4).

An example of the grid search procedure is shown in supplementary figures SI5–SI7 (available at stacks.iop.org/JNE/10/056020/mmedia). The plots show varying decoding rates for the S1 session 1 within individual parameter subspaces. The optimal parameters illustrated in figure 3 represent a single point in these spaces at peaking decoding rate. Based on these parameters figure 4 illustrates a typical result for a LFTD and HG feature sequence extracted from the shown raw ECoG signal (S1, session 2, channel 23, middle finger). All feature sequences for this trial—consisting of the entire set of channels-–are given in supplementary figure SI8.

Zoom In Zoom Out Reset image size

Using the resulting parameterization, subsequent CV runs sought to identify the optimal subset by tuning O based on the training set of the current fold. Figure 5 illustrates the rate at which a particular channel was selected for each subject. It evaluates the optimal subsets O across all sessions, folds and several CV repetitions. Results are exemplarily shown for high gamma features. Full detail on subject-specific grid placement and channel maps for LFTD and HG features can be found in supplementary figures SI1–SI4 (available at stacks.iop.org/JNE/10/056020/mmedia).

Zoom In Zoom Out Reset image size

The channel maps reveal localized regions of highly informative channels (dark red) comprising only a small subset within the full grid. While these localizations are distinct and focused for high gamma features, informative channels spread wider across the electrode grid for LFTD features. This finding is in line with an earlier study done by Kubanek et al (2009) and implies that selected subsets for LFTD features tend to be less stable across CV runs. Combined subsets were selected from both feature spaces. The selections indicated a preference for a higher proportion of LFTD features (table 3). This optimal proportion was determined by grid search.

The optimal absolute subset size O was selected by successively adding channels to the subset as defined by our selection procedure and identifying the peak performance. For the high gamma case, figure 6 already shows a sharp increase in decoding performance for a subset size below five. Performance then increases less rapidly until it reaches a maximum usually in the range of 6–8 channels for the HG features (see table 3 for detailed information on the results for the other features). This behaviour again reflects the fact that a lot of information is found in a small number of HG sequences.

Zoom In Zoom Out Reset image size

Increasing the subset size beyond this optimal size leads to a decrease in performance for both SVM and HMM. This coincides with an increase of the number of free parameters in the classifier model depending quadratically on the number of channels for HMMs (figure 6, bottom).

Based on the configuration of the feature spaces and feature selection, a combination space containing both LFTD and HG features yielded the best decoding accuracy across all sessions and subjects. Figure 7 summarizes the highest HMM performances for each subject averaged across all sessions and compares them with the corresponding SVM accuracies. Except for subject 2, these accuracies were close to or even exceeded 90% correct classifications. However, figure 8 shows that accuracies obtained only from the HG space come very close to this combined case, implying that LFTD features contribute limited new information. Using time domain features exclusively decreased performance in most cases. An exception is given by subject 2 where LFTD features outperformed HG features. This may be explained by the wider spread of LFTD features across the cortex and the suboptimal grid placement (little coverage of motor cortex, supplementary figure SI2 available at stacks.iop.org/JNE/10/056020/mmedia). The channel maps show informative channels at the very margin of the grid. This suggests that grid placement was not always optimal for capturing the most informative electrodes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Differences in decoding accuracy compared to the SVM were found to be small (table 4) and revealed no clear superiority of any classifier (figure 9). Excluding subject 2, the mean performance difference (HMM-SVM) was –0.53% for LFTD, 0.1% for HG and –0.23% for combined features. Accuracy trends for individual subjects are similar across sessions and vary less than across subjects (figure 9). Note that, irrespective of the suboptimal grid placement, subject 2 took part in four experimental sessions. Here, subject-specific variables, e.g. grid placement, were important determinants of decoding success. Thus we interpret our results at the subject level and conclude that, in three out of four subjects, SVM and HMM decoding of finger movements led to comparable results.

Zoom In Zoom Out Reset image size

A clear superiority for SVM classification was only found for subject 2, for whom the absolute decoding performances were significantly worse than those of the other three subjects.

For the single spaces, LFTD and HG, table 4 shows the mean performance loss for the unleashed HMM. That refers to a case wherein the Bakis model constraints are relaxed to the full model. For all except two sessions of subject 2 (the full model had a higher performance of 0.2% for session 3 and 0.7% for session 4), the Bakis model outperformed the full model. The results indicate a higher benefit for the LFTD features with an up to 8.8% performance gain for subject 3, session 2. Higher accuracies were consistently achieved for feature subsets deviating from the optimal size. An example is illustrated for LFTD channel selection (subject 3, session 2) in supplementary figure SI9 (available at stacks.iop.org/JNE/10/056020/mmedia). The plot compares the Bakis and full model similar to the procedure in figure 6. While table 4 lists only the improvements for the optimal number of channels, supplementary figure SI9 indicates that the performance gain may be even higher when deviating from the optimal settings.

Finally, figure 10 illustrates the accuracy decompositions into true positive rates for each finger. The results reveal similar trends across subjects and feature spaces. While the thumb and little finger provided the fewest training samples, their true positive rates tend to be the highest among all fingers. An exception is given for LFTD features for subject 2. Equivalent plots for the SVM are shown in supplementary figure SI10 (available at stacks.iop.org/JNE/10/056020/mmedia).

Zoom In Zoom Out Reset image size

For the classification of one feature sequence including all T samples, the HMM Matlab implementation took on average (3.522 ± 0.036) ms and the SVM implementation in C (0.0281 ± 0.0034) ms. The mean computation time for a feature sequence was 0.2 ms ± 2.1µs (LFTD), 1.72 ms ± 37.9 µs (HG) and 1.79 ms ± 20.0 µs (combined space). The results were obtained with an Intel Core i5-2500 K @ 3,4 GHz, 16 GB RAM, Win 7 Pro 64-Bit, Matlab 2012b for ten selected channels.