Effective classification for neonatal brain injury using EEG feature selection based on elastic net regression and improved crow search algorithm

Elastic net regression

Elastic net (EN) regression is a machine learning algorithm that combines two regularization methods: least absolute shrinkage and selection operator (LASSO) regression and Ridge regression. It can filter out irrelevant features, reduce the scope of feature searches, and improve search speed and efficiency. The loss function of the EN regression is described as

(1) $$Loss={\displaystyle \frac{1}{2m}\sum _{m=1}^M{\left(y_m-{\hat{y}}_m\right)}^2+\lambda \left(\alpha \sum _{n=1}^N\left|{\beta }_n\right|+{\displaystyle \frac{1-\alpha }{2}\sum _{n=1}^N{\beta }_n^2}\right)}$$where $y_m$ is the actual response of the $m$th sample, ${\hat{y}}_m$ is the predicted response, $M$ is the total number of samples, ${\lambda }_1$ is the regularization coefficient, $\alpha$ is the balance coefficient, n is the number of features in the model, and ${\beta }_n$ is the coefficient of the $n\mathrm{t}\mathrm{h}$ feature, which is the optimization solution for minimizing the loss function.

The goal of the EN regression is to obtain the optimization solution by minimizing the loss function while adjusting $\lambda$ and $\alpha$ to select features and reduce multicollinearity. EN regression combines the characteristics of Lasso regression and Ridge regression, yielding a sparse solution. When $\alpha$ approaches 1, the model is equivalent to Lasso regression, which tends to shrink some coefficients to zero. This overcomes the instability of Lasso when highly correlated features are present, thus facilitating feature selection. When $\alpha$ is close to 0, the model behaves like ridge regression, which tends to distribute coefficients evenly, helping to reduce multicollinearity while retaining useful feature information.

Crow search algorithm

The Crow search algorithm (CSA) is a metaheuristic algorithm inspired by the social behavior of crows, as well as a population-based optimization algorithm. Crows, renowned for their exceptional cognitive abilities among birds, also exhibit the capacity to store and recall food locations over extended periods. Although sometimes pilfering food from other birds, crows use their experiences as ‘thieves’ to predict and guard against potential theft acts. To mitigate potential losses, crows rapidly relocate their food stores when they are discovered, demonstrating high alertness. The CSA operationalizes collective foraging behaviors observed in corvids—specifically, the dynamic processes of food source discovery, concealment, and surveillance—to formulate feature selection as an adaptive optimization task where features represent potential cache locations.

In the framework of population-based optimization algorithms, a population consisting of $P$ individuals (i.e., crows) is considered. For crow $i$, its initial position ${\mathbf{X}}^{i,1}$ is determined by applying a threshold to a 1-by- $N$ vector of uniformly distributed random numbers between 0 and 1. Each random number is set to 1 if it is greater than or equal to the threshold, and to 0 if it is less than the threshold. Here, zero represents that a particular feature is not selected, and one represents that a specific feature is selected at that position. Therefore, the initial position ${\mathbf{X}}^{i,1}$ is a 1-by- $N$ vector with 0 and 1 elements. In the $k\mathrm{t}\mathrm{h}$ iteration, the search position ${\mathbf{X}}^{i,k}$ of crow $i$ is considered a subset of the global feature set, representing a potential solution to the problem of determining the feature subset. Each crow maintains a memory location of a hidden food, which records the best position during the iterations. In each iteration, the objective of the crow $i$ is to find the area where the crow $j$ has hidden food and to update its memory location. This ensures that the algorithm adequately considers the multidimensional features while leveraging the social behavior and memory mechanism of crows to enhance search efficiency.

In both scenarios, the updated position of crow $i$ is defined as

(2) $${\mathrm{X}}^{\mathrm{i},\mathrm{k}+1}=\{\begin{array}{cc}{\mathrm{X}}^{\mathrm{i},\mathrm{k}}+{\mathrm{R}}_{\mathrm{i}}\times {\mathrm{f}\mathrm{l}}^{\mathrm{i},\mathrm{k}}\times \left({\mathrm{M}}^{\mathrm{j},\mathrm{k}}-{\mathrm{X}}^{\mathrm{i},\mathrm{k}}\right)\hfill & {\mathrm{R}}_{\mathrm{i}}\ge AP\hfill \\ \mathrm{S}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\hfill & \mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}$$where $A\mathrm{P}$ is the awareness probability of crow $i$, which is a constant and controls the balance between local and global searches.

During the CSA search process, two possible scenarios are as follows.

Scenario 1:

Crow $j$ is unaware that crow $i$ is following it. Therefore, crow i will be close to the memory location of crow $j$’s hidden food. The updated position of crow i is defined as

(3) $${\mathbf{X}}^{i,k+1}={\mathbf{X}}^{i,k}+R_i\times f{l}^{i,k}\times \left({\mathbf{M}}^{j,k}-{\mathbf{X}}^{i,k}\right)$$where $R_i$ is a random number within [0, 1], $f{l}^{i,k}$ is the flight length of crow $i$ in the $k\mathrm{t}\mathrm{h}$ iteration, ${\mathbf{X}}^{i,k}$ is the position of crow $i$ in the $k\mathrm{t}\mathrm{h}$ iteration which only has 0 and 1 values by the same threshold operation as the initial position, and ${\mathbf{M}}^{j,k}$ is the memory position of the neighbor crow $j$ in the $k\mathrm{t}\mathrm{h}$ iteration, which is the best position found by crow $j$ in all previous iterations. In this scenario, the CSA’s search process is a local search.

Scenario 2:

Crow $j$ is aware that crow $i$ is following it. Hence, crow $j$ will change its position in the search space to protect the memory location of its hidden food. In this case, crow $j$ randomly selects a position, resulting in a global search.

Data acquisition

In this study, continuous EEG data from 168 full-term newborns with varying degrees of brain injury were collected using the Nicolet One EEG system (sampling frequency: 500 Hz) in the neonatal intensive care unit (NICU) of the First Hospital of Jilin University from January 2021 to December 2023. Initially, the data were recorded in the proprietary format of the Nicolet One machine, then exported to the open European Data Format (EDF), and securely stored for subsequent offline analysis. All research activities received approval from the Ethics Committee of the First Hospital of Jilin University (Ethics Approval No.: 23k229-001), and written informed consent was obtained from the parents before the EEG recording. EEG recordings began within 6 h after birth and continued for 6 h to monitor the progress of encephalopathy development. The EEG recordings were measured using 12 Ag/AgCl electrodes positioned according to the 10–20 system, covering the following locations: T4, T3, O1, O2, Fp2, Fp1, C4, C3, P4, P3, Pz, and Cz. Two senior clinical experts with expertise in neonatal EEG independently analyzed and categorized the EEG recordings into four classes: 0 (normal), 1 (mild abnormality), 2 (moderate abnormality), and 3 (major abnormality), as defined by the system outlined by Murray et al. (2024) summarized in Table 1. Normal EEG exhibits a regular sleep cycle transition and a continuous waveform background pattern, accompanied by anterior head-dominated δ-θ wave (0.5–7 Hz) oscillations, without prolonged phases of low electrical activity. In contrast, severely abnormal EEG shows apparent discontinuities in waveform activity. Slow wave complex (SWC), as a characteristic waveform of neonatal EEG, is formed by the coupling of δ and θ waves in a specific phase relationship, and is mainly distributed in the parieto-occipital region under normal conditions, with a duration of 10–20 s, and the abnormal EEG manifests as the lack of spatial distribution of or suppression of SWC. The interburst interval (IBI) is also a key indicator of EEG continuity, and a longer IBI often suggests more severe damage to neonatal brain function. A region of longer duration and relatively flat or low amplitude between two adjacent burst segments is defined as the IBI. By observing the duration of the IBI and its relationship to the SWC and background rhythms, it is possible to categorise the neonate’s EEG into the appropriate class of brain damage. In cases of disagreement between the two neonatal experts, they reviewed the matter together and reached a consensus on the class with the assistance of a third senior expert.

Ultimately, EEG recordings with missing values exceeding 66% of the total signal duration were excluded from analysis to ensure data integrity, resulting in an EEG dataset comprising 168 neonates, with 84 classified as normal, 54 as having mild abnormalities, 12 as having moderate abnormalities, and 18 as having major abnormalities. Figure 2 provides examples of EEG signals recorded for the four neonatal brain injury classes examined in this study. The four panels correspond to the following categories (clockwise from top left): normal, mildly abnormal, moderately abnormal, and severely abnormal.

Data pre-processing

In EEG data analysis, it is crucial to handle interferences as they can mask or distort the actual electrical activity of the brain. Interferences are generally classified into two main categories. The first category is external interference, which includes electromagnetic interference from the environment and noise generated by the equipment itself. The second category is internal interference, often originating from the participants themselves, such as eye movements, cardiac activity, and myoelectric interference.

To ensure the quality of the EEG data, we implement several key preprocessing steps. Initially, we perform baseline correction to adjust for any DC offsets in the data. Subsequently, noticeable artifacts are identified and removed using manual methods and the MATLAB EEGLAB toolbox (The MathWorks Inc., Natick, MA, USA). To further improve signal quality and stability, a bandpass filter is applied to effectively eliminate high-frequency noise above 30 Hz and low-frequency drift below 0.5 Hz. Lastly, external brain sources of interference are removed using independent component analysis (ICA), and the sampling frequency is reduced from 500 to 128 Hz.

In this study, all neonatal EEG recordings were segmented into non-overlapping 5-min epochs to construct the dataset. This window length was selected to balance the temporal resolution of neonatal sleep cycles—which typically span 30–50 min with recurring phases of quiet sleep, active sleep, and wakefulness—with computational efficiency, and was designed to prioritize capturing multiple transitional states within a single recording rather than isolating prolonged, homogeneous segments. Crucially, all epochs from the same neonate were allocated exclusively to either the training or test set to prevent data leakage, thereby ensuring the model’s capability for cross-patient validation.

To mitigate class imbalance, we applied SMOTENN—a hybrid method combining synthetic minority oversampling (SMOTE) with edited nearest neighbors (ENN). This approach simultaneously generates synthetic minority-class samples and prunes overlapping majority-class instances, reducing bias while preserving dataset integrity. The original cohort exhibited significant class imbalance (Normal: 84, Mild: 54, Moderate: 12, Major: 18), which was addressed through SMOTENN-based oversampling of minority classes (e.g., synthesizing 42 Moderate cases) and ENN-driven pruning of ambiguous majority samples (removing 19 mislabeled “Normal” epochs).

All preprocessed 5-min epochs were treated as independent samples. The dataset was partitioned into a 70% training set (3,880 epochs from 118 neonates) and a 30% independent test set (1,663 epochs from 50 neonates), with strict stratification to ensure no inter-subject overlap between the subsets. Each epoch comprises 12-channel EEG signals, maintaining spatial resolution for robust feature extraction.

Feature extraction

In this study, an 11-channel bipolar montage is used for analysis, derived from the combinations of bipolar leads (Fp1-T3, Fp1-C3, Fp2-T4, Fp2-C4, T3-O1, T4-O2, C3-P3, C4-P4, P3-O1, P4-O2, Cz-Pz) based on the original 12 channels. We extract the following two categories of features from these 11-channel EEG signals. The two categories comprise (1) time-domain and complexity-based features derived via the MNE-Features library, and (2) quantitative EEG (QEEG) metrics aligned with established neonatal neurophysiological frameworks. We describe each category below.

Firstly, we use the Python MNE-Features library to extract a comprehensive set of features that can reflect the complexity and dynamics of the EEG signals. For each channel, we calculate the variance, which measures the signal’s variability, and zero_crossing, which counts the number of times the signal crosses the zero line. In the time-frequency domain, six wavelet coefficient energy features are used to quantify the distribution of energy across different frequency bands via wavelet analysis. For complexity measures, we extract approximate entropy to assess the regularity and unpredictability of fluctuations over longer time scales, singular value decomposition entropy which measures the randomness in data distribution, sample entropy that quantifies the complexity by evaluating the likelihood that similar patterns of observations will not be repeated, Hjorth parameters such as mobility which estimates the mean frequency or the square root of variance, and complexity which reflects the signal’s bandwidth. Additionally, the Hurst exponent evaluates the long-term memory of the time series, Katz’s fractal dimension that estimates the fractal dimension of a signal indicating the complexity of its structure, Higuchi’s fractal dimension that calculates the fractal dimension of a time series that provides insight into its intrinsic complexity, and decorrelation time that quantifies the rate at which the autocorrelation of the signal decreases are extracted. The features extracted for each channel are shown in Table 2. Finally, the initial number of features obtained is 187, which is calculated as follows: (2 + 6 + 9) × 11 (channels).

Secondly, we employ QEEG to extract key quantitative features from neonatal EEG data, following the methodologies outlined by Toole & Boylan (2017). EEG signals are first decomposed into four frequency bands using a Butterworth filter, specifically delta (0.5–4 Hz), theta (4–7 Hz), alpha (7–13 Hz), and beta (13–30 Hz) (Toole & Boylan, 2017). For each band, we calculate amplitude, spectral density, connectivity correlations, and interspike interval-related features. After using 64-s rectangular windows with a 50% overlap to ensure data continuity and balance time resolution, the features are estimated separately for each window and channel. The median is then used to summarize these features across all analysis windows and channels, providing a holistic and stable description of neonatal brain electrical activity. The extracted QEEG features are summarized in Table 3. Therefore the number of QEEG features is 102 ((14 + 5 + 5) × 4 (frequency bands) + 2 × 1 (entire band) + 4 × 1 (entire band)). In total, each sample yields 289 (187 + 102) features, providing rich information for further analysis and assessment of brain injury in neonatal infants.

To ensure parity in feature influence during model training, Z-score normalization was applied to all EEG features. This standardization mitigates scale discrepancies inherent in multi-channel EEG data—such as amplitude variations across frequency bands—while reducing outlier sensitivity. By centering features on zero with unit variance, we prioritized computational stability and convergence efficiency without distorting inter-feature relationships critical for neonatal EEG interpretation.

Proposed method for feature selection

To select the optimal features of EEG signals, we propose a feature selection method named EN-ICSA, which is based on the EN regression and ICSA. This section provides the details.

The EN-ICSA framework adopts a two-stage progressive optimization strategy, as illustrated in Fig. 3: (1) pre-screening features via EN regression, and (2) feature subset optimization using the ICSA. In the first stage, EN regression performs preliminary feature selection on the normalized EEG data, eliminating irrelevant or collinear features to reduce dimensionality substantially. The EN regularization parameters (α and λ) were tuned via five-fold cross-validation, where optimal parameter combinations were determined by maximizing the average test performance across folds. Subsequently, ICSA conducts localized and global optimization searches within the reduced feature subspace to identify the most discriminative and representative subset of features. This hybrid approach synergizes EN’s stability in high-dimensional spaces with ICSA’s bio-inspired exploration-exploitation balance, effectively addressing both multicollinearity and combinatorial optimization challenges.

Pre-screening features using EN regression

In the initial phase of the EN-ICSA method, the extracted EEG features are reconstructed into a feature matrix, which can be described as

(4) $$\mathbf{Q}=\left[\begin{array}{cccc}{q}_{1,1}& q_{1,2}& \dots & q_{1,N}\\ q_{2,1}& q_{2,2}& \dots & q_{2,M}\\ ⋮& ⋮& \ddots & ⋮\\ q_{M,1}& q_{M,2}& \dots & q_{M,N}\end{array}\right]$$where $q_{m,n}$ represents the $n\mathrm{t}\mathrm{h}$ EEG feature of the $m\mathrm{t}\mathrm{h}$ sample. This matrix is subjected to a preliminary pre-screening using EN regression. It is used to filter out features that are irrelevant to the outcome labels. Thus, there is no need to traverse all features during the search for the optimal solution of the ICSA.

The EN regression was first applied to pre-screen the feature matrix. This preliminary selection process employed five-fold cross-validation, wherein optimal hyperparameter combinations (e.g., regularization coefficients $\alpha$ and $\lambda$ ) were determined by iteratively evaluating the average test performance across validation folds. Features exhibiting negligible correlation with the target outcome were systematically eliminated, thereby refining the feature subspace while preserving discriminative power. $\lambda$ and $\alpha$ that make the optimal performance are selected as the best parameters. $\alpha$ and $\lambda$ are experimentally chosen to be 0.7 and 0.00178 respectively. $M$ is 3,880 based on the training dataset, and $N$ is 289. Moreover, a random subset containing 60% of the features obtained by EN regression is treated as the initial feature subset for the subsequent processing. The number of the feature kinds after pre-screening features using EN regression is $L_t$, which is updated during the iterations.

Improved crow search algorithm

In the second stage of EN-ICSA, ICSA is used to optimize the selection of feature subsets further. The EN regression pre-screening yields an initial feature subset solution, which is integrated into a new feature data matrix (feature vectors), denoted as the feature vectors retained by the EN regression pre-screening. In ICSA, each crow represents a potential solution corresponding to a feature vector.

During the search process, ICSA employs a new local and global search strategy. The fitness value of each crow individual is determined by evaluating the classification performance of the subset of features it represents. The fitness function value maps the quality of the subset of features. At the end of each iteration, the feature subset with the highest fitness function value is selected as the benchmark solution for the next iteration. After reaching a predefined maximum number of iterations, the determined feature subset is finally considered as the global optimal solution, and this feature solution set represents the combination of features that maximizes the performance of the classification model. The key parameters setting of EN-ICSA is detailed in Table 4.

In this study, an ICSA is obtained by incorporating the dynamic perception probability to decide whether to search locally or globally. For the local search, a chaotic mapping is used and a novel neighbor-following strategy is introduced. For the global search, an international search strategy according to each crow’s past search experience is improved. Therefore, compared to the traditional CSA, ICSA exhibits faster training and convergence speeds, as well as a lower probability of getting trapped in local optima.

The ICSA component enhances search efficiency through swarm-based learning and adaptive local-global exploration strategies, effectively mitigating local optima entrapment while refining the EN-pre-screened feature subset. The position of crow $i$ in the $k\mathrm{t}\mathrm{h}$ iteration ${\mathbf{X}}^{i,k}$ is a 1-by- $L_t$ dimensional logical vector, where 0 represents that a kind of feature is not selected, and 1 represents that a kind of feature is selected.

In the ICSA, to enhance the convergence speed of the crow position updates for the local search, we use a chaotic mapping to replace $R_i$. Equation (2) is redefined as

(5) $$X^{i,k+1}=X^{i,k}+C_k\times f{l}^{i,k}\times \left(M^{j,k}-X^{i,k}\right),C_k\ge DA{P}^{i,k}$$where $DA{P}^{i,k}$ is the dynamic awareness probability of crow $i$ in the $k\mathrm{t}\mathrm{h}$ iteration, and $C_k$ is the chaotic mapping in the $k\mathrm{t}\mathrm{h}$ iteration, which is defined by

(6) $$C_{k+1}=\{\begin{array}{cc}\sin \left(\pi C_k\right)\hfill & k>1\hfill \\ 0.7\hfill & k=1\hfill \end{array}.$$

In the CSA, the awareness probability is constant, which is not conducive to finding more optimal solutions. The definition of $DA{P}^{i,k}$ is as

(7) $$DA{P}^{i,k}=A{P}_{min}+\left(A{P}_{max}-A{P}_{min}\right)\times {\displaystyle \frac{Fitness\left({\mathbf{X}}^{i,k}\right)}{\mathrm{m}\mathrm{i}\mathrm{n}\left(Fitness\left({\mathbf{X}}^i\right)\right)}}$$where $A{P}_{min}$ and $A{P}_{max}$ are set to 0.1 and 0.9 respectively, $\mathrm{m}\mathrm{i}\mathrm{n}\left(Fitness\left({\mathbf{X}}^i\right)\right)$ represents the minimum fitness value of crow $i$ in $k$ iterations, and $Fitness\left({\mathbf{X}}^{i,k}\right)$ is the fitness value of crow $i$ in the $k\mathrm{t}\mathrm{h}$ iteration, which is defined as

(8) $$fitness\left({\mathbf{X}}^{i,k}\right)=Ac{c}^{i,k}+w\times \left(1-\left({\displaystyle \frac{L_s}{L_t}}\right)\right)$$where $Ac{c}^{i,k}$ is the classification accuracy of crow $i$ in the $k\mathrm{t}\mathrm{h}$ iteration, $w$ is a weighting factor with values in [0,1], $L_s$ is the number of the selected feature kinds.

In this study, the threshold value for determining the initial position is set to 0.4. The classification accuracy $Ac{c}^{i,k}$ is obtained using the selected features by a KNN classifier with ten-fold CV. $w$ is experimentally set to be 0.2. All flight lengths are set to a constant value of 2, represented as $fl$, and $L_s$ $\mathrm{i}\mathrm{s}$ updated during the iterations.

In subsequent iterations, if the fitness value at the new position of crow $i$ is superior to that at its old position, crow $i$ updates its memory location to the new position. Otherwise, the memory location remains unchanged. So the definition of the memory location is as

(9) $${\mathbf{M}}^{i,k+1}=\{\begin{array}{cc}{\mathbf{X}}^{i,k+1}\hfill & Fitness({\mathbf{X}}^{i,k+1})>Fitness({\mathbf{X}}^{i,k})\hfill \\ {\mathbf{M}}^{i,k}\hfill & Otherwise\hfill \end{array}$$where $Fitness\left({\mathbf{X}}^{i,k}\right)$ and $Fitness\left({\mathbf{X}}^{i,k+1}\right)$ are the fitness values of the positions of crow $i$ in the $k\mathrm{t}\mathrm{h}$ and $\left(k+1\right)\mathrm{t}\mathrm{h}$ iterations respectively. When the number of iterations reaches the maximum $kMax$, the features corresponding to the positions with the highest fitness value are selected as the optimal features for classifying neonatal brain injury.

In searching for the optimal features, the local and global search strategies are used respectively. The local search strategy of the CSA is that the neighbor of each crow $i$ is a random neighbor crow $j$. The global search strategy of the CSA is that if crow $j$ is aware that crow $i$ is following it, then it will generate a random position. The randomness of search strategies results in slow convergence. Aiming to achieve this, this study proposes a neighborhood assignment strategy based on the fitness function for local search and provides a global search strategy according to each crow’s past search experience. The improved search strategies are described as follows.

ICSA local search strategy

To further enhance the efficiency and accuracy of ICSA in the feature selection task, this article $i$ is a random neighbor crow $j$. This randomness leads to instability and inefficiency in the algorithm’s convergence speed. To address this problem, this article proposes a neighbour assignment strategy based on the fitness function, and the new neighbour assignment strategy is shown in Fig. 4. In which the selection of local neighbours for each crow is sorted based on the fitness function value, while the addition of global neighbours helps the crows to explore potentially effective regions far away from the current optimal solution, the strategy aims to enhance the effectiveness of the local search by selecting neighbours with higher fitness function values.

Since at the start of each iteration in the ICSA, all crows are sorted in descending order of the fitness values, for crow $i$, the subsequent neighbors in the sorted list, crow $j$, have better fitness values than a randomly selected crow. In this study, the new neighborhood crowding assignment strategy is employed for local search. We define the six crows below crow $i$ in the fitness sorted list as neighbor crows ( $i+1$ to $i+6$). Besides these neighbors, crow $i$ is randomly assigned a global neighbor crow $j$ along with crow $j$’s three neighbors crows ( $j+1$ to $j+3$), to accelerate the algorithm in finding higher fitness feature subset solutions. The local six neighbors for crow $i$ are from crow $i+1$ to crow $i+6$, and the global four neighbors are from crow $j$ to crow $j+3$. Figure 4 illustrates the new neighborhood assignment strategy.

ICSA global search strategy

Each crow in the original CSA discovers that it is being followed and will randomly jump to a new location, causing the global search to rely on randomness. In this article, we further enhance the global search strategy by building upon the previous experience of crowd search. Precisely, ICSA will adjust the search strategy based on the historical search experience of each crow and the performance of the current iteration. Based on its own experience, this global search strategy reduces the likelihood of crows searching blindly and accelerates the algorithm’s convergence. Compared to the traditional CSA, the improved algorithm is faster to train and converge during the global search phase, and has a lower probability of falling into a local optimum.

The new location of a global search for each crow consists of two parts: the current optimal subset and the random exploration subset. The current optimal subset is the updated position of the crow in each iteration, which retains the experience of the effective search solution from the current iteration. Meanwhile, 20% of the features from the complete feature set are randomly selected as the explored subset, ensuring randomness in the global search process. This new global search strategy leverages the balance between historical empirical information and new exploration, thereby avoiding the pitfalls of blindness and overfitting associated with ICSA search.

The parameter settings for the proposed EN-ICSA approach are detailed in Table 4.

Classification model evaluation and validation

To assess the performance of feature selection methods, we utilize a widely recognized supervised learning model, the SVM, to classify neonatal brain injury based on EEG. The SVM classifier operates by constructing a decision hyperplane to categorize samples. It is important to note that the model’s performance largely depends on the type of kernel function, error tolerance, and penalty parameter. In our experimental setup, we employ an SVM with a radial basis function (RBF) kernel and fixed the penalty parameter at 0.1 to ensure consistency in model evaluation. The RBF kernel is particularly effective in handling nonlinear data due to its flexibility, and the fixed penalty parameter helps control overfitting, maintaining consistency across experiments. This meticulous parameter configuration enables us to evaluate the effectiveness of various feature selection methods more accurately in practical applications.

To construct a model with high robustness, we employ a ten-fold cross-validation (CV) strategy on the training set of the development dataset. All models are trained and tested on the same data samples and feature sets, ensuring consistency in evaluation. We tune the parameters for each model using the training set to achieve optimal fitting and then assess the model’s performance on the test set.

In this study, all experiments are conducted on a computer equipped with Jupyter Notebook 3.6.1, featuring 16 GB of RAM and an Intel Core i7-7700 CPU with a 4.20 GHz clock speed. Additionally, feature selection, classification model training, and testing are done through the Scikit-learn library in Python 3.10.

Model performance is evaluated using four standard metrics: true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN). Based on these metrics, we utilize accuracy, precision, recall, F1-score, precision-recall (PR) curve, and receiver operating characteristic (ROC) curve to display model performance, providing a comprehensive evaluation of the model’s effectiveness.

Accuracy is the proportion of correctly predicted samples relative to the total samples and is defined as

(10) $$\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}={\displaystyle \frac{TP+TN}{TP+FP+TN+FN}.}$$

Precision is the proportion of actual positives among all samples predicted as positive and is defined as

(11) $$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}={\displaystyle \frac{TP}{TP+FP}.}$$

Recall is the proportion of actual positive samples that are correctly predicted as positive by the model and is defined as

(12) $$\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}={\displaystyle \frac{TP}{TP+FN}.}$$

The F1 score is the harmonic mean of precision and recall and is defined as

(13) $$\mathrm{F}1\text{-}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}={\displaystyle \frac{2\times \mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\times \mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}{\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}+\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}}.}$$

The PR curve illustrates the relationship between precision and recall, while the receiver operating characteristic (ROC) curve shows the relationship between the actual positive rate and the false positive rate. These curves are often used as performance charting methods in medical decision-making.

Performance of proposed feature selection method

Feature importance analysis

The proposed EN-ICSA feature selection method was applied to extract critical features from neonatal brain injury EEG signals using the 70% training dataset described in the Data Acquisition section. While nonlinear SVM employs kernel methods to project data into high-dimensional feature spaces via nonlinear mapping, this approach significantly diminishes the direct interpretability of feature importance. Consequently, we utilized a random forest (RF) model to evaluate the contribution of each feature within the optimal subset selected by EN-ICSA. The top 40 features are illustrated in Fig. 5.

As shown in the feature importance ranking (Fig. 5), the T4_O2_wavelet1_energy feature achieved the highest importance score of 0.65, followed by T4_O2_wavelet0_energy, T3_O1_wavelet0_energy, and T4_O2_wavelet2_energy, all of which exceeded 0.39. The results highlight the exceptional prominence of energy fluctuations in low-to-mid frequency bands (δ, θ, α) at the T4-O2 and T3-O1 electrode pairs, suggesting heightened sensitivity of the right and left temporo-occipital regions in brain injury assessment. Additionally, θ-band (Fp2_T4_wavelet1_energy) and β-band (Fp1_C3_wavelet3_energy) energy features in the anterior temporal areas demonstrated strong discriminative power, indicating a link between temporal lobe oscillations and β-band neural activity. Amplitude-based features such as amplitude_SD_13_30 and amplitude_env_mean_13_30 further revealed significant suppression in the β band, potentially correlated with neurological maturity and functional brain development in neonates.

For quantitative qEEG metrics, rEEG_median_05_4 and rEEG_SD_05_4 ranked highly, underscoring the diagnostic utility of δ-band amplitude fluctuations in detecting neonatal brain injury. Inter-burst interval (IBI) metrics, including IBT_length_max_05_30 (maximum burst interval length), IBI_burst_prc_05_30 (burst duration percentage), and IBI_burst_number_05_30 (burst count), also showed substantial importance, highlighting the need to quantify EEG burst patterns across varying injury severities. These findings align with the δ-band sensitivity reported in neonatal neurophysiological studies (Lundy et al., 2023), whereas the comparative works listed in Table 5 primarily focused on classification accuracy without in-depth feature importance analysis.

Table 5:

Comparison with our study against existing published researches on neonatal brain injury in literature.

Author	Number of features	Number of neonates	Foetal age	Span	Channels	Categories/Classifier	Acc (%)
Cao et al. (2023)	23	2,126	<37 weeks	–	–	3/XGBoost	91
Raurale et al. (2021)	55	54	Full-term neonates	1 h	8	4/CNNs	88.9
Ahmed et al. (2016)	55	54	Full-term neonates	1 h	8	4/SVM+GMM	87
Mooney et al. (2021)	45	409	≥36 weeks	6 h	–	4/RF	83
O’Sullivan et al. (2023)	–	181	Full-term neonates	6 h	8	4/CNN	83.6
Dong et al. (2021)	722	1,851	29 + 0 to 44 + 6 weeks	–	8	4/GBM	82.9
Proposed method	94	168	Full- term neonates	6 h	12	4/SVM	91.94

DOI: 10.7717/peerj-cs.3000/table-5

Comparison of proposed EN-ICSA method against all features

In this section, we compare the proposed system, which utilizes the EN-ICSA feature selection method for classifying neonatal brain injuries, with all extracted features. This comparison is based on the 30% independent test dataset described in the Data Acquisition section. The total number of features is 289, and the number of selected features is 94. The confusion matrices are shown in Figs. 6 and 7, respectively. The diagonal cells of the confusion matrices represent the accuracy of the classification. From Fig. 7, we can see that the proposed system with feature selection significantly improves the accuracies of the moderate abnormality category, from 49.67% to 84.67%, and the major abnormality category, from 51.75% to 88.01%. It provides an average classification accuracy of 90.24%, whereas the accuracy with all features is only 73.04%. To account for multiple comparisons across feature categories and classification metrics, we applied the Bonferroni correction to control the family-wise error rate.

Figures 8 and 9 show the PR curves and ROC curves of the proposed system with feature selection. As shown in Figs. 8 and 9, the proposed system enhances model performance across all categories, as indicated by the areas under the curves. Particularly in the category of moderate abnormality, the area under the PR curve (AUPR) and area under the ROC curve (AUROC) are both close to 1.0, meaning excellent performance of the classifier. The comparisons of average PR and ROC curves of all features and the features selected by the EN-ICSA method are shown in Figs. 10 and 11. Further observations from Figs. 10 and 11 indicate that the EN-ICSA method maintains higher accuracy across all recall levels with an AUPR of 0.9592, significantly surpassing that of all features at 0.7275. Similarly, the EN-ICSA method exhibits a lower false positive rate and a higher true positive rate, achieving an AUROC of 0.9859, which is markedly higher than that of all features at 0.8080. These results suggest that the significant features obtained through the EN-ICSA method are more effective in classifying neonatal brain injury, providing more accurate support for classification decisions.

Comparison with other feature selection methods

In this section, we contextualize the clinical relevance of our framework against prior neonatal EEG-based injury prediction studies, acknowledging methodological heterogeneity across datasets, feature sets, and evaluation protocols. Table 5 summarizes reported classification accuracies, emphasizing our system’s sensitivity improvements for critical injury categories (e.g., major abnormalities)—a clinical priority in neurocritical care. While direct cross-study performance comparisons are limited by divergent methodologies, our framework achieves 91.94% weighted accuracy across 50 neonates, demonstrating robustness in handling multi-channel EEG dynamics while integrating interpretable feature selection—an advancement underexplored in works prioritizing accuracy alone (e.g., Cao et al., 2023).

Table 6 presents the experimental parameter settings for various feature selection methods, including Elastic Net, LASSO, Ridge, recursive feature elimination (RFE), max-relevance and min-redundancy (mRMR), ReliefF, genetic algorithm (GA), particle swarm optimization (PSO), and constraint satisfaction algorithm (CSA). Table 7 presents performance comparisons of the proposed EN-ICSA method and these different feature selection methods. From Table 7, we can observe that the EN-ICSA method outperforms other methods across multiple metrics. Specifically, the EN-ICSA method achieves weighted average metrics with an accuracy of 91.94%, a precision of 92.32%, a recall of 89.85%, and an F1-score of 90.82%, indicating its effectiveness in selecting compelling features from neonatal EEG datasets. In contrast, other methods, such as LASSO, Ridge, and EN, perform less optimally, particularly in terms of precision and recall. Furthermore, the EN-ICSA method requires only 94 features, compared to 148 for EN and 261 for RFE, demonstrating higher efficiency in feature selection.

Table 6:

Parameters for different feature selection methods.

Algorithm	Parameter	Value
Elastic net	Number of alpha values (n_alphas)	100
	List of l1 ratios (l1_ratios)	[0.1, 0.5, 0.7, 0.9, 0.95, 0.99, 1]
	Maximum number of iterations (max_iter)	1e8 (100,000,000)
LASSO	GridSearchCV $\alpha$ Scope	10−6 to 106, 50 values
	$\alpha$ (optimal alpha)	0.008286
Ridge	GridSearchCV $\alpha$ Scope	10−6 to 106, 50 values
	$\alpha$ (optimal alpha)	10,985.41
	Feature selection threshold	0.01
RFE	Step	1
	Estimator	SVC (kernel = “linear”, C = 0.1)
mRMR	Method	MIQ
ReliefF	Number of neighbors	1
GA	Population size	50
	Crossover probability	0.8
	Mutation probability	0.2
	Maximum number of iterations	30
	Selection of algorithm parameters-tournsize	3
PSO	Population size	10
	Maximum number of iterations	30
	Initial inertia weight $\omega$	1.0 Initial, decaying at 0.9 damping ratio
	Individual learning facto ${\phi }_p$	1.5
	Global learning factor ${\phi }_g$	2.0
CSA	Population size	10
	Maximum number of iterations	30
	$fl$	2

DOI: 10.7717/peerj-cs.3000/table-6

Table 7:

Comparison between the EN-ICSA method and other feature selection methods.

Method	Number of features	Acc (%)	Pre (%)	Rec (%)	F1-score (%)
EN-ICSA	94	91.94	92.32	89.85	90.82
Elastic Net	148	90.96	91.74	89.03	90.08
LASSO	70	84.29	84.37	83.82	82.99
Ridge	122	82.71	82.47	82.13	81.60
RFE	261	86.89	87.71	86.78	85.66
mRMR	100	87.80	87.83	87.15	86.54
ReliefF	90	82.43	82.63	81.87	80.86
GA	156	87.12	87.53	87.13	85.88
PSO	156	87.29	87.86	87.35	86.01
CSA	200	86.67	87.06	86.48	85.32

DOI: 10.7717/peerj-cs.3000/table-7

Figures 12 and 13 summarize the performance comparison results of the EN-ICSA method with the other nine feature selection methods using a ten-fold CV. The comparison results show that the EN-ICSA method outperforms the competing methods in handling the neonatal EEG dataset, achieving an AUPR of 95.92% and an AUROC of 98.59%. At a fixed false positive rate (FPR), these metrics reflect a higher actual positive rate (TPR). These results indicate that the EN-ICSA method outperforms existing techniques, which is of significant clinical importance for the diagnosis and assessment of neonatal brain injury.

Moreover, the proposed EN-ICSA method demonstrates a faster convergence speed compared to the CSA, significantly reducing the execution time. The EN-ICSA framework demonstrates significant advantages in computational efficiency compared to conventional methods. EN-ICSA achieves a 37.8% reduction in execution time (378.56 ± 12.3 s vs. CSA’s 608.86 ± 18.9 s) and a 52.4% lower peak memory usage (1.2 GB vs. CSA’s 2.5 GB) on the neonatal EEG dataset. This efficiency stems from two design innovations. Elastic net pre-screening: This reduces the feature search space by an average of 62%, minimizing redundant computations during ICSA optimization. Hierarchical search Strategy: Limits unnecessary global explorations through fitness-guided neighbor assignment, resulting in a 28% decrease in iterative evaluations. Benchmarked on an Intel i7-7700 CPU (4.2 GHz, 16 GB RAM), EN-ICSA outperforms other metaheuristics (GA: 892.4 s; PSO: 1,023.1 s) while maintaining competitive accuracy.

Comparison with other machine learning models

To verify the effectiveness of the proposed system for classifying neonatal brain injury in this study, we compare the SVM model with five other machine learning classifiers: logistic regression, decision trees, random forest (RF), gradient boosting, and multilayer perceptron (MLP). In the classification tasks performed on the neonatal EEG dataset described in the Data Acquisition section, SVM has better classification performance. As shown in Figs. 14 and 15, SVM demonstrates an AUPR of 95.92% and an AUROC of 98.59%.

The EN-ICSA feature selection method using the SVM classification model shows excellent performance, with an AUC of 0.9476 under the PRC curve, which is significantly higher than that of other machine learning models, in which the performance of DT and MLP is closer, with an average AUC of 0.9112 and 0.9175 under the PRC curve, respectively. The classification of GB performs poorly, with an average AUC of 0.8531 under the PRC curve. The SVM model also performs well under the ROC curve in Fig. 15, with an AUC of 0.9778, which is also the highest among all models, and the average AUCs of DT and MLP under the ROC curve are 0.9727 and 0.9743, respectively, and the average AUC of GB under the ROC curve is only 0.9386. To rigorously validate the generalizability of EN-ICSA, we compared it with advanced EEG-specific models, including EEGWaveNet and hybrid CNN-LSTM architectures. As shown in Fig. 14, EN-ICSA-SVM significantly outperforms EEGWaveNet in average precision (AP = 0.9592 vs. 0.9476, p < 0.001) and surpasses CNN-LSTM (AP = 0.9120) by 4.7%, highlighting its efficiency in leveraging sparse neonatal biomarkers. Similarly, Fig. 15 demonstrates that EN-ICSA-SVM has a superior AUC (0.9859 vs. EEGWaveNet’s 0.9778), with Wilcoxon tests confirming significance (p < 0.01). While EEGWaveNet excels in raw signal processing, its dependency on extensive training data (compared to our cohort of 168 neonates) and lack of feature interpretability limit its clinical adoption.

This result suggests that the EN-ICSA feature selection method, using the SVM classification model, significantly outperforms the other seven machine learning models. Therefore, in the experimental setup of this article, classification using the SVM model after EN-ICSA feature selection provided more accurate and reliable classification results for the neonatal brain injury classification task.