Enhanced swarm optimization for feature selection in electroencephalogram classification: investigating visibility graph and persistent homology-based features

Data source and preprocessing

This study analyzed both non-medical and medical datasets. The primary focus of this study is the non-medical EEG data collected in the previously published work by Liew et al. (2023), which includes recordings from a group of 45 subjects (25 males and 20 females selected based on predefined requirements, such as gender, health, age, and vision). A total of 15 participants were recruited for each of the three age groups (i.e., 18–25, 26–35, and 36+ years). All the subjects had a normal vision or corrected normal vision. Participants identified an image, focused on subsequent visual stimuli, and responded by clicking the mouse when the target image appeared. The experiment consisted of 150 trials, during which the target image was shown 60 times. During data collection, electrodes were placed according to the International 10–20 system, with the right earlobe as the reference and the right wrist as the ground. The sampling rate was 512 Hz, and no filtering was applied to avoid information loss. EEG signals were recorded under three simulated conditions: (a) quiet, (b) low distraction, and (c) high distraction, designed to mimic different real-world environments. An audio clip of regular office noise (≈55 dB) was played to simulate a low-distraction condition, while irregular office sounds (printer, stamping, phone ringing, etc.) were played at a sound level of 70 dB to simulate high-distraction condition. These sound levels were chosen to prevent auditory injury.

After recording, filtering, segmentation, and artifact rejection were performed to eliminate unwanted signals. A finite impulse response bandpass filter with cut-off frequencies of 8–13 Hz and 13–30 Hz was applied to extract alpha and beta bands. The signals were then segmented by trial, and any trials with amplitudes exceeding 100 μV—indicative of body movements or other artifacts—were excluded. Normal EEG amplitudes typically range from 0.5 to 100 μV (Teplan, 2002). While the dataset contains recordings from 21 electrodes, this study focused on five—T5, T6, O1, O2, and OZ—due to their relevance in auditory and visual processing.

Larger datasets generally enhance classification performance, while smaller ones risk overfitting (Althnian et al., 2021). A sliding window technique increases data volume and improves accuracy effectively (Zhao & Liu, 2023). In this study, we applied a non-overlapping sliding window of 250 Hz in one of our experiments to divide the original EEG time series with 512 data points into two segments, discarding the last 12 data points. This technique doubles the sample size for our classification and aids in capturing the transient patterns or anomalies that might be diluted in the entire signal. We presented all the non-medical EEG analysis results in Experiments I to V.

For the medical dataset, a publicly available epilepsy EEG dataset (namely the Bonn dataset) was used to benchmark the proposed methodology. Collected by the Department of Epileptology at the University of Bonn, Germany (Andrzejak et al., 2001), this dataset consists of five subsets labeled as Z, O, N, F, and S. Each subset contains 100 segments of single-channel EEG recordings. Each segment represents a time series signal with a length of 23.6 s (4,097 data points) sampled at 173.61 Hz. Subsets Z and O comprise the EEG recordings from five healthy awake individuals with their eyes open (Z) and closed (O), respectively. Meanwhile, subsets N, F, and S include the EEG signals recorded from five epilepsy patients. EEG signals were measured during seizure-free intervals in subsets N and F, while subset S contains signals measured during seizures. We focused on classifying EEG data from subset S against other subsets, excluding subset O, as in other relevant studies. We standardized the dataset before feature extraction and trimmed each segment to 5.9 s, a segment size used by Zhen, Yue & Liu (2024), to reduce the computational time while ensuring high EEG classification performance. We compared our findings with those of similar classification studies in Comparative Analysis section.

Ethics approval and consent to participate

For the non-medical EEG data collected in the previously published work by Liew et al. (2023), the ethical approval had been obtained from the Medical Research and Ethics Committee from Ministry of Health Malaysia. All data collection performed in Liew et al. (2023) were in accordance with guidelines registered in NMRR-14-1180-21356. Written informed consent was obtained from all participants after a full explanation of the study objectives and procedures.

Phase space reconstruction

Phase space reconstruction (PSR) captures complex brain dynamics, aiding in the extraction of distinguishing EEG features (Kaur et al., 2020). PSR transforms one-dimensional time series into a higher-dimensional data cloud using time delay embedding. Time delay embedding includes two key parameters: the time delay (τ) and the embedding dimension (m). The selection of the optimal values of these parameters has been well-studied, and various strategies have been proposed (Tan et al., 2023). In this study, we determined the optimal τ and m for each EEG time series using the average mutual information and false nearest neighbor methods. We then utilized the embedded EEG data to construct the PDs for feature extraction.

Visibility graph

A visibility graph is a time series analysis technique that converts time series data into a scale-free graph, capturing long-term temporal dependencies and chaotic time series properties (Bhandari et al., 2024). In a visibility graph, each time series point corresponds to a node, with edges forming between nodes that meet the visibility criteria. Specifically, two nodes, $y_i$ at $t_i$ and $y_j$ at $t_j$, are connected if all intermediate nodes $y_k$ satisfy the condition in Eq. (1).

(1) $${\displaystyle \frac{y_i-y_k}{t_k-t_i}>{\displaystyle \frac{y_i-y_j}{t_j-t_i}.}}$$

Different VGs have specific features and applications. In this study, we employed the natural and horizontal VGs with variations (weighted and limited penetrable) to convert our EEG time series into a complex network and extract significant network features for classification.

Natural and weighted visibility graph

A natural visibility graph (NVG) is constructed by mapping consecutive time series points onto the time axis, representing each value as a vertical bar (Azizi & Sulaimany, 2024). An edge is formed between two points whenever the heads of their vertical bars are visible from each other despite the presence of the intermediate bar.

While the NVG is an undirected and unweighted graph, the weighted visibility graph (WVG) extends it by assigning weights to edges based on the degree of visibility between time series points. WVG is often applied to discriminate various EEG signal types, such as detecting epileptic EEG (Supriya et al., 2016). In constructing WVG, different values, such as the Euclidean and horizontal distance, can be assigned as the edge weight. In this study, we used the absolute slope, as defined in Eq. (2), as our edge weight.

(2) $$\left|{\displaystyle \frac{y_j-y_i}{t_j-t_i}}\right|.$$

Figure 2 illustrates the NVG and WVG for a small synthetic EEG signal. For NVG in Fig. 2A, all the edges are assigned the same value regardless of their absolute slope value. In contrast, the larger the slope between two bars, the bigger the weight of the edge between them in the WVG, as shown in Fig. 2B.

Horizontal and limited penetrable visibility graph

Real-world VGs can be extremely large, making them computationally complex and impractical to analyze. To address this, Luque et al. (2009) introduced the horizontal visibility graph (HVG), which reduces connections while preserving the temporal structure. Unlike NVG, HVG connects nodes based solely on horizontal lines of sight. The visibility condition between two nodes $\left(t_i,y_i\right)$ and $\left(t_j,y_j\right)$ with the presence of a middle rod $\left(t_n,y_n\right)$ is defined in Eq. (3).

(3) $$y_i\cdot y_j>y_n,\mathrm{\forall }i<n<j.$$

HVG can also be directed and weighted, but we sticked to undirected and unweighted HVG in this study for simplicity.

Rather than reducing graph edges, the limited penetrable visibility graph (LPVG), which increases connectivity compared to VG, was used in Wang et al. (2016). LPVG enhances the VG algorithm by mapping a time series onto a graph while preserving its temporal characteristics. Moreover, LPVG demonstrates greater noise tolerance than VG. Its connectivity criterion allows a specified number of obstacles to be ignored, as defined in Eq. (4).

(4) $$y_{i+l}<y_j+\left(y_i-y_j\right){\displaystyle \frac{t_j-\left(t_i+t_l\right)}{t_j-t_i},l<j-i.}$$

The parameter l represents the permeability limit, permitting a fixed number of bar penetrations to reduce the impact of noise. Specifically, we applied the limited penetrable horizontal visibility graph (LPHVG) with a penetrable limit of two to extract network features from our EEG time series.

Figure 3 elucidates the HVG and LPHVG of the same synthetic signal as in Fig. 2. It is clear that the number of connections in HVG is fewer (Fig. 3A) compared to its corresponding VG (Fig. 2). Also, every bar in an LPHVG (Fig. 3B) has denser connectivity than its corresponding HVG (Fig. 3A). The pink lines in Fig. 3B are the additional connections between bars when we granted a permeability limit of two to the HVG in Fig. 3A.

Feature extraction

After the preprocessing stage, we extracted complex network features from the constructed VGs. Then, we computed PH for our EEG using different filtrations that work with distinct EEG time series representations, including raw signals, point cloud, and VG. We extracted various topological features, including persistence statistics and other PH representations, from the resulting PDs.

Specifically, PH is a TDA method that captures structural changes in data across multiple spatial resolutions (Phang et al., 2024). PH computation involves representing the data cloud as simplicial complexes by applying filtration to track changes in geometric characteristics of the data and summarizing these changes with a persistence barcode or diagram. Different filtration methods analyze the geometric information of data in distinct ways. Hence, the following subsections elaborate on the critical steps involved in the PH-based approach, including types of filtrations, a brief description of the persistence statistics, and various PH representations.

Types of persistent homology filtrations

In this study, we investigated three types of filtrations: the Vietoris-Rips (VR), sublevel set, and graph filtration. In VR filtration, a distance parameter ε varies, forming an n-simplex when the pairwise distance between all points in the simplex is at most ε (Schweidtmann et al., 2022). At the early stage of the VR filtration, ε is very small, and only 0-simplexes (i.e., single points) exist. As ε increases, edges and filled triangles between data points emerge, creating 1- and 2-simplexes. Once the dataset is encoded as a nested VR complex, one tracks the birth and death of topological features within the simplexes (n-dimension homology classes, H_n) as the filtration progresses. In this study, we analyzed the lifespan of connected components (H₀) and holes (H₁), summarizing them in PD.

In sublevel set filtration, the evolution of the connected components in the raw EEG signal is tracked by focusing on their local critical points. Topological features extracted from our EEG time series through this filtration are limited to a single topological dimension (i.e., H₀) rather than spanning multiple dimensions (for instance, H₀ and H₁ in VR filtration). Instead of a distance parameter used in VR filtration, sublevel set filtration employs a horizontal reference line that moves from the EEG signal’s minimum point to its maximum value as the threshold. As the threshold line ascends, we recorded the time whenever the line reaches the local minima and maxima of the signal, marking the birth and death of a connected component. Specifically, the most recent minimum point is paired with the following maximum point, leaving the oldest minimum point in the sublevel set for the pairing with a later maximum point (Wang, Ombao & Chung, 2018).

The last filtration we employed is the graph filtration presented in the study by Zhen, Yue & Liu (2024). While VR filtration is computed on the two-dimensional point cloud obtained from the phase space reconstruction of the EEG signal, and sublevel set filtration is on the raw EEG signal, graph filtration is carried out on continuous EEG time series transformed into complex networks. Hence, for graph filtration in this study, we first converted our continuous EEG data into discrete VG and computed its adjacency matrix. The adjacency matrix $A=(a_{ij}{)}_{n\times n}$ contains either a zero value or the weight of the edge between two data points (i and j). We then transformed the adjacency matrix into a distance matrix $D=(d_{ij}{)}_{n\times n}$ using Eq. (5).

(5) $$d_{ij}=\{\begin{array}{cc}{A}_{max}-a_{ij},& \mathrm{i}\mathrm{f}i\ne j\\ 0,& \mathrm{i}\mathrm{f}i=j\end{array}.$$

The distance matrix measures the difference between each element and the maximum element $A_{max}$ in A. This transformation ensures that nodes with stronger connections appear closer. The resulting distance matrix becomes the input for generating PDs.

As a preliminary example to showcase the different traits of these three filtrations, we presented the PDs of the T5 signal recorded under quiet condition from Subject One in Fig. 4. The sublevel set filtration will typically yield PD with broader ranges of persistence values, possibly including negative values (see Fig. 4B), as it captures geometric information directly based on signal amplitude rather than pairwise distances between data points as in VR and graph filtrations. When examining the distance between the points and the diagonal line in the PDs, the PDs derived from VR and sublevel set filtrations show points with long lifespans. In contrast, the points in the PD constructed through graph filtration tend to have relatively shorter lifespans overall.

Persistence statistics

After constructing PDs from the EEG data, we extracted multiple persistence statistics features for the classification. We computed the maximum, sum, mean, median, standard deviation, range, interquartile range, and coefficient of covariance of the lifespan (l), birth (b), and death (d) time of H₀ components and H₁ holes in each PD. Notably, we excluded the statistics of H₀ components’ birth as all H₀ have a birth time of zero. The lifetime of the H₀ and H₁ is crucial as longer-persisting events often indicate meaningful structures, whereas shorter-lived ones are more likely to represent noise.

In addition to conventional statistical metrics, we recorded the total count of H₀ and H₁ in each PD. We also counted the number of relevant H₀ and H₁, defined as those that live at least half as long as the most persisting H₀ and H₁. Furthermore, we calculated the Wasserstein (W) and Bottleneck (B) distance between our PD (D) and an empty PD with diagonal only (D₀) using the formulas in Eqs. (6) and (7).

(6) $$W\left(D,D_0\right)={\displaystyle \frac{1}{\sqrt{2}}\sum _{i\in D}\left(l_i\right)}$$

(7) $$B\left(D,D_0\right)={\displaystyle \frac{1}{\sqrt{2}}\underset{i\in D}{sup}\left(l_i\right).}$$

These additional features are essential as they consider the overall structure of the PD rather than individual values, hence reducing the influence of outliers with extreme lifespans.

Persistent homology representations

While PDs encapsulate valuable information about the “shape” of a dataset, their direct application in machine learning and statistical contexts is challenging (Berry et al., 2020). Shape-based features from these diagrams may lack stability, as slight perturbations in the dataset can introduce or eliminate numerous points with short lifetimes (Karan & Kaygun, 2021). To address this, we transformed PDs into one- or two-dimensional functional summaries and vectorized representations, including persistence landscape, entropy, image, silhouette, and Betti curve, to extract additional features. Except for persistence entropy, most PH representations produce finite-dimensional value vectors. Therefore, we reduced them into single-value features through L¹ and L² normalization for the classification. Table 1 provides brief descriptions of each PH representation used in this study.

Visibility graph based features

Based on the VGs generated from the EEG signal, we can directly extract network features for EEG classification, including the average (weighted) degree, average (weighted) shortest path length, modularity, density, transitivity, local and global efficiency, average clustering coefficient, and degree assortativity. Table 2 gives concise definitions and descriptions of each of these features.

Category of features

As this study utilized 71 EEG features obtained from persistent homology and complex network approach, and within the persistent homological based features, numerous PH representations and persistence statistics were employed, we grouped the features into three broad categories, as listed in Table 3. It is restressed that the grouped features from f₃₃ to f₆₂ in persistence statistics category do not include the birth for H₀ as all H₀ have a birth time of zero.

Feature selection

After extracting 71 features from the EEG data, feature selection is necessary to identify the most relevant subset for classification. We applied binary particle swarm optimization (BPSO) algorithm to search for the optimal subset of features from the 71 available that yields the highest classification accuracy. Moreover, we introduced several strategies to address the standard BPSO limitations and enhance its efficiency.

Binary particle swarm optimization

PSO is a population-based algorithm for solving continuous and discrete optimization problems. Inspired by the collective behavior of birds and fish (Qiao et al., 2024), PSO enables individual particles to explore a solution space, adjusting their trajectories based on shared group knowledge to converge toward the global optimum. Each particle i is defined by two attributes: a velocity vector $V_i=\left(v_{i1},v_{i2},v_{i3},\dots ,v_{ij}\right)$ and a position vector $X_i=\left(x_{i1},x_{i2},x_{i3},\dots ,x_{ij}\right)$, where j represents feature j in the solution space. The velocity vector dictates the movement direction, while the position vector corresponds to a candidate solution. The velocity and position updates follow the Eqs. (8) and (9).

(8) $$v_{ij}\left(k+1\right)=\omega \times v_{ij}\left(k\right)+c_1{r}_1\times \left(Pbes{t}_i^k-x_{ij}\left(k\right)\right)+c_2{r}_2\times \left(Gbest-x_{ij}\left(k\right)\right)$$

(9) $$x_{ij}\left(k+1\right)=x_{ij}\left(k\right)+v_{ij}\left(k+1\right).$$

In the equation, $Pbes{t}_i^k$ denotes the best position previously attained by particle i at iteration k, while Gbest represents the best position found by the entire swarm. We used 100 particles and set the maximum iteration to 200 in this study. The acceleration constants $c_1$ and $c_2$ regulate the balance between local exploitation and global exploration. Adjusting the values of $c_1$ and $c_2$ alters the influence of personal and global best on updating the particles’ velocities, facilitating the exploitation of particles around local best, or the exploration of particles in a broader area. The terms $r_1$ and $r_2$ are random values uniformly distributed in $\left[0,1\right]$ to enhance search randomness. The inertia weight parameter $\omega$ controls the search range of the swarm, helping balance global and local search tendencies. In standard PSO, $\omega$ is dynamic, as defined in Eq. (10), to optimize exploration and exploitation.

(10) $$\omega ={\omega }_{max}-k\times \left({\omega }_{max}-{\omega }_{min}\right)/Mk,$$where ${\omega }_{max}$ and ${\omega }_{min}$ represent the maximum and minimum inertia weight, while Mk denotes the maximum iterations. The $c_1$, $c_2$, ${\omega }_{max}$, and ${\omega }_{min}$ are set to 2, 2, 0.9, and 0.4, respectively, in the standard PSO.

PSO was initially designed for continuous optimization, but Kennedy & Eberhart (1997) later introduced BPSO for discrete problems. Unlike standard PSO, where Pbest, Gbest, and position $x_{ij}$ take continuous values, BPSO represents them as binary values (0 or 1) to indicate feature selection in a given solution. Besides, BPSO employs a sigmoid function to transform velocity and update particle positions, as defined in Eq. (11).

(11) $$x_{ij}=\{\begin{array}{cc}1,& \mathrm{i}\mathrm{f}rand<{\displaystyle \frac{1}{1+e^{-v_{ij}}}}\\ 0,& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array},$$where rand is a random number within the range of $\left[0,1\right]$.

PSO is widely used in optimization due to its simplicity, fast convergence, effectiveness, and strong generalization (Xie et al., 2021). However, it has notable limitations, including susceptibility to local optima and inefficient fine-tuning. Since PSO relies solely on the global best solution for coevolution, particles do not exchange information and lack diversity. Moreover, the algorithm often revisits previously explored regions, as particles strictly follow their historical best experiences, leading to premature convergence and early stagnation (Shaqarin & Noack, 2023).

To address these limitations, researchers have developed various enhanced PSO algorithms for different applications. Common improvements involve adjusting parameter distributions, modifying position update strategies, optimizing swarm initialization, and integrating other intelligent algorithms (Qiao et al., 2024; Zhao et al., 2024). However, many PSO variants incorporate complex improvement strategies yet still struggle with local optimum traps and slow convergence (Emambocus et al., 2021). Hence, we proposed relatively simple approaches that enhance the BPSO feature selection performance, especially with feature input from numerous categories.

Improvement strategy 1: swarm division based on feature categories

The first step of BPSO implementation is initiating a population. In this step, we divided our population into multiple sub-swarms based on our extracted feature categories, as in Table 3. We separated our 100 particles into five groups of 20 particles. The first three groups search for the optimal solution in the solution spaces where a particular feature category is excluded (i.e., groups 1, 2, and 3 search within spaces that have no graph, persistence statistics, and PH representations features, respectively). Furthermore, the fourth group explores the area with all features but not H₀ topological features to investigate the discriminative power of 0-dimensional homology. The last group is free to select any features as the optimal solution. Figure 5 visualizes each solution space and the sub-swarm groups that explore the corresponding space as a Venn diagram. This strategy helps prevent particles from sticking to redundant features falsely selected in their historical best search, forcibly pulling them out from the suboptimal solution.

Improvement strategy 2: modified velocities update and clamping

After defining the population and sub-swarm sizes, the particles search for an optimal solution with the given velocity. We refined the velocity update strategy of our BPSO accordingly. Initially, we updated the particles’ velocity based on their personal best position and the global best position discovered by the whole swarm, as in Eq. (8). With the sub-swarm approach, we produced two global best positions, one by each sub-swarm (SSGbest), and one by the whole swarm as usual (SGbest). Thus, we modified the velocity update formula (Eq. (8)) to include the new sub-swarm global best, as shown in Eq. (12). This update method introduces greater diversity to the search for optimal solutions by particles. Moreover, we clamped the velocity of the particles within the range of $\left[-1,000,1,000\right]$ to avoid overly rapid particle convergence.

(12) $$v_{ij}\left(k+1\right)=\omega \times v_{ij}\left(k\right)+c_1{r}_1\times \left(Pbes{t}_i^k-x_{ij}\left(k\right)\right)+c_2{r}_2\times \left(SSGbest-x_{ij}\left(k\right)\right)+c_3{r}_3\times \left(SGbest-x_{ij}\left(k\right)\right).$$

Improvement strategy 3: dynamic inertia weight and acceleration coefficient

As in Eq. (10), we implemented a dynamic inertia weight that gradually decreases from 0.9 to 0.4 by iterations in our BPSO velocity update mechanism. This dynamic weight facilitates better particle search by letting the particles explore a broader area during early iterations and narrowing down the search area afterward. Similarly, we used dynamic acceleration coefficients in this study instead of the fixed coefficients in standard BPSO. We determined the value of our acceleration coefficients in every iteration with Eqs. (13) and (14). We set our $c_{max}$ to 2 and $c_{min}$ to 1. Higher $c_1$ at the beginning of the search enables the exploration of extensive space, while higher $c_2$ and $c_3$ in the subsequent search allow the particles to converge to the optimum.

(13) $$c_1=c_{max}-\left(c_{max}-c_{min}\right)\times {\displaystyle \frac{k}{Mk}}$$

(14) $$c_n=c_{max}+\left(c_{max}-c_{min}\right)\times {\displaystyle \frac{k}{Mk},n=2,3.}$$

Improvement strategy 4: worst solution enhancement with mutation and concatenation

All particles update their position (or solution) after the search in every iteration. Thus, our final improvement strategy is to further update the worst solutions with the mutated sub-swarm global best concatenation. In every three iterations, we identified the two weakest and one best particles in each sub-swarm, which are the particles that achieved the lowest and highest accuracy at that round. Then, we updated the position of the weakest particles with the concatenation of the best position in each sub-swarm. We enhanced a total of 10 particles (10% of the population) each time, and we further mutated the best position concatenation to ensure all 10 particles had different new positions.

We set our BPSO to stop the search if the global best score increases by less than 1% after 10 iterations. Moreover, we implemented a fitness cache that records the fitness score of each solution, allowing the algorithm to retrieve scores for repeated solutions without reclassification. This approach reduces computational costs and enhances the efficiency of our BPSO. We performed five independent BPSO runs for each classification and obtained the feature subset that yielded the highest accuracy. We illustrated the details of our enhanced BPSO in Algorithm 1 and present a flowchart for the algorithm in Fig. 6.

We empirically determined the values of the hyperparameters for our enhanced BPSO, as shown in Table S1 in the Supplemental Information. We considered the appropriate trade-off between computational efficiency and the solution search effectiveness when selecting the hyperparameter values. We also performed a hyperparameter sensitivity analysis to evaluate how different hyperparameter values influence the classification performance and runtime.

Classification

In this study, we conducted binary EEG classification to distinguish between EEG recorded under quiet and distraction conditions by assuming that the effect of low and high sound distraction on EEG patterns has little to no differences. We employed four classifiers, namely logistic regression (LR), support vector machine (SVM), random forest (RF), and k-nearest neighbor (KNN). LR is a popular binary classifier that classifies input with a sigmoid function. SVM with radial basis kernel functions is powerful in separating non-linear data into classes. RF is a decision tree-based classifier that effectively handles large datasets with higher dimensionality. KNN classifies data based on their distance from their neighbors. We implemented these classifiers through the Python sklearn library (Pedregosa et al., 2011) and maintained the parameters’ default values except for LR maximum iterations (1,000) and KNN number of neighbors (3).

Before the classification, we standardized all extracted features by removing the mean and scaling to unit variance using the Python StandardScaler. We then conducted multiple experiments on single-trial EEG classification with various settings, as summarized in Table 4. Most experiments utilized VR filtration and WVG to extract features from alpha EEG signals. However, we alternated the type of VG and PH filtration methods in Experiments II and III respectively to compare their effectiveness. We applied the sliding window technique in Experiment IV to examine its impact on classification performance. Finally, we replaced alpha EEG with beta EEG in Experiment V to study how EEG in different frequency bands responds to surrounding noise.

We employed the accuracy, precision, recall, and F1-score metrics to evaluate our classification results. We also applied the 5-fold cross-validation in our classification to ensure robust results. We implemented the PH and VG techniques using the Python Ripser (Tralie, Saul & Bar-On, 2018), Persim (Saul & Tralie, 2019), Gudhi, NetworkX (Hagberg, Swart & Schult, 2008), and ts2vg libraries.

To assess the significance of the features selected by our enhanced BPSO, we conducted a permutation test using the optimal feature sets that achieved the highest accuracy for each classifier. As a nonparametric method requiring minimal assumptions, permutation testing is particularly appropriate for topological descriptors (Sathyanarayana, Manjunath & Perea, 2025). Under the null hypothesis that the features used have no association with the class labels, we performed a 1,000-permutation test by randomly shuffling the class labels in each iteration. We chose the classification accuracy score as our test statistic and recorded the accuracy scores across the 1,000 classifications with permuted labeling. We then calculated the p-value of our test by dividing the number of significant permutations (i.e., the permutation that yielded a higher accuracy score than the original dataset) by the total number of permutations.

On the other hand, we also performed classification using the standard BPSO to compare its results with those from the enhanced BPSO classification. To assess the statistical significance of their differences, we first applied the Shapiro–Wilk test to examine the normality of the performance differences observed between standard and enhanced BPSO classifications across various classifiers. In the Shapiro–Wilk test, a p-value smaller than 0.05 indicates that the null hypothesis is rejected, where the result is not normally distributed. Then, we employed a parametric paired t-test for scenarios that show normality and a non-parametric Wilcoxon signed-rank test for the opposite. We selected these statistical test methods due to their strong statistical power, robustness, and effectiveness in the presence of outliers (Yavas et al., 2025).

Experiment I: comparison of classification with complete and reduced features set

In our first experiment, we compared the performance of EEG classification using the complete feature set (i.e., the combination of PH representations, persistence statistics, and WVG features) against the reduced set selected through our enhanced BPSO. The classification was conducted separately for the five selected channels using single-trial alpha band EEG signals. From each of the 135 EEG time series (45 subjects × 3 conditions), we extracted 71 features to form a 135 × 71 feature matrix as input for the classification models.

We first explored the list of features selected by our enhanced BPSO for each classification (4 classifiers × 5 channels) in Table 5. For each channel, if none of the features in a particular category (PH representations, persistence statistics, or WVG features) are selected by a specific classifier, then the relevant row of that feature category is filled with “none”. Several features from PH representations are constantly selected, with exceptions for SVM with O1 and LR with O2 signals. Similarly, at least one of the WVG features is included, except for LR and RF with respective T5 and OZ signals. However, the persistence statistics are omitted in half of the classifications (see the row filled with “none” in the column of persistence statistics in Table 5). These findings indicate that PH representation and WVG effectively capture important information from nonlinear and nonstationary EEG signals better than persistence statistics.

To examine the usefulness of persistence statistics in EEG classification, we analyzed the eight cases (i.e., LR and RF with T5 and OZ, SVM with OZ, and KNN with T6, O1, and O2) that include both PH representations and persistence statistics features in the classification by calculating their feature selection proportions (i.e., the ratio of the number of features selected to the total number of features in a particular category). The proportion facilitates the comparison of the number of features chosen across different categories, as we have different numbers of features in each category (i.e., 22 PH representations, 40 persistence statistics, and nine network features, as shown in Table 3).

Figure 7 shows that persistence statistics have a higher selection proportion, indicating that relatively more persistence statistics are chosen than the PH representation features in most of the eight cases (except KNN with O1 and SVM with OZ). This finding suggests that these statistics still offer strong discriminative power in the classification despite being omitted in many cases. On the other hand, some PH representation features might be redundant and carry less essential topological information in particular cases. However, it is noteworthy that the combination of PH representations and WVG features is generally sufficient for the classification models to reach optimal performance in the remaining cases. Thus, the persistence statistics act more like a complementary feature in this experiment.

After applying feature selection, we compared the classification results between complete and reduced feature sets across all channels and classifiers in Table 6. Accompanied by Table 6, Fig. 8 gives the accuracy and F1-score increments achieved in Table 6 across different classifiers. The highest accuracy and F1-score improvements obtained from LR are 16.30% and 14.83% with O1 and T6 signals, respectively. SVM yields the highest accuracy (3.7%) and F1-score (2.06%) gains with OZ signals. Meanwhile, RF achieves the highest accuracy (15.56%) and F1-score (11.16%) increments with OZ signals. KNN with O2 signals gives the highest accuracy (23.71%) and F1-score (17.77%) increments. These suggest that our enhanced BPSO reliably excludes redundant and weakly informative features, resulting in improved classification performance across all channels and classifiers.

Although Fig. 8 indicates that EEG signals from channel T5 do not yield the highest accuracy and F1-score improvement under any classifier, Table 6 reveals that T5 constantly gives the highest accuracy for reduced feature sets across classifiers, except for KNN. A possible explanation for the modest accuracy gain with T5 is that its performance with the complete feature set already outperformed that of the other channels, leaving limited room for further improvement. When focusing on the reduced feature set, T5 produces the highest accuracy for LR (71.85%), SVM (70.37%), and RF (74.81%), with corresponding F1-scores of 81.10%, 81.09%, and 82.66%. For KNN, T5 achieves 74.81% accuracy and an 81.85% F1-score, slightly below the highest scores of 76.30% and 83.62% obtained with O2. These results suggest that features extracted from T5 effectively distinguish EEG signals across auditory conditions, highlighting the sensitivity of the left temporal brain region to surrounding noise.

Besides the overall metric scores reported, we also computed the per-class precision, recall, F1-score, and precision-recall area under the curve (PR-AUC) for alpha T5 and O2 EEG, which yielded the highest classification performance across classifiers. We presented the results in Table S2 and constructed confusion matrices for these classifications in Figs. S1 in the Supplemental Information. Both the low recall and F1-scores for class 0 (i.e., quiet condition), along with the patterns in the confusion matrices, indicate that our classification framework is biased towards classifying distraction EEG, likely due to the class imbalance in our dataset. The results also suggest that distracted EEG with increased alpha and beta activity (Kuang et al., 2024) may generate more distinguishable topological features, leading to higher performance than quiet EEG.

We further computed the permutation importance of the features used in each T5 classification to study the feature contributions in improving the classification. Specifically, we permutated the values of a given feature and assessed the model’s performance over 30 times. Permuting features with high importance will cause a significant performance drop when compared to the model with the original feature. The feature importance plots in Fig. S2 in the Supplemental Information reveal that persistence statistics and landscape consistently rank among the top three most important across all classifiers, except KNN. For KNN classification, the top three features are all network features, including degree assortativity (f₇₁), which were frequently selected across multiple cases in Table 5. These features effectively capture the discriminative topological structure of the EEG and improve the classification results.

In addition, we compared the classifications between the standard and our enhanced BPSO in Table S3 of the Supplemental Information. Our enhanced BPSO produced accuracy and F1-score matching or exceeding the standard BPSO in most cases, except the F1-score from the classification of T5 signal with KNN. This finding suggests that our proposed BPSO outperforms the standard BPSO in classifying EEG recorded under different auditory conditions. We also reported in Table S4 the permutation test results for the four scenarios that achieved the highest accuracy scores with each classifier. The p-values acquired are lower than 0.01 for every classification, indicating that the likelihood of observing such a high accuracy with the selected features purely by chance is less than 1%.

We conducted statistical tests to determine whether the performance differences between the enhanced BPSO, the standard BPSO, and no feature selection were statistically significant across various classifiers. We first checked the normality using the Shapiro–Wilk test (see Table S5). We then performed a paired t-test for cases with a Shapiro–Wilk test result greater than 0.05 and a Wilcoxon signed-rank test for the remaining cases. Table 7 shows the results of the statistical tests. In the direct comparison between the enhanced and standard BPSO, only SVM and RF classifiers yielded p-value less than 0.05 for the accuracy score, suggesting that the observed accuracy differences between these two BPSO are statistically significant for these classifiers. However, the p-values obtained for LR and KNN classifiers exceed 0.05, indicating no significant difference. This lack of significance may be attributed to the high variability in their pairwise accuracy differences (Arévalo-Cordovilla & Peña, 2025). Although the enhanced BPSO achieved statistically significant accuracy improvements over the standard BPSO only for SVM and RF, its comparison with classification without applying any feature selection revealed consistently lower p-values across most classifiers and metrics. This finding suggests that the enhanced BPSO provides greater and more consistent improvements than the standard BPSO, even when statistical significance is not always observed in their direct comparisons.

Experiment II: comparison of classification with different visibility graph techniques

In “Materials and Methods”, we have seen that different VG techniques convert time series into graphs using distinct strategies. Although this can capture varying characteristics of EEG time series, it may impact classification performance. To investigate this, we applied three different VG techniques (namely HVG, NVG, LPHVG) to our EEG data, extracted their corresponding features, and compared the classification results with those from the first experiment using WVG features. We focused on the O2 classification as it yielded relatively high accuracy across classifiers with network features constantly selected for optimal performance (see Tables 5 and 6).

Table 8 summarizes the EEG classification results using different VG techniques with O2 signals. Features extracted from WVG constantly yield the highest classification accuracy, followed by HVG, NVG, and LPHVG, regardless of the classifier used. WVG differs from the other VGs by assigning weight to each edge, capturing both the data points’ connectivity and amplitude dynamics, and allowing the classifier to distinguish between patterns that might appear similar when connectivity is the only consideration. Despite being computationally efficient, HVG, NVG, and LPHVG focus solely on connectivity, leading to a loss of finer quantitative details in the signal.

Apart from comparing the classification performance across different VG techniques, we also examined which network features listed in Table 2 are worth considering using our enhanced BPSO. Table 9 enumerates the features selected in each VG technique across different classifiers for O2 signal classification. The most frequently chosen feature is modularity (10 times), followed by transitivity (eight times), with density and global efficiency selected six times. In contrast, the average weighted degree and clustering coefficient appear least often (three times). This variation suggests that some network features have greater discriminative power for EEG classification. Features such as average weighted shortest path length and modularity that respectively capture the flow of information and global community structure effectively summarize the complex spatiotemporal dynamics inherent in EEG signals. Conversely, metrics like average weighted degree and clustering coefficient focus more on local connectivity and may not differentiate the underlying brain states as robustly.

Experiment III: comparison of classification with different filtrations in persistent homology

Besides the variety in VG construction, different filtrations exist for constructing PH. In our third experiment, we constructed PH for our EEG data through different filtrations, including the VR, sublevel set, and graph filtration, and compared their classification results. Each filtration generates distinct PDs and topological features, potentially influencing classification accuracy. We focused on the classifications with OZ since PH representations and persistence statistics features are constantly selected in these cases (see Table 5).

Table 10 presents the classification results of OZ signals using topological features from different PH filtrations. For LR and SVM, VR filtration features achieve the highest accuracy, while sublevel and graph filtration features yield similarly lower scores. In contrast, RF and KNN perform better with sublevel filtration. Generally, VR and sublevel filtration features constantly produce the highest accuracy, whereas graph filtration performs less well than the others.

As VR filtration constructs simplicial complexes based on distance thresholds in a time-delay embedded data cloud, it effectively captures the signal’s global structure and periodicity. Their associated features may be more linearly separable, and therefore, linear classifiers like LR and SVM can achieve higher accuracy with these features. In contrast, sublevel set filtration extracts topological features directly from the EEG amplitude, which is unavoidably sensitive to nonlinearity and transient local fluctuations in the signal. As a result, non-linear classifiers such as RF and KNN perform better with sublevel filtration features. Meanwhile, graph filtration may lose some discriminative detail during the discretization or matrix transformation, leading to relatively lower classification accuracy for most of the classifiers in Experiment III. Additionally, graph filtration primarily captures local connectivity and lacks the flexibility to reflect dynamic brain network changes over time (Bhattacharya et al., 2025).

Apart from comparing the classification accuracy across different filtrations, we also investigated which PH representations are most relevant for OZ classification using our enhanced BPSO. Table 11 shows the topological features of H₀ and H₁ selected by our enhanced BPSO for each PH representation in OZ classification across different classifiers and filtrations. Notably, sublevel set filtration differs from the other filtrations as it only produces H₀ topological features (see the dashes in all H₁ columns for sublevel rows in Table 11). As for PH representation features, persistence landscape shows strong discriminative power and is chosen most frequently, followed by persistence silhouette, persistence image, and Betti curve features. However, persistence entropy features are not being selected at all. This suggests that persistence entropy features might have minimal influence on classification for our EEG time series. The importance of persistence entropy features is even less than persistence statistics. Similar to Experiment I, nearly half of the persistence statistics features are omitted in this OZ classification. Even so, persistence statistics are the sole features being selected in the graph filtration for LR and RF classifiers. Results from Table 11 also suggest that H₀ topological features carry valuable information comparable to H₁ topological features except they are consistently excluded in VR filtration using SVM classifier.

Experiment IV: comparison of classification with and without segmentation

In our fourth experiment, we repeated the EEG classification by applying a 250 Hz non-overlapping sliding window to investigate its impact on classification performance. We focused on the T5 and O2 classifications, which provide the highest accuracy score for each classifier among all channels. With sliding window segmentation, we doubled our data size from 135 to 270 time series.

Figure 9 compares T5 and O2 classification results with and without sliding windows. For T5, the classification without a sliding window consistently achieves higher accuracy and F1-score across all classifiers. Similarly, for O2, the classification without a sliding window outperforms the classification with a sliding window, except for the SVM F1-score with a minimal difference (80.05% vs 80.49%). Given this negligible variation, the sliding window approach has minor impact on the classification performance with our EEG data. The findings may be attributed to our relatively short EEG time series (with only 512 Hz or 512 data points). Segmentation may disrupt the time series intrinsic pattern and reduce their temporal correlations, especially when the time series is short.

Experiment V: comparison of classification with different frequency bands EEG

In our last experiment, we conducted the EEG classification using beta-band EEG signals. For clarity, we only discuss the T5 and O2 beta EEG classification in this section and present the classification results of the other channels in Table S6 (Supplemental Information). We first showed the results of classifying T5 and O2 signals with complete and reduced feature sets, in Table 12. Similar to the results of alpha-band EEG signals in Experiment I (see Table 6), the enhanced BPSO effectively selects the optimal features and enhances classification performance for beta-band signals across classifiers, except for SVM with T5, which gives the same results regardless of the feature set.

We then compared the classification results across frequency bands to assess their discriminative power and validate our feature selection approach. Figures 10A and 10B depict the accuracy and F1-score improvements with feature selection via our enhanced BPSO for alpha and beta signal classifications. Generally, the increments of the classification scores in alpha signals are higher than or close to those of their corresponding beta signals. This finding suggests that the alpha signal may respond more actively to the surrounding noises than the beta signal. The only exceptions are in the T5 signals with KNN and O2 signals with SVM, where the respective beta signals achieve higher accuracy (17.04% and 2.22% against 14.81% and 1.48%) and F1-score (11.45% and 1.14% against 10.11% and 0.75%) increments than alpha signals. Notably, the largest performance boost from enhanced BPSO is recorded in alpha O2 signals using KNN, increasing accuracy by 23.71% and F1-score by 17.77%.

Although our enhanced BPSO performed unfavorably in distinguishing beta-band EEG signals under different auditory conditions during a light task than its alpha band counterpart, it is still worthwhile to explore the features selected by our BPSO for beta signals from T5 and O2 channels. We explored such feature selection results in multiple binary heatmaps shown in Fig. 11. Similar to Experiment I, most of the persistence statistics features are not selected (see f₂₃–f₆₂ in Figs. 11C, 11D, 11E, 11F and 11G), except in T5 classification with LR and RF, and O2 classification with LR. This may suggest that our enhanced BPSO performs the feature selection consistently across different bands of EEG signals.

The multiple binary heatmaps also reveal that a certain number of PH representations and network features are consistently being selected, as depicted in Figs. 11A, 11B, and 11G. These results indicate that regardless of the EEG frequency and channel, the PH representations and network features always have a higher selection rate by our enhanced BPSO. However, the L² norm of the H₀ average landscape (see f₆ in Fig. 11A) and the average clustering coefficient (see f₇₀ in Fig. 11G) are excluded in all cases. The H₀ landscape feature appears less informative for distinguishing beta signals, while the average clustering coefficient is omitted due to its focus on local connectivity, as observed in our Experiment II.

Comparative analysis with Bonn EEG

To demonstrate the validity of our proposed methodologies, we conducted a comparative analysis using the Bonn dataset. We conducted three classification tasks: one separating healthy from epileptic EEG signals (Z–S), and two distinguishing seizure signals from non-seizure signals recorded in different brain regions (N–S and F–S) of epileptic patients. Table 13 shows the results of each task. We observed that the LR and SVM classifiers achieved 100% accuracy and F1-score in every classification task. RF and KNN classifiers also exhibit their robustness by achieving accuracy and F1-score levels close to or exceeding 99% in all classifications. It suggests that our proposed methodologies allow multiple classifiers to acquire superior performance even with relatively short signal inputs.

Table 14 compares the results obtained from our proposed methodologies with other studies in classifying EEG signals from the Bonn dataset. In classifying seizure-free and seizure signals recorded from epilepsy patients, our proposed methodologies achieve higher accuracy scores compared to other studies. Multiple studies, including our work, achieved 100% accuracy in distinguishing between seizure and non-seizure signals. However, those studies used a longer signal length or applied a sliding window to expand the sample size, unlike our proposed methodologies, which work well even with a small data size. This finding is important since using longer signals or employing additional preprocessing steps to achieve better classification increases computational cost.

Robustness evaluation: independent runs and hyperparameters

While all reported results are based on five independent runs, we further validated the stability and generalizability of our enhanced BPSO by increasing the number of runs fourfold (i.e., to 20). Table S7 summarizes the minimum, maximum, mean, and standard deviation of the results obtained over the 20 runs for our T5 EEG classification across classifiers. The marginal deviation of the mean accuracy between the 5 and 20 runs, as well as the small standard deviation corroborates the stability and robustness of our enhanced BPSO (Zhu et al., 2025). As for the features selected over the 20 runs, Fig. S3 reveals that persistence landscape and silhouette features are constantly selected among other PH representation features, aligning with the results obtained from 5 BPSO runs (see Table 5). These findings demonstrate that increasing the number of independent runs did not significantly alter the derived conclusions, thereby supporting the robustness of the reported findings.

We then conducted a hyperparameter sensitivity analysis focusing on four distinct population sizes and maximum iterations to assess the effect of parameter optimization on our enhanced BPSO. These two parameters were chosen because their variation impacts sub-swarm sizes, dynamic magnitude of the inertia weight, and acceleration coefficients, subsequently affecting both the algorithm’s runtime and search efficiency.

Table S8 in the Supplemental Information presents the average hyperparameter sensitivity results obtained from five independent BPSO runs. We only focused on alpha T5 using LR and KNN, as SVM yielded relatively lower results while RF required a long execution time. While increasing the population size greatly extends the runtime, it only resulted in a small improvement in accuracy, particularly for LR classifier. The effect of the maximum iteration on runtime was rather irregular, though 200 iterations consistently produced the shortest execution time. Overall, our enhanced BPSO demonstrated stable performance across various population sizes and iteration limits with LR. For KNN, a larger population generally improved accuracy but incurred higher computational cost. Therefore, the use of a population size of 100 and a 200 maximum iterations throughout this study was justified, as it gives a good balance between accuracy and runtime.

Moreover, we reported the average execution times for PH and VG computations on 135 signals, as well as the average runtime of five independent BPSO searches with each classifier, using T5 and O2 EEG as examples, in Table S9 in the Supplemental Information. All computations were performed on a machine with an AMD Ryzen 7 6800H CPU (3.20 GHz), 16 GB RAM, and an NVIDIA RTX 3050 GPU. The PH and VG computation took around 18 min, while the BPSO search time ranged from 0.98 to 1.40 min for all classifiers except RF. The BPSO search time with the RF classifier exceeded 23 min due to its exhaustive decision tree strategy.