Sitting and standing intention can be decoded from scalp EEG recorded prior to movement execution

Low frequency signals recorded from non-invasive electroencephalography (EEG), in particular movement-related cortical potentials (MRPs), are associated with preparation and execution of movement and thus present a target for use in brain-machine interfaces. We investigated the ability to decode movement intent from delta-band (0.1–4 Hz) EEG recorded immediately before movement execution in healthy volunteers. We used data from epochs starting 1.5 s before movement onset to classify future movements into one of three classes: stand-up, sit-down, or quiet. We assessed classification accuracy in both externally triggered and self-paced paradigms. Movement onset was determined from electromyography (EMG) recordings synchronized with EEG signals. We employed an artifact subspace reconstruction (ASR) algorithm to eliminate high amplitude noise before building our time-embedded EEG features. We applied local Fisher's discriminant analysis to reduce the dimensionality of our spatio-temporal features and subsequently used a Gaussian mixture model classifier for our three class problem. Our results demonstrate significantly better than chance classification accuracy (chance level = 33.3%) for the self-initiated (78.0 ± 2.6%) and triggered (74.7 ± 5.7%) paradigms. Surprisingly, we found no significant difference in classification accuracy between the self-paced and cued paradigms when using the full set of non-peripheral electrodes. However, accuracy was significantly increased for self-paced movements when only electrodes over the primary motor area were used. Overall, this study demonstrates that delta-band EEG recorded immediately before movement carries discriminative information regarding movement type. Our results suggest that EEG-based classifiers could improve lower-limb neuroprostheses and neurorehabilitation techniques by providing earlier detection of movement intent, which could be used in robot-assisted strategies for motor training and recovery of function.


INTRODUCTION
Robot-assisted therapies have shown promising results, compared to traditional therapy, for functional recovery of movement after injury in the upper and lower extremities (Winchester et al., 2005;Hogan and Krebs, 2011). These neurorehabilitation paradigms could be improved by faster and more robust detection of movement intent where it originates in the brain. Incorporation of a brain machine interface (BMI) can reduce the latency between motor planning in the cortex and activation of a device to execute (or assist) the movement, thereby enhancing the opportunity for brain plasticity and motor recovery (Daly and Wolpaw, 2008). The intuitive nature of a BMI based on signals directly related to intended movement could be advantageous for rehabilitation by expediting adaptation of the brain to the BMI algorithm and the robotic device. Electroencephalography (EEG) provides a non-invasive method for imaging brain activity with enough time resolution to exert control over an assistive device. Many strategies for deploying EEG in a BMI by detecting movement intent (imagined and real) have been reported (Pfurtscheller et al., 1996Wolpaw et al., 2002;Millán et al., 2004;Qin et al., 2004;Hung et al., 2005;Morash et al., 2008). These systems typically leverage one of two phenomena to detect movement intent: event related (de)synchronization (ERD/ERS) and movement related slow cortical potentials (MRPs). ERD, a decrease of power in alpha and beta bands, is typically localized to the contralateral sensorimotor areas before movement while ERS, a power increase, has been observed after movement (Pfurtscheller and Lopes da Silva, 1999). Modulation of these sensorimotor rhythms has been employed for classification of imagined Pfurtscheller and Neuper, 2006) and executed (Morash et al., 2008) movements with some success. ERD has also shown capacity to categorize gross lower extremity tasks, including differentiation of right and left leg motor imagery (Boord et al., 2010) and identification of imagined standing (Zhong et al., 2007).
MRPs are slow negative potentials observed in EEG preceding movement. MRPs can be divided into two segments: the first begins as early as 2 s before movement onset and has been observed over the entire pre-supplementary motor area (SMA), and over the SMA and lateral premotor cortex according to somatotopic organization (Ikeda et al., 1992;Hallett, 1993;Shibasaki and Hallett, 2006;Bai et al., 2011). The second, or late, segment typically has a steeper negative slope and is observed in the contralateral primary motor cortex (M1) and lateral premotor cortex according to precise somatotopic arrangement. These potentials are well established in upper and lower extremity movements both real and imagined (Boschert and Deecke, 1986;Shibasaki and Hallett, 2006). Interestingly, MRPs recorded from EEG preceding toe, foot, and ankle movements tend to be larger on the ipsilateral side of the brain, which is the opposite of upper extremity movements that create larger MRPs on the contralateral side (Brunia and Van Den Bosch, 1984;Boschert and Deecke, 1986). This paradoxical lateralization of the MRP during foot movements may be explained by its localization along the midline deep within the precentral gyrus of the motor cortex, thereby directing electrical current from activation of these cell columns to the opposite hemisphere.
The type and sequence of movement affects MRPs recorded from EEG. MRPs appear to be more pronounced during self-initiated movements compared to triggered movements (Jahanshahi et al., 1995;Cui and MacKinnon, 2009); the difference appears to be further enhanced if the timing of the triggered movements is variable (Jankelowitz and Colebatch, 2002). In the case of finger movements, force level (Slobounov et al., 2002), finger sequence (Bortoletto et al., 2011), and task complexity (Shibasaki and Hallett, 2006) all appear to modulate the MRP. MRP amplitude was found to be highly correlated to joint torque and electromyography (EMG) amplitude during isolated elbow flexion (Siemionow et al., 2000). In the lower extremity, the rate of torque development appears to influence the late MRPs preceding isolated ankle movements (do Nascimento et al., 2006). Slow negative shifts in EEG similar to MRPs have been observed during coordinated movements of the lower extremity, including rising onto the toes (Saito et al., 1996) and self-paced forward postural sway (Slobounov et al., 2005). The direction of gait initiation and stepping has been reported to influence both the slope and magnitude of MRPs (do Nascimento et al., 2005). These previously published studies suggest that slow developing, movement related potentials observed prior to movement may contain discriminative information regarding the movement that is being performed. Further, MRPs appear to provide an appropriate measure for timing of afferent feedback to induce long term potentiation of cortical projections. As demonstrated in the tibialis anterior muscle, only peripheral stimulation delivered at the peak of the MRP increased motor evoked potentials from transcranial magnetic stimulation (TMS) targeting the ankle area of the motor cortex .
Because of their small amplitude and low frequency content, the best way to extract MRPs from EEG recording is to average across many trials of the same movement. Single trial classification of movement intention from MRPs is possible, but achieving high accuracy can be difficult. Classification typically involves several steps, including signal pre-processing, feature extraction, dimensionality reduction, and finally feature classification (Bashashati et al., 2007). Numerous approaches to these steps have resulted in application of many machine learning, feature selection, and pattern recognition techniques for classification of movement intent and direction based on EEG signals (Garrett et al., 2003;Peterson et al., 2005;Bai et al., 2007;Lotte et al., 2007). The first example of a BMI-based spelling device utilized slow cortical potentials derived from a motor imagery task to provide individuals with amyotrophic lateral sclerosis control of a cursor on a screen (Birbaumer et al., 1999). Two individuals were able to achieve accuracies greater than 75% after 327 and 288 training sessions. Recent studies have demonstrated success in utilizing MRPs extracted via low frequency or delta band EEG, including classification of finger movement (Liao et al., 2007), joystick direction (Waldert et al., 2008), wrist movement direction (Vuckovic and Sepulveda, 2008), direction of a center out reaching task (Robinson et al., 2013), and movement intention in a self-paced reaching task (Lew et al., 2012). The latter study showed higher detection accuracy using the lower delta band than alpha (7-13 Hz) or beta (13-20 Hz) bands. MRPs have also been successfully deployed for classification of lower extremity movements. At the ankle, MRPs have been used to detect movement intention in healthy subjects with average accuracy of 82.5% for movement execution, and with slightly lower accuracy for motor imagery (64.5%) and attempted movement in stroke patients (55%) (Niazi et al., 2011). Similar accuracies were reported in a study that did not incorporate an individual-specific training phase (Niazi et al., 2013), further supporting the robustness of MRP as a BMI target. In addition to movement intention, MRPs recorded during imagined plantar flexion have also been used to distinguish between two different rates of torque development (do Nascimento and Farina, 2008). Recent studies have demonstrated that MRPs recorded from EEG can be deployed in real-time BMIs. In one, MRPs preceding imagined ankle dorsiflexion were identified online to trigger electrical stimulation of the tibialis anterior . Not only did this study show feasibility of MRPs for use in a BMI, but it also demonstrated the potential benefits BMMI-based neurorehabilitation since motor evoked potentials from TMS were enhanced following the intervention in healthy individuals. Another study showed that delta band EEG could reliably ascertain ankle movement initiation in real time with a mean latency of 315 ms (Xu et al., 2014).
In addition to detecting and classifying movement type, sparse networks of low frequency EEG have also been successful in decoding kinematics and EMG activity during various movements, including decoding of hand grasping patterns (Agashe and Contreras-Vidal, 2013), hand and finger velocity (Bradberry et al., 2010;Liu et al., 2011;Paek et al., 2014), and muscle synergies during reaching (Beuchat et al., 2013). Additionally, peri-movement neural activity representative of movement direction has been observed in electrocorticographic (ECoG) signals over primary motor, premotor, posterior-parietal, and lateral prefrontal cortex (Ball et al., 2009). Action intention can also be decoded from fMRI data recorded from a wide cortical network, spanning from the parieto-occiptial sulcus through the prefrontal cortex, both preceding and during movement execution (Gallivan et al., 2011). Taken together these studies suggest non-invasive EEG recorded from large areas of the scalp immediately prior to movement execution could carry useful information about movement.
EEG has been used to examine cortical activity during gait, including studies demonstrating that intra-stride changes in spectral power are coupled to gait cycle (Gwin et al., 2011) and that level of user-involvement in robotic-assisted walking alters gaitrelated patterns of electrocortical activity (Wagner et al., 2012). Low frequency EEG also appears to carry useful information regarding walking. A recent study showed that features corresponding to frequencies less than 2 Hz were the most heavily weighted during single trial classification of walking and pointing direction (Velu and de Sa, 2013). Delta-band EEG was used to classify walking intention in one individual with paraplegia using a robotic exoskeleton with accuracy greater than 98%  and to decode lower limb kinematics during walking in healthy individuals (Presacco et al., 2011(Presacco et al., , 2012. MRPs have also been used with a matched filtering technique to detect single-trial step initiation ). An important consideration for application of low frequency EEG to the study of whole-body movements such as walking or sit/stand transition is the presence of movement-related artifacts. A recent study showed similar power spectral density patterns from an accelerometer mounted on the head and from EEG electrodes (Castermans et al., 2014). Interestingly, the patterns were similar only at higher walking speeds, while differences between the accelerometer and EEG were observed at slower speeds. The study did not compare spectral patterns from EEG during walking without the rigid plate and linkage assembly used to mount the accelerometer on the head, so the effect of its mass and inertia remains unknown. Also, the study did not employ active EEG electrodes which provide amplification at the electrode to minimize movement artifacts and increase signal-to-noise ratio. Spatial filtering techniques, such as independent component analysis (Delorme et al., 2007), may be used to isolate gait-related artifact, but the effectiveness of these techniques is still under investigation. In one study, gait-related artifact remained in many independent components of EEG, resulting in development of a template subtraction technique to clean EEG collected during walking (Gwin et al., 2010). This type of template regression would not be appropriate for studying cortical contribution to locomotion because all signals coupled to the gait cycle would likely be removed. Another technique utilizes principal component analysis to compare sliding windows of EEG to a baseline recording, thereby removing high amplitude artifacts (Mullen et al., 2013); this approach may be better suited for removing movement artifacts but has not yet been applied to gait. Thus, the feasibility of utilizing EEG to study cortical activations during whole-body movement tasks is an ongoing area of research. Nevertheless, an inherent advantage of MRPs is their presence in EEG recorded before movement, when motion artifacts are minimized.
In this study we examined the use of non-invasive EEG recorded prior to movement execution to discriminate a user's intent to perform two coordinated whole body movementsrising from a seated to standing posture and lowering from a standing to a seated posture-in a three class problem, where the third class constituted no movement or "quiet"; this class included data collected during quiet standing and quiet sitting. Based on the previous body of evidence regarding the discriminative nature of MRPs with regards to movement, we utilized delta band EEG to build our features for classification. We trained and tested our classifier using time periods before executed movements, as opposed to cue-based imagery, so we could precisely align EEG recordings with movement onset detected from EMG recordings. We studied classification accuracy during two different paradigms: a self-initiated series of stand-to-sit and sit-tostand transitions and transitions which were cued by an audio trigger. Because triggered movements are reported to produce less prominent MRPs (Jankelowitz and Colebatch, 2002;Cui and MacKinnon, 2009), this protocol allowed us to examine the effect of MRP signal to noise ratio on classification accuracy. We utilized time-embedding and concatenation of EEG channels from the time before movement execution to create a feature vector of high dimension to classify the intended movement (standup, sit-down, or quiet). Given the autoregressive nature of EEG signals (Muller et al., 2003) and the underlying neurophysiology (e.g., volume conduction), we assume that the recorded EEG originates from a system with fewer degrees of freedom than our feature vector dimensions, resulting in a manifold data structure. Recent advances in machine learning have resulted in algorithms which preserve the local structure of a manifold data set in a reduced dimensional subspace (Sugiyama, 2007;Li et al., 2012) thereby enhancing the discriminative power of the data set. Based on the observation that information pertinent to movement is contained in low frequency EEG, we hypothesized that applying a locality preserving dimensionality reduction technique to our high dimensional feature vector derived from time-embedded and spatially diverse delta band EEG would reveal its underlying discriminative structure. We coupled this supervised data reduction with a Gaussian mixture model classifier to test if we could reliably ascertain the intended movement of the user from offline analysis of EEG recordings. We believe such a classifier could eventually be deployed in a real-time BMI system to control an assistive device or as a component of a neurorehabilitation paradigm to restore motor control.

DATA COLLECTION
Ten healthy adults (6 male, 4 female) with no history of neurological disease participated in the study after giving informed consent. This study protocol was approved by the Institutional Review Board at the University of Houston. Participants completed two trials of 10 alternating sit-to-stand and stand-to-sit transitions; one trial was self-paced and one trial was cued via audio trigger. Each trial began with the participant standing quietly in an upright posture for 15 s. In the triggered trial, an audio cue (beep) was given after which point the participant initiated a transition to a seated posture. The seated posture was held for a period ranging randomly from 3 to 10 s, after which a second audio cue was given to initiate the transition from sit-to-stand. The standing posture was held for another (random) 3-10 s interval, at which point the process was repeated until 20 transitions (10 of each) were completed. The procedure for the self-paced trial was similar. After 15 s of quiet standing, the participant was instructed via verbal cue to begin the selfinitiated stand-to-sit and sit-to-stand transitions. The participant was instructed to wait for a random interval of 3-10 s before self-initiating the next transition. Finally, the participant was notified by verbal cue once he/she had completed 20 self-initiated transitions.
Time-locked EMG and EEG data were collected simultaneously using a previously developed data collection system (Bulea et al., 2013). Surface EMG (Biometrics, Ltd, Ladysmith, VA) was recorded at 1000 Hz bilaterally from the tibialis anterior, gastrocnemius, biceps femoris, and vastus lateralis. Whole scalp, active electrode, 64-channel EEG (Brain Products, GmbH, Morrisville, NC) were collected at 1000 Hz and labeled by the 10-20 international system. The impedance of each EEG electrode was maintained below 25 k for the entire data collection.

Preprocessing
All data analysis and classifier optimization and evaluation were performed off-line using custom software in Matlab (Mathworks, Natick, MA). The data processing and classification methodology is shown in Figure 1. Peripheral EEG channels susceptible to eye blinks and facial/cranial muscle activity were removed from offline analysis (all channels labeled Fp, AF, FT, T, TP, O, PO, and F5-8, P5-8) resulting in 28 channels being retained for classification. EEG signals were then high pass filtered at 0.05 Hz using a zero-phase 8th order Butterworth filter. Next, we removed transient, high-amplitude artifacts from stereotypical (e.g., eye blinks) and non-stereotypical (e.g., movement, muscle bursts) using an automated artifact rejection method termed Artifact Subspace Reconstruction (ASR) (Mullen et al., 2013) which is available as a plug-in for EEGLAB software (Delorme and Makeig, 2004). ASR uses a sliding window technique whereby each window of EEG data is decomposed via principal component analysis so it can be compared statistically with data from a clean baseline EEG recording, collected here as 1 min of EEG recorded during quiet standing. Within each sliding window the ASR algorithm identifies principal subspaces which significantly deviate from the baseline EEG and then reconstructs these subspaces using a mixing matrix computed from the baseline EEG recording. In this study, we used a sliding window of 500 ms and a threshold of 3 standard deviations to identify corrupted subspaces. After ASR, the cleaned EEG was band pass filtered with a zero phase, 3rd order Butterworth filter from 0.1 to 4 Hz to isolate the delta band activity. The EEG data were then standardized by channel by subtracting the mean and dividing by the standard deviation (z-score).
EMG recordings from the lower extremity muscles were used to determine movement onset of each stand-to-sit and sit-tostand transition. First, the Teager-Kaiser energy operator was applied to each EMG channel to enhance the signal-to-noise ratio for onset detection . Next, each EMG channel was detrended, band pass filtered (15-300 Hz), rectified, and low pass filtered at 3 Hz to compute the linear envelope. Then, the FIGURE 1 | Flow chart describing the EMG and EEG data processing for neural decoding of sitting and standing movement. A threshold of 3 standard deviations was applied to the EMG linear envelope to identify quiet periods and periods of movement (sitting and standing). Only pre-movement epochs (1.5 s before movement to movement onset) and quiet epochs (1.5 s after movement completion to 1.5 s before next movement) were retained for analysis. As a control, a separate decoding analysis using movement epochs (movement onset to 1.5 s after onset) was also performed. Artifact subspace reconstruction (ASR) algorithm, available as a plug-in for EEGLAB software (Delorme and Makeig, 2004), was applied to eliminate artifacts from EEG data during pre-processing. Note that the optimization and evaluation data sets are mutually exclusive.
linear envelope of each muscle was thresholded into a binary signal which was equal to 1 when the envelope exceeded its mean baseline value during quiet standing and sitting by more than 3 standard deviations (Hodges and Bui, 1996) and zero when it was within 3 standard deviations of baseline. The baseline period of EMG activity before each movement was identified a posteriori by visual inspection starting with the initial 15 s of rest before the first movement. The baseline period between each successive sit-to-stand and stand-to-sit transition comprised at least 2 s. Movement onset for each transition was determined when any of the 8 thresholded EMG envelopes transitioned from rest (0) to active (1). Likewise, the end of each movement was determined when all 8 channels returned to rest (0). The algorithmically determined periods of activity were visually inspected for accuracy. Using prior knowledge of the experimental protocol (i.e., the order of the stand-to-sit and sit-to-stand transitions), the periods of muscle activity were labeled as stand-to-sit or sitto-stand. Note that for some trials, gastrocnemius muscles were active during the quiet stance phase and/or biceps femoris EMG was contaminated by artifact from the leg during sitting, thereby increasing the standard deviation in these channels and limiting the ability to determine the true state using that muscle. When these periods of activity/artifact were observed visually, these muscles were removed from the trial; in this case the user activity was assessed using the remaining 6 muscles.
Next, the time-locked EEG and EMG data were downsampled to 200 Hz. EEG data were then epoched into pre-movement, postmovement and quiet periods based on the thresholded (binary) EMG signal. Each pre-movement epoch consisted of data from 1.5 s before movement onset up to movement onset. EEG data from 1.5 s after movement completion until 1.5 s before the next movement onset, with a maximum of 5 s, comprised the quiet epochs. These epochs were then concatenated into a single time series containing alternate periods of quiet and pre-movement. For control purposes, we also created a second time series of data containing concatenated quiet epochs and epochs of EEG from movement onset to 1.5 s after movement onset (post-movement epochs).
The concatenated EEG data sets comprised the three-class classification problem for each trial; each time point of the quiet epochs was labeled as class 0 (quiet) while each time point of each pre-movement epoch was labeled according to the type of movement it preceded: class 1 (stand-to-sit) or class 2 (sit-to-stand). Next, a time-embedded feature matrix was constructed for each trial. Each time point in the feature matrix was a vector composed of 10 lags, corresponding to 50 ms in the past, of EEG data. The number of lags and embedded time interval was chosen based on previous studies demonstrating accurate decoding of movement kinematics from low frequency EEG (Bradberry et al., 2010;Presacco et al., 2011). The feature vector for each time point was constructed by concatenating the 11 lags (the current time point plus the 10 prior) for each channel into a single vector of length 11 × N, where N is the number of EEG channels used for classification (for this study, N = 28). To avoid the problem of missing data, the feature matrix was buffered by starting at the 11th EEG sample of each epoch, resulting in a feature matrix of dimension [M t −L] × [11 × N] for each trial of self-initiated and triggered movements where M t is the number of time points in each trial and L is the number of past time lags multiplied by the number of epochs in each trial (for this study, L = 10 * 41 = 410). On average, there were 18,442 ± 2110 time points in each feature matrix, with exactly 2900 time points for class 1 and 2900 time points for class 2 while the remaining time points represented class 0. For all subjects, the original dimensionality of the feature space was 308 (11 × N).

Dimensionality reduction
Since our EEG-based feature vectors were of relatively high dimension and were composed of time lagged and spatially distributed samples, we assumed our original dataset to represent a manifold which may contain multimodal within-class distributions. Furthermore, we sought to classify gross motor intention and therefore had a limited number of classes (in this case there were three: quiet, stand-to-sit, and sit-to-stand). Thus, we performed dimensionality reduction on our feature matrices to eliminate any redundant features, reduce computational complexity, prevent over-fitting during classifier training and increase classification performance. Many techniques have been reported for dimensionality reduction in EEG based classifiers, including principal component analysis (PCA), linear discriminant analysis (LDA), and genetic algorithm (GA) (Bashashati et al., 2007;Lotte et al., 2007). Consideration of the task, neurophysiology and EEG recording system suggests that a supervised dimensionality reduction technique could improve feature selection for classification purposes. EEG data generally have a low signal-to-noise ratio and unsupervised linear dimensionality reduction techniques may be affected by these signal distortions. PCA reduces dimensionality by maximizing data variance in the projected subspace via a linear transformation. The transformation, dictated by the eigenvectors that correspond to the largest eigenvalues of the data covariance matrix, is unsupervised and can discard useful information for classification that is contained in the lower energy dimensions of the original data (Prasad and Bruce, 2008). In contrast, LDA is a supervised dimensionality reduction technique since it attempts to maximize between-class scatter while minimizing within-class scatter in the projected subspace. However, LDA has difficulty doing this if the original data are heteroscedastic or multimodal. Furthermore, the size of the LDA-reduced subspace is limited to c-1 (where c is the number of classes).
Local Fisher's discriminant analysis (LFDA) combines the strategy of LDA with a locality-preserving projection to provide a linear manifold learning technique that preserves the withinclass structure of the original space in the projected subspace; details of the LFDA algorithm applied in this study are provided in Sugiyama (2007). Briefly, LFDA seeks to find a transformation that preserves local neighborhood information, thereby ensuring that the underlying structure of the data distribution is preserved in the lower dimensional (size r) subspace. To accomplish this, the scatter matrices typical of LDA are scaled using an affinity matrix that measures the closeness of any two points relative to their k nnnearest neighbor. The parameters k nn and r must be optimized in concert with the classifier for each subject. LFDA has been previously deployed as a preprocessing step for classification of walking intention  and classification of expressive movement (Cruz-Garza et al., 2014) from EEG. A similar locality preserving projection was also employed for detection of ankle movement intention from low frequency EEG (Xu et al., 2014).

Classification algorithm
Once a suitable algorithm for dimensionality reduction was determined, we next identified a classification scheme to decode movement intent from our EEG-based features. Gaussian mixture model (GMM) classifiers are common in the fields of www.frontiersin.org November 2014 | Volume 8 | Article 376 | 5 biometrics and biomedical engineering because GMMs are capable of representing arbitrary statistical distributions as a weighted summation of multiple Gaussian distributions, termed components (Paalanen et al., 2006). Utilizing a GMM to compute the class-conditional probabilities in a maximum-likelihood classifier could improve performance over the traditional formulation, especially when the within-class feature set may be non-Gaussian, as could be the case for the temporally and spatially diverse EEG based features used in this study. The probability density function for a given training data set in the LFDA projected subspace, where K is the number of components, α k is the mixing weight, μ k is the mean vector, and k is the covariance matrix of the k-th component. The parameters of each GMM component K, including α k , μ k , and k , are estimated as those which maximize the log-likelihood of the training set given by: where p(x) is given in (1). Maximization of (3) is carried out using an iterative expectation-maximization (EM) algorithm (Vlassis and Likas, 2002), with the initial estimate of the parameters α k , μ k , and k established via k-means clustering (Su and Dy, 2007), until the log-likelihood reaches a predetermined threshold. The number of components K is a critical parameter for successful implementation of a GMM classifier. During training, we limited the maximum value of K to be 10 and computed the maximum log likelihood from (3) for each model with values of K ∈ {1, 2 . . . 10}. We estimated the optimal value of K as the model that minimized the Bayes information criterion, which has been reported as an effective measure for optimizing the number of GMM components (Li et al., 2012). In this manner, GMMs representing each movement class were specified for use in a maximum-likelihood classifier.
The parameters for each class-conditional GMM were computed using an optimization data set for each participant (see Classifier optimization section). The parameter space which must be explored in order to fit these mixture models can be quite large, especially if the feature dimension is large. Given the limited time and training data available during EEG studies, this learning task may be impractical, but as indicated in the previous section, LFDA has been shown to effectively reduce data dimensionality while preserving the statistical information. Thus, we applied LFDA dimensionality reduction on our EEG feature set prior to training and testing a GMM model for use in a maximum-likelihood classifier of intended motion.

Classifier optimization
The EEG feature matrix from each trial was split into two mutually exclusive sets: one for LFDA-GMM classifier optimization and one for classifier evaluation (Figure 1). The optimization data set was selected randomly from the full data set, and it comprised 400 samples (2 s) of data from each class. The optimization data set was then split into two equally sized exclusive subsets, one for training and one for testing. The parameters for the LFDA-GMM classifier (the nearest neighbor (k nn ) used in the affinity matrix, the dimensionality (r) of the projected subspace, and the number of mixture components (K) in the mixture model) were optimized for each subject and trial typeself initiated and triggered-using the optimization data set. Optimization involved three steps (Figure 1): (i) dimensionality reduction using LFDA for values of k nn and r from 1 to 249 and 1 to 250, (ii) identification of the optimal value of K for each class at each grid point in (i) using the training data from the optimization set, and (iii) computation of the accuracy of the LFDA-GMM classifier at each grid point in (i) using the testing data from the optimization set. The optimal parameters {k nn , r, K} for each subject were selected as those which produced the highest overall classification accuracy from the testing data.

Classifier performance via cross validation
The performance of the LFDA-GMM classifier with the optimal parameter set was analyzed for each subject and trial using repeated random sub-sampling cross validation (Figure 1). Repeated sub-sampling was chosen because the variable timing of the movements in each trial would result in an unequal number of samples from each class if k-fold cross validation scheme was used. The evaluation data set was randomly split into mutually exclusive training and testing data sets (Figure 1). Each of the three classes in the training set contained 600 data points representing 20% of the sit and stand classes. (Because the sit and stand classes were composed of ten 1.5 s long pre-movement epochs for each subject, their size was always equal). After training, LFDA-GMM classifier performance was analyzed using the testing data set, which contained all remaining data from the sit and stand classes, and an equal number of data points randomly selected from the quiet class. Thus, each class in the testing set contained 1900 data points. This test set structure was used to control for effects of class population size by assuring an equal number of testing samples in each class. During testing a classification decision was made for each data point, which represented a single time sample from the trial. The posterior probability of each data point was computed using the optimized GMM for each class and the data point was then assigned to the class that returned the largest value. This process yielded a classification decision for 1900 data points per trial. To avoid training bias, the random training and testing process was repeated 20 times and the average classification accuracies were reported for each subject under each condition (self-initiated and triggered movements). We performed post-hoc statistical comparisons between conditions using the non-parametric Kruskal-Wallis one-way analysis of variance.
To examine the effects of the ASR algorithm and the potential contribution of motion artifacts, we repeated the optimization and cross validation procedure using EEG data from premovement epochs pre-processed in the same manner as Figure 1 except that the ASR process was omitted. We also examined the classification accuracy using EEG epoched from movement onset to 1.5 s after movement onset both with and without the ASR algorithm. Finally, we divided the scalp into four major regions of interest (ROI) to assess the classification ability of each area individually. The ROIs included the frontal cortex (F3,F1, Fz, F2, F4, FC2, FC1, FC2, and FC4), the motor strip (C5, C3, C1, Cz, C2, C4, and C6), the parietal cortex (CP5, CP3, CP1, CPz, CP2, CP4, CP6, P3, P1, Pz, P2, and P4) and the central midline (FC1, FC2, C1, Cz, C2, CP1, CPz, and CP2). For each condition, we assessed within subject differences in accuracy across ROIs using the nonparametric Friedman test. The statistical sign test was used to assess if the difference in accuracy between self-initiated and triggered movements for each participant and ROI were significantly different from a distribution with a median of zero.

Demonstration of simulated real-time classification
We implemented a two-fold approach to demonstrate LFDA-GMM classifier performance in a simulated real-time environment using EEG data from the self-paced trial. The classifier was trained using ASR-cleaned EEG data from the first half of the trial with the optimal parameter set for each subject. Unlike during the cross-validation procedure, the time periods immediately following the movement execution were not trimmed from the data set but instead were included in the quiet class. Data from the second half of the trial, containing five transitions each of stand-to-sit and sit-to-stand, was used to test the controller in a simulated real-time manner resulting in a continuous time series of classification decisions.

OBSERVATIONAL EEG MEASURES
In addition to classification of movement intent, we computed several observational measures to help assess differences in cortical activity across the experimental conditions. We computed the MRPs from each subject during both the self-initiated and triggered conditions. To compute MRPs, each EEG channel was band pass filtered between 0.1 and 50 Hz and epoched from 2.5 s before movement onset to 1 s after onset. Each channel and epoch was baseline corrected using the mean voltage from 2.5 to 2 s before onset. Each channel was then averaged over all 20 epochs for each condition.
To ascertain differences between periods of quiet (i.e., rest between movements), pre-movement, and post-movement epochs under each condition (self-initiated and triggered) we computed the power spectral density (PSD) for each EEG channel with a frequency resolution of 0.12 Hz using the Thompson Multitaper method in Matlab with a time bandwidth product of 4. The PSD was computed after artifact removal with ASR but before band-pass filtering and standardization. EEG was common average referenced for purposes of PSD computation. The spatial distribution of alpha band (8-13 Hz) ERD was computed for the pre-movement and post-movement epochs under both conditions as was the change in power in the delta band (0.1-4 Hz). The change in power for both frequency bands was computed relative to the quiet epochs for each condition (self-initiated and triggered). We assessed statistical differences across conditions using the non-parametric Kruskal-Wallis one-way analysis of variance with a Bonferroni correction for multiple comparisons.

OBSERVATIONAL MEASURES
Standardized EEG and the linear envelope of EMG recorded during a typical trial for one subject is shown in Figure 2. EEG with and without ASR is shown, demonstrating the removal of high amplitude artifacts, especially in the time periods following movement onset. Although all 64 channels of EEG are displayed, those channels marked with an asterisk ( * ) were removed prior to classification of movement intention. The EEG PSD computed during rest (quiet standing) and the pre-movement epochs during the self-initiated and triggered trials is shown in Figure 3. The grand mean PSD across all participants and electrodes used for classification (lower inset, Figure 3) is shown. Two identifiable peaks are present in the rest condition, during which the subject was standing quietly; one in the theta band at approximately 7 Hz and one in the alpha band at approximately 11 Hz. Power in these bands were significantly greater at rest than during the premovement epochs under both conditions (p < 0.01 for both). Notably, the delta band power during the pre-movement epochs was greater than rest while the power in the theta and alpha band was greater during rest (upper inset, Figure 3). In the premovement epochs, there was significantly less power in the theta band (4-8 Hz) during self-initiated transitions compared to triggered (p = 0.004), while power in the alpha band (8-13 Hz) was not statistically different between conditions (p = 0.107). Finally, power roll-off, indicated by the slope of the PSD, was diminished in theta and alpha bands compared to surrounding delta and beta bands for the self-initiated pre-movement; however, roll-off was only decreased in the alpha band for the triggered condition.
The change in delta and alpha band power for the pre-and post-movement epochs, relative to the periods of quiet sitting and standing between movement executions, averaged over all participants is shown in Figure 4. In the delta band, we observed slightly increased power in the pre-movement epochs over all electrodes for both conditions, with slightly more delta power present in the self-initiated trials. In contrast, delta band power during the post-movement epochs was much larger, especially for the triggered trials, which showed nearly double the delta band power of the rest condition. The same level of increase was not observed over the full scalp in the self-initiated trials, although delta band power over the central midline electrodes increased by nearly 100%. Alpha band power was similar to quiet periods across most electrodes (note the difference in scale between alpha and delta power in Figure 4). Bilateral alpha band ERD was observed in both conditions; however for the triggered trials the ERD was less prominent and restricted to the central sensorimotor and parietal electrodes, while frontal and peripheral electrodes showed a slight increase in alpha power. Conversely, alpha ERD was stronger in the self-initiated condition, especially in the central-parietal areas of the scalp.
We found the presence of MRPs to be variable across subjects and conditions. In 3 subjects, MRPs were prominent across the scalp during the self-initiated movement epochs but not during  the triggered movements ( Figure 5A). For the remaining subjects, less prominent MRPs were present at some electrodes for both conditions ( Figure 5B). We examine the relationship between MRP and classification accuracy in more detail below.

CLASSIFIER VALIDATION
The LFDA-GMM classification accuracy surface followed a similar pattern for most subjects (Figure 6), rising sharply as the size of the reduced subspace (r) increased. Accuracy typically

FIGURE 5 | Example of movement related potentials (MRPs) recorded in two different subjects. (A) MRP from S5 indicating a difference between triggered (black line) and self-initiated (red line) movements. (B) MRP from S9
indicating similar, less prominent RPs for both the triggered and self-initiated trials. For each subject, MRPs were averaged across all 20 movements for each condition; movement onset is at 0 s.

FIGURE 6 | Example of a subject-specific accuracy surface created during LFDA-GMM classifier optimization.
The accuracy plotted at each point {r, k nn } on the surface is the average accuracy with the optimal number of mixture components (K ) for each class at that point.
peaked for r values between 50 and 125 before decreasing slightly, and then reaching a plateau as the value of r was further increased. Classification accuracy was generally insensitive to the k nn parameter with the exception of very low r values. The optimal parameter set for each subject and condition is provided in Table 1. Across subjects and conditions, the average dimension of the EEG-based feature space following LFDA was 88 (range 30-118), representing a significant reduction from the original size of 308. With few exceptions, the optimal accuracy was achieved using only one mixture component (K = 1) and thus, the LFDA-reduced EEG features were generally not strongly multimodal.
The mean overall classification accuracy obtained from the 20 times cross validation procedure for each subject and condition is shown in Figure 7 along with the overall mean across all subjects for each condition. The mean accuracy across subjects was 74.1 ± 5.7% for the triggered condition and 78.0 ± 2.6% for self-initiated. Testing sample size was equal across the three classes (1900 samples per class for each subject and condition). Interestingly, there was no significant difference in overall accuracy between self-initiated and triggered movements across the entire group of subjects. For subjects S2, S4, S5, and S7 decoding accuracy was significantly greater (p < 0.01) for the self-initiated sit-to-stand and stand-to-sit transitions compared to the triggered paradigm. Two subjects, S1 and S3, showed significantly better classification accuracy for the triggered movements compared to self-initiated, though with less strength (p < 0.05). The normalized confusion matrix for each condition was computed by summing the total number of predicted samples for each class across all 10 subjects and then dividing each predicted sum by the actual class sample size (Figure 8). We also computed the overall kappa coefficient (Cohen, 1968;Carletta, 1996) for each condition, resulting in κ = 0.61 for triggered and κ = 0.67 for self-initiated. For both triggered and self-initiated conditions, the quiet class was decoded with the highest accuracy and misclassifications for the quiet class were evenly distributed between the two types of movement (sit and stand). Notably, classification accuracy for all three classes was slightly, though not significantly, higher during the self-initiated trials. The majority of misclassifications for sit and stand movements were in the quiet class regardless of condition. Classifier confusion between movement types was slightly larger for the triggered paradigm, with 10.2% of sit movements misclassified as stand (as opposed to 4.2% for self-initiated) and 7.6% of stand movements misclassified as sit (compared to 3.0% for self-initiated).

FIGURE 7 | Mean accuracy (n = 20) by subject for decoding triggered and self-initiated sitting and standing from pre-movement EEG.
Error bars indicate ±1 standard deviation. Statistically significant within subject differences across conditions are indicated as follows: * p < 0.01, * * p < 0.05.
To assess the relationship between classifier accuracy and MRPs we computed the grand median area under the MRP curve for each condition and subject in a three step process. We first computed the area under the MRP of each channel for each movement epoch; a negative number for this area indicated a larger MRP presence. Next, we computed the median area under the curve for each electrode, and then we took the grand median area across all electrodes. We plotted this value against the mean classification accuracy for both the self-initiated and triggered conditions ( Figure 9A). Surprisingly, we did not find a strong correlation between area under the MRP curve and classification accuracy (R 2 = 0.09). Based on our prior observation that some subjects showed more prominent MRPs during the self-initiated movement compared to triggered, we computed the individual change in accuracy and the change in median area under the MRP curve across these conditions for each subject (Figure 9B). There was a slightly stronger, but still modest (R 2 = 0.27) correlation between individual change in accuracy and area under the MRP curve. Interestingly, the subject with the most visually prominent difference in MRP between conditions (S5, Figure 5A; blue arrow in Figure 9B) showed the second largest increase in accuracy between the self-initiated and triggered conditions. However, the subject with the largest increase in accuracy across conditions (S8, red arrow in Figure 9B) showed only a moderate increase area under the MRP curve. The two subjects with significantly greater accuracy for the triggered condition also had larger areas under the MRP curve in that condition ( Figure 9B).

CLASSIFICATION BY ROI
The mean and subject specific classification accuracy was lower for all four ROIs than with the full set of non-peripheral electrodes for both self-initiated and triggered movements (Figure 10), a result that was expected due to the lower number of electrodes used for classification. Of note, however, was that despite the differing number of electrodes within each ROI we observed few within subject significant differences in accuracy for each condition (Figures 10B,C). Similarly, when accuracy was averaged across the 10 subjects, there were no statistically significant differences between the ROIs for either condition. To assess the effect of self-initiated vs. triggered movements, we computed the within subject difference in accuracy for each ROI between these conditions ( Figure 10D). A majority of participants (8/10) showed similar or significantly greater accuracy for all four ROIs in the self-initiated condition. The two subjects (S1 and S3) who showed significantly greater accuracy for the triggered movements with the full set of electrodes also showed greater accuracy in several, but not all, ROIs in this condition. Interestingly, when the difference was averaged across subjects, only the motor strip ROI showed significantly increased classification accuracy for the self-initiated condition. Indeed, decoding accuracy of movement intent during self-initiated sitting and standing using the motor ROI was significantly greater than during triggered movement in 7/10 subjects, similar in 2/10 subjects, and decreased in only 1/10 subjects.

EFFECTS OF ARTIFACT REMOVAL
To examine the effect of the ASR artifact rejection algorithm, and the potential effect of motion or other artifacts on classification accuracy, we repeated the classifier optimization and cross-validation procedure for the self-initiated condition using three control data sets and compared those with the original preprocessing (Figure 11). The original data set is termed ASR pre in Figure 11. The first control data set was composed of the same pre-movement epochs consisting of 1.5 s of EEG data recorded immediately prior to movement onset, however, ASR was omitted from the pre-processing (Figure 1); this data set is termed Raw pre . We decoded movement intent using an equally sized epoch encompassing the 1.5 s time period immediately after movement onset. We processed these data with (ASR move ) and without (Raw move ) the ASR artifact rejection algorithm. We found that the ASR algorithm had no statistically significant affect on accuracy when using the pre-movement epochs to decode movement intent (Figure 11). This result was consistent for every subject and when accuracy was averaged across all subjects. When movement type was classified with EEG from epochs immediately after movement onset, a statistically significant increase in accuracy was observed in every subject when the data were not cleaned with ASR (Raw move ). Application of the ASR algorithm (ASR move ) resulted in a statistically significant drop in accuracy for decoding with the post-movement epochs in 9/10 subjects. When averaged across participants, no significant difference in accuracy was observed between ASR cleaned pre-and post-movement epochs, while accuracy was significantly higher for decoding with raw post-movement data.

SIMULATED REAL-TIME CLASSIFICATION
The results of simulated real-time decoding using cleaned EEG data are shown in Figure 12. Class-wise accuracy in this demonstration was different than observed from the cross-validation (Figure 8) an effect caused by the training sample bias inherent to the two-fold procedure used for the demonstration. The quiet class (0) contains a larger number of samples than either stand-to-sit (class 1) or sit-to-stand (class 2) resulting in very high accuracies during quiet periods. Confusion between classes 1 and 2 was present during most transitions; the low number of transitions used in this demonstration likely contributed to this confusion. Errors at the beginning and end of the movement periods skewed toward class 0 (quiet).

CLASSIFICATION OF SELF-INITIATED AND TRIGGERED MOVEMENT FROM PRE-MOVEMENT EEG
Our results demonstrate successful, high accuracy classification of movement intent in healthy individuals from delta-band EEG recorded before movement execution. We framed our experiment into a three-class problem where each time point was classified into one of three states: quiet, stand-to-sit transition, or sit-tostand transition. It is important to note that we trimmed the The mean difference in pre-movement decoding accuracy between the self-initiated and triggered conditions for each subject ±1 standard deviation. Asterisks ( * ) indicate differences which were statistically significant (p < 0.05) from a distribution with a median of zero based on the sign test.
time periods of actual movement execution-as determined from EMG activity-from our EEG recordings. Thus, our classifier was trained and tested using mutually exclusive EEG datasets recorded during either quiet standing or quiet sitting but when subjects presumably were preparing for the incoming action. We labeled each time point in the 1.5 s epoch before movement onset according to the type of movement that was executed in the future: stand-to-sit or sit-to-stand. All other time points were placed into a single quiet class. Classification ability was assessed in two different movement execution paradigms, one that was cued by an audio signal (triggered) and one that was self-paced (selfinitiated). Interestingly, we observed no statistically significant difference in classification accuracy between these two conditions, though average accuracy across the 10 subjects was slightly higher for the self-initiated condition (78.0 ± 2.6%) compared to triggered (74.7 ± 5.7%) and both of these were significantly better than chance accuracy of 33.3%. Prominent MRPs were not visible in all subjects ( Figure 5) and we found almost no correlation between median area under the FIGURE 11 | Classification accuracy using pre-and post-movement epochs with and without ASR pre-processing. The classifier was trained and tested for the self-initiated case using pre-movement epochs with the original pre-processing pipeline (ASR pre , green) and using pre-movement epochs omitting ASR from pre-processing (Raw pre , red). As a control, the classifier was also trained and tested using equally sized epochs (1.5 s) immediately following movement onset that were pre-processed with (ASR move , gray) and without (Raw move ) ASR for artifact rejection.
MRP curve and classification accuracy ( Figure 9A). For within subject comparisons between conditions, we observed significantly better accuracy in four of ten subjects during the selfinitiated compared to triggered paradigm, while two subjects had higher accuracy for triggered standing and sitting. When examining subject specific changes in accuracy across the two different paradigms, we found a slightly stronger correlation between increased accuracy and area under the MRP curve. And the two individuals that showed a decrease in accuracy in the selfinitiated vs. triggered trials also showed an increased area under MRP curve, indicating less prominent MRPs. These results appear to contradict previous examples which indicated that MRPs may be more prominent in self-paced vs. cued movement paradigms (Jahanshahi et al., 1995;Jankelowitz and Colebatch, 2002;Cui and MacKinnon, 2009). There are several possible explanations. First, our experimental paradigm included a relatively low number of epochs (n = 20) for each condition, compared to traditional studies of MRPs which typically utilize close to 100 (Shibasaki and Hallett, 2006). This low number of epochs may be the reason for the large variability in the presence of MRPs ( Figure 5). Additionally, in the self-paced experiment, participants were instructed to pause 3-10 s between each movement though they were also instructed not to count the seconds between each movement. As a result, participants rarely waited 10 s between self-paced movements; most periods of quiet lasted 5 s or less. Previous studies have observed trial-to-trial variation in timing and power of MRPs relating to self-paced left and right hand movements, making classification of those movements using low frequency features more difficult . Another study found that while they were present for most-but not allsubjects and movements, low frequency features were less critical than ERD/ERS in classifying four different types of movement from EEG (Morash et al., 2008). The latter study utilized the contingent negative variation (CNV), which is a low frequency,

Frontiers in Neuroscience | Neuroprosthetics
November 2014 | Volume 8 | Article 376 | 14 event related-potential entailing a widespread negative shift in EEG observed in paradigms involving conditional and imperative stimuli (Walter et al., 1964). While our paradigm did not involve dual stimuli, it is possible that some participants experienced a similar effect due to the alternating nature of the movements. That is, completing the previous maneuver (sitting or standing) may have created a conditional response in which the subject then began to prepare for the next movement, which would be the opposite of the prior one. This conditional response may be another reason that we did not observe prominent MRPs in some subjects. Indeed, trial-to-trial variation in CNV amplitude has been described previously and this variation may be representative of anticipated events and/or fluctuations in attention to the task (Scheibe et al., 2010). The observed variation in MRPs may also be responsible for the skewed misclassification of sit and stand movement intentions as quiet (Figure 8). Note that while the full time series of EEG data contained more samples in the "quiet" class than in the "sit" and "stand" class, an equal amount of data from each class was used for cross-validation, and thus, this pattern of misclassification was not a result of training bias. Variable timing of movement execution and conditional response may have affected the prominence of MRPs, but it did not hinder classification accuracy. One reason for this may be the time-embedding of our classification features which encompassed information from up to 50 ms before the current time point, helping to alleviate previously reported MRP-based feature variability . Low frequency EEG has been shown to contain information regarding intention (Lew et al., 2012), direction (Liao et al., 2007;Vuckovic and Sepulveda, 2008;Waldert et al., 2008;Robinson et al., 2013), velocity (Bradberry et al., 2010), and type (Agashe and Contreras-Vidal, 2013) of hand movement. In the lower extremity, the ability to detect voluntary ankle dorsiflexion movement from MRPs with accuracies up to 80% has been reported (Niazi et al., 2011;Xu et al., 2014). During walking, intra-stride changes in electrocortical activity coupled to gait phase have been observed at frequencies as low as 3 Hz (Gwin et al., 2011) and inter-limb and intra-limb kinematics (Presacco et al., 2011(Presacco et al., , 2012 as well as the intention to start and stop walking  have been decoded using delta band EEG. In another recent study, features extracted from the delta band were the most heavily weighted for single trial classification of walking movement intention from EEG recorded prior to movement (Velu and de Sa, 2013). Our results, which classified lower extremity movement type using pre-movement EEG, corroborate these findings and provide further evidence that low frequency EEG contains discriminative information pertaining to lower extremity movement intent.

CLASSIFICATION BY REGION OF INTEREST
The results from our ROI analysis (Figure 10) support the hypothesis that stand-to-sit and sit-to-stand transitions are preceded by event-related activity across a distributed, sparse cortical network. As expected due to the reduced number of electrodes, no ROI reached the classification accuracy attained when all electrodes were included in the classifier. When averaged across subjects, there were no statistically significant differences in classification accuracy between the ROIs for either condition, despite the difference in number of electrodes. The ROI analysis also revealed a statistically significant increase in accuracy for within subject differences across conditions (self-initiated vs. triggered) when using only the electrodes over the motor area. A similar difference was not found for any other ROI or for the entire scalp. This result suggests that the primary motor cortex (M1) region contains more discriminative information for identification of standing and sitting intention when the movements are self-initiated compared to cued. This finding is supported by previous work indicating MRPs from this region differ when the motor task emphasized sequence initiation compared to rhythm (Bortoletto et al., 2011). EEG recorded from these electrodes has also been demonstrated to most accurately track movement initiation using other frequency bands such as mu/alpha ERD and beta ERS (Wolpaw et al., 2002).

ARTIFACT SUBSPACE RECONSTRUCTION
This study, along with previously mentioned work, establishes compelling evidence for neural correlates of movement within EEG signals recorded immediately prior to movement execution; however, it is important to address the possible role of artifacts, both physiological such as muscle and eye and non-physiological, such as movement. Our signal processing approach for classifier training and evaluation (Figure 1) was designed to minimize the effect of artifacts in several ways. First, we eliminated frontal, temporal, and occipital electrodes which can be contaminated by EMG and/or EOG artifacts. Second, we trimmed all EEG that was recorded during periods of movement as indicated by lower extremity EMG from our data set, leaving only EEG recorded during periods of quiet sitting and standing for classification. Third, we applied a PCA-based artifact rejection algorithm (ASR) that was designed to eliminate high amplitude and high variance artifacts, such as those from movement or muscle, from EEG (Mullen et al., 2013). Our pre-processing analysis demonstrated similar power spectral density between rest (quiet standing) and premovement periods under both conditions (Figure 3), suggesting that our pre-processing steps were effective in removing artifacts from EEG. We also observed alpha ERDs in the period immediately following movement onset (Figure 4), especially during self-initiated trials, an observation that would have been unlikely if muscle activity had remained in the cleaned-EEG signals since EMG tends to have power in this frequency band.
To further elucidate the possible role of artifacts and these steps to mitigate them, we compared the LFDA-GMM classifier performance when it was trained and tested with three different control data sets with our original processing pipeline (Figure 11). This analysis showed no statistically significant difference in accuracy, regardless of whether the pre-movement EEG was cleaned with ASR or not, suggesting that artifacts were not present and therefore did not affect classification using the pre-movement epochs. We did observe a significant increase in accuracy when the premovement epochs were replaced with equally sized epochs immediately following movement onset that had not been cleaned using ASR. After ASR cleaning, classification accuracy was commensurate with pre-movement epochs, although with a slightly larger standard deviation across subjects. The increased accuracy using post-movement epochs without ASR suggests that artifacts may have been present during this time and these artifacts may have enhanced decoding accuracy. The decreased accuracy following ASR suggests that this algorithm is effective at removing high amplitude artifacts from EEG data. This conclusion is further supported by the simulated real-time demonstration using ASR-cleaned data. The time periods after movement onset were included in the quiet class during training and were decoded with high accuracy during testing (Figure 12). But, caution should be exercised regarding the conclusion that ASR completely eliminates low frequency, high amplitude artifacts. We note that while we did observe alpha ERD in ASR-cleaned post-movement epochs, we also observed enhanced power in the delta band across the scalp, particularly in the triggered condition (Figure 4). One possible explanation for the post-movement increase in delta band power in the triggered trials could be residual head movement and/or muscle artifacts as the participant reacted to the audio cue to stand or sit. Further spectral, topographical, and temporal analysis should be undertaken to parse movement related artifacts from true electrocortical sources recording during the actual sitting and standing movements. In particular, the parameters of the ASR algorithm can be optimized to more aggressively remove artifacts at the expense of potentially removing true EEG. We emphasize that our primary analysis involved only EEG from pre-movement and quiet periods, thereby limiting the contribution of these potential artifactual components as indicated by the above analysis.

EEG USE IN REHABILITATION AND RESTORATION OF MOVEMENT
To our knowledge, this is the first study that classifies this type of gross, full lower extremity movement intention-sit-down, stand-up, or quiet-from non-invasive EEG signals. Previously, surface EMG from leg muscles has been used with an LDA classifier to identify standing and sitting transition in amputees with accuracies greater than 99% (Zhang et al., 2012). Achievement of these high accuracies required the use of a post-processing majority voting step, which resulted in a decision delay of up to 400 ms. Another approach has deployed center of pressure to detect sitting and standing transition in individuals with paraplegia (Quintero et al., 2011). Classification of sitting and standing using EEG offers advantages over these approaches. On average, we were able to achieve 78% accuracy using features extracted from the pre-movement epochs with no post-processing required, thereby minimizing delay between movement intention and classification. It should be noted that our classification accuracy was assessed using single time points that were randomly selected from each trial. This conservative approach was necessary to prevent model over-fitting during training and to assure an equal number of data points in each class during testing due to the relatively low number of movements executed (20 per condition) for each subject. An example of the LFDA-GMM algorithm in a simulated realtime environment is shown in Figure 12. We note that classifier training was not optimal for this demonstration; only 5 stand-tosit and sit-to-stand transitions were employed. Further, clinical deployment of the classifier as a component of a BMI could be significantly improved by addition of an aggregate post-processing step-such as requiring a number of consecutive time points to be predicted as the same movement type or a sliding window moving average with a threshold-to trigger a change in state. The parameters of this post-processing step need to be tuned for each subject and application to maximize accuracy and minimize false positives. Future studies will investigate this possibility and the tradeoff between gains in accuracy and increased classification latency from post-processing.
One drawback of utilizing GMM based classifiers is the size of the parameter space which must be learned, which is given by K * (1 + d * (d − 1)/2) + K * d, where K is the number of Gaussian components in the mixture, and d is the dimensionality of the data to be fit (Li et al., 2012). To fit a GMM to our timeembedded EEG-based feature data set, which includes data from 28 channels of EEG at 11 time points and a maximum of K = 10 components for a given class, requires learning a parameter space of dimension 4.76 × 10 5 . Our results demonstrate that LFDA is a powerful dimensionality reduction technique; the median dimension of the reduced subspace was 96 (Table 1), representing a median reduction of 69% across subjects. LFDA reduced the size of the GMM parameter by an order of magnitude, resulting in a large decrease of computation time to fit the models of the classifier. Classifier optimization and training was performed using custom software developed in Matlab®, including the parallel processing toolbox, run on a dual core PC (2.40 GHz, 24 GB RAM). On average, optimization across the full LFDA-GMM parameter space was complete in less than 15 min per subject, and training of the optimized LFDA-GMM classifier in less than 5. If deployed for control of an assistive device, LFDA-GMM classifier optimization and training may be required before each session of use; these results suggest this is feasible. Examination of the optimization surface (Figure 6) shows that gains in accuracy level-off at moderate values of r while accuracy is relatively insensitive to k nn . The same trend is observed in all subjects, with some showing decreases in accuracy for increasing r-values, while in others there is no difference in accuracy as the parameter values are increased. Thus, these parameters could be limited to smaller values, thereby reducing the parameter space to be searched during LFDA-GMM optimization. However, the optimal parameter set is expected to vary with the task and also with the ability of the subject to learn how to operate the BMI over time, and so caution should be exercised when determining the upper limits. Also, full covariance matrices ( k ) were deployed for each component of the GMMs; however, if the subspace of the data following LFDA dimensionality reduction was large, employing diagonal covariance matrices could be used as a way to speed classifier training.
The LFDA-GMM classifier presented here could be incorporated into a closed loop BMI system with an exoskeleton to restore function to individuals with paralysis. Such a system would be comprised of a shared control paradigm, whereby the gross motor instruction (in this case, the intention to sit-down or stand-up) is extracted from the user's EEG and the commands to execute the movement are performed autonomously by the exoskeleton. In this setup, the exoskeleton would be triggered at the first time point in which the BMI detected a change in class; a process that would likely include a post-processing step requiring a sequence of consistent classifier decisions to trigger a change in state. The decoding algorithm would then be blanked so that no state changes could be triggered during the execution of a movement.
Our observed accuracy of 78% in self-paced movements would need to be improved for clinical viability. However, the data used in this study were purely observational, while operation of a BMI is a learned skill that incorporates feedback to the user regarding performance; thus accuracy of the BMI may increase as the user gains additional experience with the device. In the future, EEG and EMG could be combined to create a comprehensive neural-machine interface for control of advanced prosthetics. The combined EEG-EMG interface could provide intuitive control of artificial limbs while minimizing delay between detection of voluntary movement intention and its execution. Our classification approach could also be used in an intervention to treat phantom limb pain, whereby a descending motor command is determined from EEG and a motorized prosthesis executes the movement providing afferent feedback which could obviate maladaptive cortical reorganization following amputation. EEG-based classification of movement intent could also be incorporated into a neurorehabilitation protocol to recover more normal motor function in individuals with neurologic impairments. For example, the EEG based classifier would activate a device to assist movement, thereby creating more normal afferent feedback, which could enhance brain plasticity and speed motor recovery (Daly and Wolpaw, 2008). Such a strategy requires extraction of motion intent from the motor impaired population; in this study only healthy able-bodied individuals were tested. Future studies will examine the ability to apply LFDA-GMM classification to individuals with central nervous systems deficits with an aim toward neurorehabilitation strategies.