Neural decoders can classify neural signal to characterize or predict responses from the brain, spinal cord, or peripheral nervous system. Their role is fundamental for building a brain-machine interfaces (BMI) which have clinical potential to assist or treat neurological or psychiatric patients.
In this post, we’ll be building an event-related potential (ERP) decoder that detects electroencephalogram (EEG) error-related potentials (ErrP) using Matlab. The goal of this ERP decoder is to characterize ErrP signals that are generated when wrong actions are perceived. With minor modifications, this decoder can predict whether the person has noticed an error or not in real-time. We will run the ERP and classification analysis based on the signal processing methods mentioned in Iwane et al., 20161 and the Analyzing Neural Time Series Data book2.
Source codes in the post are available upon request: minsuzhang@gmail.com
1. Data curation
- We will first load the EEG recording data and check if there is any problem with the structure of the data. The data was recorded at 500 Hz sampling rate from a subject while he was performing 4 blocks of the error monitoring task mentioned in Iwane et al., 2016. The data from the 4 blocks were concatenated and saved in the matrix size of 604314 x 10.
- Here, the number of rows (604314) is the number of timepoints. This corresponds to around 20 min of recorded data, which makes sense since each of the block takes a around 5 min. The number of columns is 10 which corresponds to numbers of EEG channels (8) + EOG channel (1) + event marker (1).
- Let’s also check the summary of the recorded event markers. Below table summarizes the event markers and the count of the markers. In this experiment, markers 9 and 13 were used to mark the onset of Correct trials while marker 10 was used to for the onset of Error trials. The error rate set for the experiment was 0.3 which is roughly the value of #(error trials) / #(total trials).
2. Temporal filtering
- We are then going bandpass filter EEG channels between 1 to 10 Hz in order to increase the signal-to-noise ratio of ErrP. The 8 EEG channels [1-10] Hz bandpass filtered signal is shown below. Each channel is separated by 50 μV for easier visualization.
- The signals from each electrode are not overlapping much, meaning that most of the amplitude is within 50 μV range, which is an indication for clean EEG data.
3. Epoch segmentation
- An epoch is a data structure that is in the 3-dimensional form of time samples x electrodes x trials (1400 x 8 x 404). We will extract epochs from our data by segmenting a piece of data in our time window (-1.3 s to 1.5 s) and concatenating the matrices of electrodes x trials.
- We will also save the labels (Correct vs. Error) of this epochs in an array. Let’s also save the block number the trials were extracted from in blockID variable to examine if there were any effect of blocks. While segmenting epochs, we will do trial-wise baseline normalization to remove task irrelevant attenuations of activity.
4. Trial rejection
- Let’s use automatic trial rejection to reject any trials that have absolute value of above our threshold (50 μV). Notice how only the number of trials in the epoch structure was reduced (404→400).
5. Butterfly plot
- A butterfly plot overlays ERP from all of the electrodes on a same figure. It is useful for detecting bad electrodes or electrodes contaminated with artifacts.
6. ERP image
- ERP image is a 2-dimensional representation of the EEG data from a single electrode (FCz). Single-trial EEG traces are stacked and color coded to show changes in amplitude as colors. They are useful for inspecting trial amplitudes and behaviors in the time domain.
7. Event-related potential (ERP)
- We will then average the ERP images with respect to their trial condition to obtain an ERP.
- The ERP above shows error-related negativity around 300 ms and error-related positivity around 500 ms which is in line with ErrP literature.
9. Feature engineering
- Feature extraction
- Time samples from all 8 channels.
- Power spectral density (PSD) between 2 to 10 Hz.
- Feature normalization
- Dimension reduction using principal component analysis (PCA)
- Classification using diagonal linear discriminant analysis (LDA)
10. Classification analysis
- For non-balanced data set like ErrP (30% Error + 70% Correct trials), we will used the area under the receiver operating characteristics (ROC) curve to evaluate the classification accuracy.
11. References