A brain–computer interface for the continuous, real-time monitoring of working memory load in real-world environments

Abstract

We developed a brain–computer interface (BCI) able to continuously monitor working memory (WM) load in real-time (considering the last 2.5 s of brain activity). The BCI is based on biomarkers derived from spectral properties of non-invasive electroencephalography (EEG), subsequently classified by a linear discriminant analysis classifier. The BCI was trained on a visual WM task, tested in a real-time visual WM task, and further validated in a real-time cross task (mental arithmetic). Throughout each trial of the cross task, subjects were given real or sham feedback about their WM load. At the end of the trial, subjects were asked whether the feedback provided was real or sham. The high rate of correct answers provided by the subjects validated not only the global behaviour of the WM-load feedback, but also its real-time dynamics. On average, subjects were able to provide a correct answer 82% of the time, with one subject having 100% accuracy. Possible cognitive and motor confounding factors were disentangled to support the claim that our EEG-based markers correspond indeed to WM.

Appendices

Appendix 1: Images used for Task 1

Figures used to determine the memory span

Geometric Shapes	Fruits	Landscape
Circle	Apple	House
Pentagon	Banana	Building
Square	Orange	Church
Rhombus	Pear	Castle
Cross	Grape	Bridge
Star	Watermelon	Tower
Triangle	Pineapple	Tree

Figures used for the test

Animals	Vehicles	Supplies	Clothes
Cat	Plane	Book	Trousers
Deer	Train	Scissors	Shirt
Dog	Skateboard	Pen	Hat
Elephant	Truck	Ruler	Shoes
Penguin	Car	Backpack	Socks
Snake	Ship	Compass	Belt
Turtle	Bicycle	Set square	Tie

Appendix 2: Technical aspects of the BCI system

Offline study: designing the BCI

The offline study was a feasibility study, in which the constrains that would be encountered online were taken into account and reproduced. Examples of the online constraints include low latency of the online system (and therefore only a low number of features can be afforded for a quick analysis), limited calibration data for short calibration sessions, and the inability to remove artifacts such as eye-blinks. We called these parameters $P_1, P_2, P_3, \ldots$, and an acceptable range of each of them was determined. For different combinations of the acceptable parameters, the following pipeline was performed:

• (1) 20 subjects performed Task 1 while wearing the EEG set. Frequencies below 1 Hz and above 45 Hz were removed from the EEG signal with a 3rd order Butterworth filter. The EEG data were segmented into epochs of $P_1$ seconds.

• (2) A subject $s_i$ was removed from the dataset. The data of subject $s_i$ were divided into two subsets, one subset with $P_2$ epochs that will simulate the calibration session (“Online tests: the calibration session” section), and another subset with the remaining epochs. The calibration data and the data from the remaining subjects were used as training data. The data from the subject that were not part of the calibration data were used as testing data.

• (3) For the training data, each epoch was visually inspected and all the epochs contaminated with noise or muscular artifacts were rejected. In particular, epochs with eye-blinks or with arousal flags were rejected. An arousal flag was placed on an epoch if ether a distracter or the target was displayed during its course. During preliminary tests, subjects had reported an arousal effect due to the appearance of distracters or targets. The testing data were not cleaned.

• (4) For both the training and the testing data, spectral features were extracted using the Matlab p-Welch function, with a Hamming window of 0.5 s. The spectral features were absolute and relative power in the following bands: delta (1–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), lower beta (12–20 Hz), upper beta (20–30 Hz) and lower gamma (30–45 Hz). The relative power in a band is the fraction of the total power in that band. The latter being normalised has the advantage of reducing inter-subject variability. Having 16 channels, 2 features per band, and 6 bands, we obtained 192 features for each epoch.

• (5) The calibration data were expanded by adding noisy copies. Enough noisy copies were created so that the number of epochs of the expanded calibration data divided by the number of epochs of the total training dataset equalled $P_3$. The noise added was Gaussian noise with zero mean and standard deviation equal to $P_4$ times the standard deviation of the feature.

• (6) To select relevant features, orthogonal forward regression (OFR) (Stoppiglia et al. 2003) was performed on the above matrix, the best $P_5$ features were kept. OFR is a linear regression technique that can be used as a supervised feature selection technique. In the first step of OFR, features are ranked in order of decreasing correlation to the classifier output; the first selected feature is the top-ranking feature. In the second step, all remaining features, and the output, are orthogonalized with respect to the first selected feature, thereby discarding the part of the output that was explained by that feature; the projected features are ranked in order of decreasing correlation to the projected output, and the top-ranking feature is selected. Orthogonalization, ranking and selection are iterated until $P_5$ features are selected. OFR was used instead of other common procedures such as Common Spatial Pattern (CSP) (Koles et al. 1990; Cheng et al. 2017) in order to minimise the number of required sensors for future implementations. In addition, the weighted average of features from different sensors is similar to the spatial average provided by CSP.

• (7) A linear discriminant analysis [LDA (Fisher 1936)] classifier was trained with the selected features and their corresponding labels (high or low-WM load). The output of the classifier is the posterior probability that the EEG segment belongs to the high-WM-load condition. A typically high-WM-load segment then would have a high posterior probability, therefore we defined the WMLE as the posterior probability.

• (8) We performed cross-validation by iterating over all possible subjects $s_i$ for a given set of possible parameters ($P_1$, $P_2$, $P_3$, $P_4$, $P_5$), the performance of the classifier measured as the AUC was computed as a function of said parameters. The set of parameters that maximised the performance of the classifier was chosen, in other words, the parameters were optimised by cross-validation. Cross-validation is a family of techniques to asses a model’s ability to generalise when faced with new data, in this case, the testing data that were removed from the dataset. The reader interested in the topic can consult (Kohavi 1995) for a more in-depth discussion.

The set of parameters that represented the best trade-off between classification performance and feasibility was chosen to build the actual BCI. The values are shown in Table 2.

It is important to notice that although optimising by cross-validation might overestimate the performance results, we did this only to design the classifier, and all the results presented in this work are online results. Online there was no hypothesis testing (parameter selection), therefore, as the BCI design (set of parameters) was fixed for the online tests, p value corrections are not necessary. Furthermore, even though the number of features was optimised by cross-validation, the set of features itself (step 6) was chosen inside every cross-validation fold (step 8).

Online analysis: from EEG recordings to a WMLE

The online experiments followed the procedure outlined in the previous section. As before, recordings started with a calibration session, calibration data were cleaned and then added to the cleaned data from the 20 subjects recorded offline. This time however the set of parameters ($P_1$, $P_2$, $P_3$, $P_4$, $P_5$) was already set, and the pipeline described in the previous section was performed with the values of the Table 2.

Figure 8 shows how the information from the previous section (offline study) is integrated with that in this section (calibration) to select a good set of features and train the classifier that was designed offline.

Fig. 8

Design methodology

Afterwards, the BCI is able to provide a WMLE. A continuous stream of data is analysed. A sliding window including the last 2.5 s of EEG is used as input.

Appendix 3: Arithmetic operations in Task 2

A random sequence of digits ${d_1, d_2, \ldots , d_n}$ was presented to the subjects in each trial. Three possible ways of manipulating the digits were suggested to the subjects, who were asked to choose the one that felt more resource demanding for them:

• Progressive multiplication. Multiply $d_1 d_2 \ldots d_i$ until time is over

• Pairwise multiplication and successive addition. Multiply $d_1$ and $d_2$ and store the result. Add the result to the product of $d_3$ and $d_4$, replace the result. Add the result to the product of $d_5$ and $d_6$, replace the result. Continue until time is over.

• Free choice. Subjects comfortable with their arithmetic skills were left to choose the structure of the operations, provided they maintained a high level of use of their mental resources.

Appendix 4: Estimation of statistical significance

We developed a method to estimate analytically the statistical significance of the performance of a two-class classifier. The null hypothesis is that the results come from a random classifier, whereas the alternative hypothesis is that the classifier is based on informative features. The first step towards estimating significance, then, is to choose a random classifier and to determine its success rate. Let us assume that our dataset consists of N examples, with $N_1$ examples of class 1 and $N-N_1$ examples of class 2, with $N_1 \ge N/2$. The best that a random classifier can do is to take into account the prior probabilities of the classes. Denoting by q the prior probability of class 1, assumed to be larger than 0.5, and estimated by $N_1/N$, a possible random classification rule is to assign any object to class 1 with probability q. The probability of correct classification of this classifier (which can be estimated by its rate of correct classification) is given by : $c_0 = q^2 + (1-q)^2$. Let us define a random variable whose realisation $z_i$, for example i, is

$$\begin{aligned} z_i = {\left\{ \begin{array}{ll} 1 &{}\text {if the random classifier classified example } i \text { correctly} \\ 0 &{}\text {otherwise} \end{array}\right. } \end{aligned}$$

The total number of successes, $Z = \sum _{i = 1}^N z_i$, follows a binomial distribution $Z \sim B(N,c_0)$, and hence the probability of obtaining exactly k successes is

$$\begin{aligned} Pr(Z = k) = \left( {\begin{array}{c}N\\ k\end{array}}\right) c_0^k(1-c_0)^{N-k} \end{aligned}$$

The p value is, by definition, the probability of obtaining results at least as extreme as the observed ones, assuming that the null hypothesis is true. Our goal is to compare a random classifier with a particular classifier that yielded c correct answers. In this case, “results as extreme” means observing at least c correct answers in a random classifier. Therefore, the p value associated with the null hypothesis defined above can be computed as

$$\begin{aligned} p = \sum _{k=c}^N Pr(Z = k) \end{aligned}$$

In general, the use of any other random classifier would lead to a different $c_0$. In particular, the most efficient classification rule under complete lack of informative predictors is the zero classifier, that assigns all the objects to the largest class. In the above notation, the rate of correct classification of a zero classifier is $c_{0z} = q$. As, by definition, $0.5< q < 1$, it is easy to show that $c_{0z} > c_0$ for all q. However, although the zero classifier is the best classification rule when no relevant predictors are available, for a zero classifier $Pr(Z = k) =0$ if $k\ne N_1$ (by definition a zero classifier can correctly predict only $N_1$ objects). Therefore, for any classifier with a number of correct predictions larger than $N_1$, p would be zero. It is a good practice to compare classification results with those of a zero classifier when facing imbalanced datasets. However, a zero classifier is not useful for computing p values.