Enhancing EEG-based cross-day mental workload classification using periodic component of power spectrum

We conducted an experiment to collect EEG data and evaluated the day-to-day stability as well as the cross-day classification performance of different features. In addition, we validated the classification results of various features using a public dataset.

2.1.1. Self-collected dataset

The study involved 19 graduate students (12 females), who received financial compensation for their participation. Participants ranged in age from 21 to 27 years (M = 23.8 years, S.D. = 1.18). None of the individuals were pregnant, on medication, or reported known fatigue-related problems. Participants were asked not to consume alcohol the night before the experiment and to abstain from caffeine on the morning of the test. The study was conducted in accordance with the principles embodied in the Declaration of Helsinki. All participants provided their written informed consent before the experiment. The experimental procedure was approved by the Ethics Committee of Tianjin University.

The experimental task for this study utilized the Multi-Attribute Task Battery II (MATB-II) software (Comstock Jr and Arnegard 1992). Three subtasks, namely system monitoring, tracking, and resource management, were employed as components of the MATB. The difficulty of these three subtasks was adjusted to modulate MWL at two levels (low/high). All participants underwent MATB-II training before the formal experiment until performance reached asymptote with minimal error. In this study, participants were asked to complete two sessions (Day 1 and Day 2) separated by more than 48 h. Each session consisted of four randomly assigned blocks, including two blocks of low-load tasks and two blocks of high-load tasks. Each block included a 3-minute rest period followed by a 5-minute MATB task, as depicted in figure 1.

Zoom In Zoom Out Reset image size

The EEG signals were recorded at a sampling rate of 1000 Hz throughout the experiment using the NeuroScan NuAmps 64-channel amplifier. We recorded 60 channels of EEG data using the international 10-20 electrode system. Impedance levels were maintained below 5 kΩ to ensure accurate signal acquisition. The right mastoid electrode served as the online reference. All data processing and analyses were performed using EEGLAB Toolbox functions (Delorme and Makeig 2004) and additional scripts developed in MATLAB 2020b^® (The MathWorks, Natick, MA). The EEG data were re-referenced to the average of the left and right mastoids, band-pass filtered between 0.1 and 100 Hz, and downsampled to 250 Hz. Visual inspection was used to eliminate muscle artifacts, followed by an independent component analysis (Jung et al 2001) to remove blink and muscle artifacts.

2.1.2. Public dataset

The continuous EEG was split into 10 s segments with a 1 s overlap. PSDs were then calculated for each segment using the Welch method, which involved applying Hann-windows of 2 s with a 50% overlap. To parameterize the raw PSD, we used the FOOOF method (Donoghue et al 2020) to decompose the power spectrum into periodic and aperiodic components. Thus, we obtain the raw spectrum and the spectra of two components with a frequency resolution of 1 Hz. Then, the power of the two components and the original spectrum were averaged over the canonical frequency bands, specifically theta (4–7 Hz), alpha (8–13 Hz), low beta (14–20 Hz), high beta (21–30 Hz), and gamma (31–40 Hz) (Newson and Thiagarajan 2018), as shown in figure 2.

Zoom In Zoom Out Reset image size

We computed grand-averaged PSD brain topographic maps for three types of PSD features, including raw PSD, periodic components, and aperiodic components. These maps were averaged across all participants per condition, which consisted of three MWL levels and two days. In addition, we show the spectra of three types of features on different days and MWL states in typical channels FZ and OZ. After considering previous research and the statistical results of our study, we identified the region of interest for each frequency band associated with the MWL. We then calculated the mean PSD across channels for each condition in these regions of interest.

Two-way repeated measure ANOVAs (2 d × 3 MWLs) were used to examine changes in three types of PSD features across days for each channel. False discovery rate correction (Benjamini and Hochberg 1995) was applied to all p-values to account for multiple tests. In addition, we performed two-way repeated measures ANOVAs (2 d × 3 MWLs) followed by a post hoc paired t-test for the regions of interest associated with each frequency band. All statistical analyses were performed using the EEGLAB Matlab toolbox. We employed three different p-value thresholds: p < 0.05, p < 0.01, and p < 0.001. This approach allowed us to assess the degree of evidence against the null hypothesis at varying levels of rigor.

Next, we examined the test-retest reliability of the features to evaluate day-to-day stability. As defined by (Elliott et al 2020), test-retest reliability refers to the reproducibility of results from a given measure under nearly equivalent conditions. A common statistic used to assess test-retest reliability is the ICC, typically defined as the ratio of the variance of interest to the sum of the variance of interest plus error (Shrout and Fleiss 1979, McGraw and Wong 1996). We calculated the ICC to assess the sensor-level day-to-day stability of three types of features (raw PSD, periodic components, and aperiodic components) of each frequency band (60 channels). In other words, ICC calculations were implemented for each channel in the PSD brain topography to obtain spatially specific reliability estimates. Furthermore, we calculated the ICC and their 95% confidence intervals of the averaged PSD across channels for the region of interest in each frequency band. Additionally, (Koo and Li 2016) recommend using a two-way combined-effects model and an absolute agreement definition for test-retest reliability studies, which is consistent with our experimental design. The ICC calculation formula was as follows,

$$\begin{align}{\text{ICC}} = \frac{{{\text{M}}{{\text{S}}_{\text{R}}} - {\text{M}}{{\text{S}}_{\text{E}}}}}{{{\text{M}}{{\text{S}}_{\text{R}}} + \left(k - 1\right){\text{M}}{{\text{S}}_{\text{E}}} + \frac{k}{n}\left({\text{M}}{{\text{S}}_{\text{C}}} - {\text{M}}{{\text{S}}_{\text{E}}}\right)}},\end{align} \tag{ 1 }$$

where ${\text{M}}{{\text{S}}_{\text{R}}}$ represents the mean square for rows, ${\text{M}}{{\text{S}}_{\text{C}}}$ represents the mean square for columns, and ${\text{M}}{{\text{S}}_{\text{E}}}$ represents the mean square error. We used ICC statistics with thresholds of ICC > 0.90 indicating excellent stability, ICC > 0.75 indicating good stability, and ICC > 0.50 indicating moderate stability levels (Koo and Li 2016).

We assessed the performance of four types of PSD-based features in cross-day classification. These included the raw PSD, the periodic feature, the aperiodic feature, and a combination of periodic and aperiodic features. Unlike statistical analysis that focuses on regions of interest, we used data from all channels to construct our classification model. The raw PSD, periodic, and aperiodic features had 300 dimensions (60 channels × 5 frequency bands). In contrast, the combined periodic and aperiodic features encompassed 600 dimensions (60 channels × 5 frequency bands × 2). Each one-dimensional component represented the power in a specific frequency band of one electrode. Periodic features came from the periodic components in the PSD, while aperiodic features came from the aperiodic components. We performed binary-class (rest vs. task) and 3-class (rest vs. low MWL vs. high MWL) classifications using the support vector machine method, implemented in the Matlab machine learning toolbox. The RBF kernel was used, and we optimized the parameters through grid optimization (Cortes and Vapnik 1995).

In contrast to the cross-validation method that uses samples from all days, we did cross-day cross-validation using one day's data as the training set and the other day's data as the test set. The performance was evaluated using the classification accuracy, the area under the curve of the receiver operating characteristic (ROC AUC), and the confusion matrix. The ROC AUC results were considered excellent for AUC values between 0.9 and 1, good for AUC values between 0.8 and 0.9, fair for AUC values between 0.7 and 0.8, poor for AUC values between 0.6 and 0.7 and failed for AUC values between 0.5 and 0.6 (Obuchowski 2003). Additionally, to investigate the classification performance using other classifiers, we reported classification results using linear discriminant analysis, naive Bayes, decision tree, and Adaboost.

To quantitively measure both the stability across days and the separability across mental states, we defined a CFSI. The formula was as follows, where the numerator describes the separability of the features and the denominator expresses the day-to-day variability of the features,

$$\begin{align}{\text{CFSI}} = \frac{{{\text{Feature}}{\text{ }}{\text{Separability}}{\text{ }}{\text{Index}}}}{{{\text{Day-to-day}}{\text{ }}{\text{Variability}}{\text{ }}{\text{Index}}}}.\end{align} \tag{ 2 }$$

Here, we assumed that the data set $D = \left\{ {{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}, \ldots ,{{\boldsymbol{x}}_k}} \right\}$ contained $k$ samples with labels, ${{\boldsymbol{x}}_k} \in {\mathbb{R}^l}$. Each sample was an $l$-dimensional feature vector. Dividing the sample set $D$ into $m \times n$ clusters according to the labels ${ }\left\{ {\left. {{C_{ij}}} \right|1 \unicode{x2A7D} {\text{i}} \unicode{x2A7D} m,1 \unicode{x2A7D} j \unicode{x2A7D} n} \right\}$, where $m$ was the number of days and $n$ was the number of classes of features. We defined the feature separability index (FSI) as:

$$\begin{align}{\text{FSI}} &= \frac{1}{m}\mathop \sum \limits_{i = 1}^m \left[ \frac{2}{{n\left( {n - 1} \right)}}\sum\nolimits_{1 \unicode{x2A7D} {j_a} < {j_b} \unicode{x2A7D} n}\right.\nonumber\\ &\left.\quad \times {\left( {\frac{{{\text{dist}}\left( {{{\boldsymbol{u}}_{{j_a}i}},{{\boldsymbol{u}}_{{j_b}i}}} \right)}}{{{\text{avg}}\left( {{C_{{j_a}i}}} \right) + {\text{avg}}\left( {{C_{{j_b}i}}} \right)}}} \right)} \right],\end{align} \tag{ 3 }$$

where ${{\boldsymbol{u}}_{ij}}$ represented the centroid of cluster ${C_{ij}}$, ${\text{dist}}\left( { \cdot , \cdot } \right)$ represented the Euclidean distance between two samples, ${\text{avg}}\left( {{C_{ij}}} \right)$ represented the average Euclidean distance of all samples within cluster ${C_{ij}}$ from the centroid ${{\boldsymbol{u}}_{ij}}$. The centroid is defined as the arithmetic mean position of all the points in the cluster. The FSI indicated the average distance between all the different classes within the same day. A larger FSI value meant a greater separability of the features. In addition, we defined the day-to-day variability index (DDVI) as:

$$\begin{align}{\text{DDVI}}& = \frac{1}{n}\mathop \sum \limits_{j = 1}^n \left[ \frac{2}{{m\left( {m - 1} \right)}}\sum\nolimits_{1 \unicode{x2A7D} {i_a} < {i_b} \unicode{x2A7D} m}\right.\nonumber\\ &\left.\quad \times {\left( {\frac{{{\text{dist}}\left( {{{\boldsymbol{u}}_{{i_a}j}},{{\boldsymbol{u}}_{{i_b}j}}} \right)}}{{{\text{avg}}\left( {{C_{{i_a}j}}} \right) + {\text{avg}}\left( {{C_{{i_b}j}}} \right)}}} \right)} \right],\end{align} \tag{ 4 }$$

it indicated the average distance between days for the same class, and the smaller the DDVI value, the lower the day-to-day variability of the feature. A higher index can be obtained when feature separability is high and day-to-day variability is low.

The averaged PSD of all participants in each condition was shown in supplementary figure 1. The results of two-way repeated measures ANOVAs (2 d × 3 MWL levels) on the raw PSD, periodic components, and aperiodic components were presented in figure 3. For the raw PSD, significant differences in theta band power in the prefrontal midline region and alpha band power in the left and right parietal regions were observed only at the MWL level. Additionally, significant differences in low beta, high beta, and gamma band powers were observed at both time and MWL levels, but without any interaction. Similar to the raw PSD results, the theta and alpha powers for the aperiodic components showed significant differences only at the MWL level, while differences in high beta and gamma powers at the time level were concentrated in the left and right parietal and occipital regions. In contrast, for the periodic components, significant differences in band powers were observed only at the MWL level, with no significant differences at the time level and no interaction. Furthermore, to compare the differences between the three types of PSD-based features on different days, the power spectra of typical channels FZ and OZ are shown in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 was mainly intended to show the spatial distribution of the statistical results for the three PSD-based features at the channel level, and for further analysis, we selected regions with p < 0.01 (F value > 8.29) and typical regions associated with MWL (Gevins and Smith 2000, Antonenko et al 2010) as regions of interest for the analysis (theta: mid-frontal region, alpha: bilateral parieto-occipital region, low beta: central parietal region, high beta: central parietal region, gamma: parieto-occipital region). Tables 1–3 show the results of the two-way repeated-measures ANOVA for the regions of interest, and figure 5 shows the results of the post-hoc tests with the channels included in each region.

Zoom In Zoom Out Reset image size

For raw PSD, the theta and low beta bands were only affected by the MWL factor and not by the time factor, the alpha band was not significantly different under either factor, and the high beta and gamma bands were significantly different under both factors, but there was no interaction effect (table 1). For the periodic component, the power of each band differed significantly only at the level of MWL and did not differ significantly across days (table 2). For the aperiodic components, both theta, alpha, and low beta bands were influenced by MWL factors only, while high beta and gamma were influenced by time factors (table 3).

The ICC of the theta band in the prefrontal theta, high beta, and gamma bands in the parieto-occipital region of the periodic component was larger than that of the raw PSD and the aperiodic component (figure 6). Table 4 illustrates the ICC of the average PSD of the region of interest in five different frequency bands (theta: mid-frontal region, alpha: bilateral parieto-occipital region, low beta: central parietal region, high beta: central parietal region, gamma: parieto-occipital region). The ICCs of the five bands of the region of interest for the periodic component were all greater than 0.7, except for the alpha band in the low (ICC = 0.62, CI₉₅ = 0.25–0.83) and high (ICC = 0.61, CI₉₅ = 0.24-0.83) MWL states.

Zoom In Zoom Out Reset image size

In our data, periodic features achieved the highest classification accuracy in binary (84.2% ± 11.0%) and three-class (58.9% ± 11.5%) classification tasks. This performance was significantly superior to that of raw PSD features, which yielded accuracy rates of (69.9% ± 18.5%, t₁₈ = 4.064, p < 0.001, Cohen's d = 0.932) and (52.7% ± 14.0%, t₁₈ = 2.161, p = 0.044, Cohen's d = 0.496) in binary and three-class classifications, respectively. Additionally, consistent results were obtained in the public dataset. Periodic features exhibited the highest classification accuracy of (93.8% ± 4.6%), significantly outperforming raw PSD features at (73.4% ± 25.4%, t₁₈ = 4.471, p < 0.001, Cohen's d = 0.830). Figure 7 shows the averaged classification accuracy across all subjects.

Zoom In Zoom Out Reset image size

Moreover, to examine the performance of the classification model based on periodic components versus other models, figure 8 shows the results of the ROC AUC analysis. Classifiers trained by periodic feature had the highest ROC AUC in binary and 3-class classification. There was no significant difference in the ROC AUC for the four types of features, and all of them reached the level of 'good' or 'excellent'. To provide more details on the classification results than just the accuracy, as shown in figure 9, we plotted a confusion matrix to visualize the performance of different classification models. The periodic component feature had better performance than other features in identifying the rest and task states. However, in the 3-class classification results, the periodic component, like the other features, performed poorly in distinguishing between high and low MWL. The ROC AUC and confusion matrix results for the public dataset are shown in supplementary figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further investigate the reason for the superior performance of periodic components in cross-day classification, we proposed the CFSI. It evaluated the cross-day separability of features by calculating the distance between the centroids of sample sets with different MWLs on the same day divided by the distance between the centroids of sample sets with the same MWLs on different days. A higher index can be obtained when feature separability is high and day-to-day variability is low. As shown in figure 10, the CFSI of the periodic component feature was higher than that of the raw PSD feature (t₁₈ = 3.664, p = 0.001, Cohen's d = 0.840) and the aperiodic component feature (t₁₈ = 4.276, p < 0.001, Cohen's d = 0.981) in the results of binary classification. Consistent with the binary classification results, the periodic component feature outperformed the other features in the 3-class classification (raw PSD vs. periodic component: t₁₈ = 3.687, p = 0.001, Cohen's d = 0.846; aperiodic component vs. periodic component: t₁₈ = 4.984, p < 0.001, Cohen's d = 1.143). In addition, we also did a statistical analysis of the numerator (FSI) and denominator (DDVI) of the CFSI separately, and the results are shown in supplementary figure 4.

Zoom In Zoom Out Reset image size

Aperiodic activity in the power spectrum varied from day to day, as demonstrated by changes in the slope of the 1/f curve. Figures 3 and 4 show that aperiodic activity differed between two days across all levels of MWL. Specifically, variations from day to day were noticeable in the beta and gamma bands, whereas the theta and alpha bands showed minimal fluctuations. The region most affected by day-to-day variability was the parieto-occipital region. To our knowledge, the reason why aperiodic activity varies from day to day is unclear yet. Conventionally, the aperiodic activity was considered noise in the frequency domain and often required correction. However, recent studies, such as (Gyurkovics et al 2021) and (Ouyang et al 2020), have suggested that it might have functional significance. Gao et al proposed the aperiodic activity as an estimator of excitation-inhibition (E-I) balance using computational modeling (Gao et al 2017). Furthermore, the aperiodic activity has been related to the integration of the underlying synaptic currents, which have a stereotyped double-exponential shape in the time domain that naturally gives rise to the 1/f-like nature of the PSD (Buzsaki et al 2012). Such day-to-day variability may caused by changes in the excitation/inhibition balance (Turrigiano 2011), continuous interaction and competition among the neuronal populations (Kelly et al 2008), and distributed neuromodulation effects (Marder 2012).

In contrast, the periodic component was stable from day to day. The two-way repeated measures ANOVA (see figure 3) indicated that the periodic components in each frequency band did not show significant variations across days. Conversely, the high beta and gamma bands of the raw PSD and aperiodic components were significantly affected by both time and MWL. The ICC analysis results shown in figure 6 further support these findings. The periodic component's consistency was generally higher than that of the aperiodic and raw PSD components. Notably, the FOOOF method extends the PSD-based MWL feature from the low-frequency band to the high-frequency band. The theta and alpha bands of raw PSD did not differ significantly between days (figure 5). This was consistent with previous findings that selected theta and alpha bands of the raw PSD as stable MWL features rather than high-frequency bands (Brouwer et al 2012, Arico et al 2016). Due to the limitations of the narrow-band PSD extraction method, PSDs in high-frequency bands have rarely been used as MWL features in previous studies. We found that the high beta and gamma bands of the periodic component have good day-to-day stability and can distinguish different MWLs because they strip out the aperiodic component (figure 5).

The observed spatial distribution of the raw PSD differs from that of the periodic components (supplementary figure 1). This suggests that the spatial distribution of the raw PSD may be driven, at least in part, by aperiodic activity, as oscillations above the aperiodic component are not separated. Notably, variations in the spatial distribution of theta and alpha bands in the periodic components under different MWL states recapitulate well-established spatial patterns, with theta power concentrated in the frontal midline and alpha power predominantly distributed over posterior areas. However, the spatial distribution of low beta, high beta, and gamma bands of raw PSD differs significantly from that of the periodic components. For the periodic component, the low-beta band was concentrated in the parietal-occipital midline and the high-beta band over posterior areas, with power decreasing as MWL increased. The power of the gamma band in the parieto-occipital region was almost negligible at rest but increased as MWL levels went up. Furthermore, the spatial distribution of the periodic components varies less across days than the raw PSD.

As shown in figure 7, the cross-day classification performance of the periodic component feature was significantly better than the raw PSD feature. Periodic features achieved the highest cross-day classification accuracy in our data and public datasets. This result was consistent with our expectations. Previous results show that MWL can be easily classified using data from the same day (Borghini et al 2014, Aghajani et al 2017). However, the accuracy of the classification model drops when we use test data collected at a different time from the training data. Christensen et al demonstrated that the decline in classification accuracy appears to level off at the hours level and was maintained above chance from one day to up to two weeks later (2012). The variations in PSD-based EEG features from day to day primarily stem from the aperiodic component. Since the raw PSD includes the aperiodic component, training the classifier with raw PSD features leads the classifier to learn not just MWL-related information but also time-dependent patterns. This is a plausible explanation for the lower accuracy observed when using raw PSDs as features for cross-day classification. Moreover, we evaluated the performance of four different feature types on other classifiers (see supplementary tables 1 and 2) and obtained similar results.

The CFSI results (figure 10) show that the periodic component feature outperformed the raw PSD and aperiodic features. Furthermore, the results of the FSI and DDVI can be found in supplementary figure 4. We demonstrated that the day-to-day variability of the periodic feature is significantly lower than that of the raw PSD features (binary-class: t₁₈ = 4.748, p < 0.001, Cohen's d = 1.089; 3-class: t₁₈ = 5.151, p < 0.001, Cohen's d = 1.181), and in contrast, the feature separability is significantly higher than that of the raw PSD (binary-class: t₁₈ = 3.664, p = 0.001, Cohen's d = 0.840; 3-class: t₁₈= 3.687, p = 0.001, Cohen's d = 0.846). Furthermore, our hypothesis was supported by the finding that the day-to-day variability of aperiodic features was significantly higher than that of the raw PSD feature (binary-class: t₁₈ = 4.189, p <0.001, Cohen's d = 0.961; 3-class: t₁₈ = 4.828, p <0.001, Cohen's d = 1.107).

Notably, when we analyzed the day-to-day variability of PSD-based features, we focused on specific regions of interest for our statistical analysis. However, we used data from all channels to construct our classification model. We did this because every channel might have valuable characteristics that can help with classification. Additionally, preprocessing is essential before applying the FOOOF method to extract features. FOOOF extracts features from the PSD curves. It fits Gaussian shapes to the peaks in the spectrum without any prior knowledge of what those peaks represent. Unfortunately, factors like blinking and muscle artifacts can introduce noise into the PSD curve. This noise can lead to unreliable results when extracting periodic components.

The day-to-day variability of PSD-based features has been a significant challenge in the development and practical application of pBCIs. This study suggests that improving the day-to-day stability of features is crucial to avoid BCI calibration. The decoding of MWL states in BCIs is strongly affected by the day-to-day variability of EEG features. Existing classification methods typically require recalibration of the model with new EEG data for each subject, which is time-consuming and inconvenient for users. Compared with the raw PSD features, periodic features for cross-day binary classification improved the accuracy by 14.3% in our dataset and 20.4% in a publicly available dataset, achieving 84.2% and 93.8% accuracy, respectively. This approach can enhance the practicality of PSD-based BCIs without requiring calibration, making them more user-friendly. Notably, the proposed approach of using periodic component features to improve cross-day classification accuracy is not limited to MWL discrimination. Other PSD feature-based classifications, such as emotion recognition, may also benefit from this approach. However, further research is necessary to explore the generalizability of this approach to other brain states.

This study has a few limitations to keep in mind. First, the process of obtaining periodic features involves complex preprocessing, which can make it computationally intensive. This could affect its practical use in real-time classification. Second, EEG data reflects the ever-changing activity of the brain, which can be influenced by various non-cognitive and physiological factors. Despite our best efforts, it is challenging to eliminate the impact of these variables. So, when interpreting our results in the context of MWL assessment literature, it is essential to be cautious.