The HHO is a population-based metaheuristic inspired by the cooperative hunting strategies and surprise pounce behaviors of Harris' hawks in nature. The algorithm mathematically models the dynamic interaction between a swarm of hawks, which act as search agents, and a prey, which represents the current best solution. The central mechanism governing the transition between exploration and exploitation is a non-linear escape energy factor that decreases over the optimization, mirroring the prey's physical exhaustion as it attempts to escape [19]. This design is particularly attractive for high-dimensional FS because it provides a principled, automatic schedule for shifting from broad global search to focused local refinement without requiring hand-tuned decay parameters.

The escape energy E is defined as E = 2 * E0 * (1 - t / T), where E0 is the initial energy state drawn uniformly from the interval [-1, 1] at the beginning of each iteration, t is the current iteration index, and T is the maximum number of iterations. The absolute value of E determines the prey's remaining capacity to escape: when |E| is greater than or equal to 0.5, the prey retains sufficient energy to evade, which triggers soft besiege or soft exploration modes, whereas when |E| is less than 0.5, the prey is exhausted, which triggers hard besiege or hard exploration modes. A second random number r drawn uniformly from [0, 1] further distinguishes besiege behavior from progressive rapid dives, producing four condition-dependent update rules that together encode a rich, multi-mode search strategy.

In the soft besiege regime, which occurs when |E| >= 0.5 and r >= 0.5, the hawks slowly encircle the prey. At the same time, it still has energy to escape, and the position update is given by X(t+1) = X_rabbit(t) - E * |J * X_rabbit(t) - X(t)|, where J = 2 * (1 - r1) is the random jump strength of the rabbit and r1 is uniformly distributed in [0, 1]. In the hard besiege regime, which occurs when |E| < 0.5 and r >= 0.5, the prey is exhausted and the hawks perform a tight encirclement around the current best position, with the update rule X(t+1) = X_rabbit(t) - E * |X_rabbit(t) - X_mean(t)|, where X_mean denotes the mean position of the hawk swarm. Together, these two regimes implement the algorithm's exploitative component.

The two exploration regimes incorporate Lévy flights to model the hawks' progressive rapid dives when the prey escapes. In the soft exploration regime, which occurs when |E| >= 0.5 and r < 0.5, the candidate positions Y and Z are generated as Y = X_rabbit(t) - E * |J * X_rabbit(t) - X(t)| and Z = Y + S * LF(D), where D is the problem dimension, S is a random vector of size 1 x D, and LF denotes the Lévy flight function; the next position is set to Y if its fitness is lower than the current position's, otherwise to Z. In the hard exploration regime, which occurs when |E| < 0.5 and r < 0.5, the mean swarm position replaces the individual hawk position in the computation of Y, which sharpens the search around the swarm consensus. This multi-stage, condition-dependent updating mechanism renders HHO highly adaptive to complex, multimodal search spaces such as those encountered in gene selection.

To provide a comprehensive and historically layered baseline, HHO is compared against ten algorithms spanning three evolutionary eras. The classical evolutionary algorithms include GA, which evolves populations through selection, crossover, and mutation [6], and DE, which generates candidate solutions by adding the weighted difference of two population vectors to a third vector, followed by crossover and greedy selection [8]. The classical swarm intelligence algorithms include PSO, which updates agent velocities based on the personal best and global best positions [7], and ABC, which simulates the foraging behavior of honeybees through employed, onlooker, and scout bee phases [12]. These four algorithms have served as the standard baselines in FS research for nearly three decades and represent the canonical reference points for any new optimizer.

The mid-era hierarchy and behavior-based SA include GWO, which mimics the strict social hierarchy of alpha, beta, delta, and omega wolves and updates positions based on the encircling behavior of the three leading wolves [9]; WOA, which simulates the bubble-net hunting mechanism of humpback whales through spiral and shrinking-circle updates [10]; and SSA, which models the chain-like swarming and jet propulsion of salps in the ocean, with the leader salp guiding the chain and the followers adjusting toward the leader and each other [11]. These three algorithms represent the wave of swarm methods that dominated the metaheuristic literature in the mid-2010s and remain widely used in applied FS studies.

The recent high-novelty SA include MPA, which models the foraging strategies of marine predators through a combination of Lévy and Brownian movements across simulated ecosystem gradients and incorporates a memory-saving mechanism that retains the best historical position of each agent [13]; SMA, which simulates the venous network oscillations of slime mould during food foraging, with adaptive weights that regulate the trade-off between positive and negative feedback during the search [14]; and AO, which mimics the swooping, gliding, and grab-and-carry hunting mechanics of Aquila birds across four distinct foraging modes [15]. These three algorithms represent the current frontier of swarm intelligence research. They are frequently cited as state-of-the-art optimizers, which justifies their inclusion as HHO's most challenging contemporaries. Table 2 consolidates the key characteristics and tuned parameters of all eleven algorithms used in this study.

Table 2 documents the eleven algorithms and their foundational parameters, set to the author-recommended values reported in the original publications. The deliberate use of default parameters, rather than algorithm-specific tuning, reflects a realistic deployment scenario in which a practitioner applies off-the-shelf implementations to a new dataset without extensive hyperparameter optimization. This design choice eliminates the confounding effect of unequal tuning effort, which many comparative studies have identified as a source of bias. The parameter spread also illustrates the diversity of mechanisms under comparison: classical algorithms such as GA and PSO rely on fixed rates, whereas modern algorithms such as GWO and SSA introduce time-varying convergence factors, and the most recent algorithms such as MPA and AO incorporate memory and stochastic perturbation terms. This heterogeneity ensures that the comparative evaluation stress-tests HHO against a representative sample of the design philosophies that currently dominate the metaheuristic literature.

The study uses the seminal Leukemia gene expression dataset introduced by Golub et al. [33], widely regarded as a gold standard in computational biology for evaluating FS algorithms because of its high dimensionality, modest sample size, and clinical significance. The dataset includes 72 bone marrow samples characterized by Affymetrix oligonucleotide microarrays with 7,129 probe sets, partitioned into 47 cases of Acute Lymphoblastic Leukemia (ALL) and 25 cases of Acute Myeloid Leukemia (AML). Distinguishing between ALL and AML is clinically critical because the two subtypes originate from different hematopoietic lineages and require drastically different treatment protocols, with misdiagnosis carrying severe prognostic consequences. The 47-to-25 class ratio also introduces moderate class imbalance, making the dataset well suited for evaluating FS algorithms under realistic clinical conditions, since accuracy alone is a misleading metric on imbalanced data. Table 3 summarizes the dataset's principal characteristics and the preprocessing pipeline applied before optimization.

Table 3 highlights the dual challenge motivating this study: an original feature space of 7,129 genes measured across only 72 samples, yielding a feature-to-sample ratio of approximately 99 to 1 and placing the problem firmly in the small-n, large-p regime where overfitting is the dominant risk. The ANOVA F-test pre-filter reduces the search space from 2^7,129 to 2^1,000 candidate subsets, making wrapper optimization tractable while retaining the biologically relevant signal, because the most discriminative leukemia markers, such as CD33 and ZYX, are known to reside within the top 1,000 genes ranked by F-score. The z-score normalization step is essential because the raw expression values span several orders of magnitude and would otherwise cause features with large numerical ranges to dominate the SVM kernel computations. Stratified 10-fold cross-validation preserves the 47-to-25 class ratio across all folds, ensuring each training and testing partition remains representative of the overall class distribution and that reported performance metrics are not biased by incidental class imbalance in individual folds.

Metaheuristic algorithms are inherently designed for continuous optimization, so they require a transfer function that maps continuous agent positions to binary feature-selection vectors X in {0, 1}^d, where 1 indicates the corresponding gene is selected and 0 indicates it is discarded. The study employs the standard S-shaped sigmoid transfer function, which has been shown to preserve population diversity in high-dimensional spaces more effectively than V-shaped functions because it probabilistically flips bits in proportion to the magnitude of the continuous position rather than deterministically thresholding it [10]. Specifically, the sigmoid value is computed as S(x) = 1 / (1 + exp(-x)) for each dimension, and the binary decision is drawn as x_binary = 1 if a uniform random number in [0, 1] is less than S(x) and 0 otherwise. This stochastic binarization introduces additional exploratory pressure that is particularly valuable in the early iterations of the search, when the algorithm should sample diverse regions of the combinatorial space.

We used an SVM with a Radial Basis Function (RBF) kernel as the wrapper evaluator because of its well-documented robustness in high-dimensional, small-sample scenarios and its ability to capture non-linear class boundaries through implicit feature space expansion [24]. The SVM hyperparameters, namely the regularization constant C and the kernel width gamma, were optimized via an internal grid search within the training folds only to prevent information leakage from tuning on the test fold. The fitness function is designed to maximize classification accuracy while simultaneously penalizing large feature subsets to ensure biological sparsity, and is defined as Fitness(X) = alpha * ErrorRate(X) + (1 - alpha) * (|S| / Total_Features), where ErrorRate(X) = 1 - Accuracy(X), |S| is the cardinality of the selected subset, Total_Features is 1,000, and alpha is a weighting parameter. After rigorous preliminary tuning over alpha in {0.90, 0.95, 0.99}, we set alpha to 0.99, which heavily prioritizes classification accuracy while still providing sufficient gradient information to reduce subset sizes. This bi-objective formulation balances predictive performance and biological interpretability without resorting to a full Pareto-based multi-objective treatment, which would have complicated the comparative analysis.

To ensure absolute fairness in the comparative study, all eleven algorithms share strictly identical computational budgets and environmental conditions. The entire experimental framework was developed in Python, utilizing scikit-learn [34] for the SVM classifier and evaluation metrics and NumPy for the algorithm implementations, and was executed on a standardized computational node equipped with an Intel Xeon E5- 2680 v4 processor and 64 GB of RAM running a 64-bit Linux operating system. We set the population size to 30 agents, the maximum number of iterations to 100, and each algorithm was independently executed for 30 runs with different random seeds to account for stochastic variance. We estimated performance using stratified 10-fold cross-validation, ensuring the ALL-to-AML ratio was preserved across all folds. We set all internal parameters of the eleven algorithms to the standard author-recommended values reported in their foundational publications, as consolidated in Table 2. We permitted no algorithm-specific tuning to simulate a real-world scenario in which a practitioner applies standard out-of-the-box algorithms to a new dataset. Fig. 1 schematizes the end-to-end experimental workflow, from raw microarray input to biological validation. Fig. 1 encodes the study's methodological logic as a sequential pipeline in which each stage constrains and informs the next. The ANOVA pre-filter at stage (b) is deliberately positioned upstream of the optimizer so that all eleven algorithms navigate an identical search space, eliminating any confounding effect of unequal preprocessing. The wrapper loop at stages (d) and (e) is the framework's computational core: the binary optimizer proposes candidate gene subsets, and the SVM-RBF classifier returns a fitness value that drives the next iteration. The stratified 10-fold cross-validation at stage (f) ensures that reported performance reflects generalization to unseen patients rather than in-sample fit. Finally, the bifurcation at stages (g), (h), and (i) reflects the study's dual evaluation philosophy: statistical significance testing establishes that performance differences are not attributable to random chance. In contrast, biological enrichment analysis establishes that the selected genes correspond to genuine disease mechanisms. This integrated design distinguishes the present work from purely computational FS studies that omit either statistical or biological validation.

Table 4 presents the comprehensive classification results averaged over the 30 independent runs. The metrics include accuracy, F1-score, MCC, area under the ROC curve AUC, the average number of selected genes, and the standard deviation of accuracy across runs. We included MCC because it is widely regarded as the most informative single metric for imbalanced binary classification, as it aggregates true and false positives and negatives into a single value in the range [-1, 1] that is insensitive to class prevalence [21]. The F1-score complements MCC by emphasizing the harmonic mean of precision and recall, while AUC summarizes the ranking quality of the classifier across all decision thresholds.

Table 4 demonstrates that HHO decisively outperforms all competitors across every metric, achieving an accuracy of 98.61% with a remarkably low standard deviation of 1.02%, which indicates exceptional stability and robustness against initial random seed variations. Most notably, HHO achieves this peak performance by selecting an average of only 12.4 genes, the smallest subset of any algorithm in the comparison. The classical algorithms PSO and GA require four to five times as many genes to achieve significantly lower accuracy, suggesting their search mechanisms are less effective at isolating the transcriptome's minimal discriminative core. Among the recent algorithms, AO and GWO emerge as the closest competitors. However, they still lag behind HHO by more than 1.4 percentage points in accuracy while requiring more features, which indicates an inability to refine the search space as effectively as HHO's dynamic escape energy mechanism. The MCC values follow the same ranking as accuracy, confirming that HHO's superior performance is not an artifact of class imbalance but reflects genuinely better discrimination between ALL and AML samples.

The failure of classical evolutionary algorithms is particularly pronounced in this binary search space. GA's crossover operator, while effective in continuous spaces, frequently disrupts highly fit building blocks, namely synergistic gene combinations, in a 1,000-dimensional binary vector, which prevents the algorithm from consolidating promising partial solutions. ABC overemphasizes exploration; its random scout bee mechanism is too disruptive in high dimensions and prevents the algorithm from converging on a stable, minimal gene subset. These observations align with prior reports that classical optimizers struggle to balance exploration and exploitation once the binary search space exceeds a few hundred dimensions [4]. Fig. 2 visualizes the convergence behavior of all eleven algorithms, offering further insight into the mechanisms underlying the performance gap documented in Table 4.

Fig. 2 reveals three qualitatively distinct convergence regimes that map onto the algorithm taxonomy established in Section 3. The first regime, premature convergence, is exhibited by PSO and WOA, which converge rapidly by iteration 35 but immediately stagnate; in high-dimensional gene spaces, PSO agents are drawn to the global best too quickly, which depletes population diversity and traps the swarm in local optima formed by redundant gene combinations

while WOA's spiral update lacks the stochastic perturbation needed to escape sharp local optima in discrete binary spaces. The second regime, slow exploitation, is exhibited by ABC and DE, which maintain high diversity but converge painfully slowly and do not reach the accuracy levels HHO achieves by iteration 50, even by iteration 100. The third regime, adaptive balance, is exhibited by HHO, MPA, and AO. Still, only HHO exhibits the characteristic secondary fitness drop around iteration 65 that signals a successful escape from a local plateau via its Lévy-flight-driven soft and hard exploration modes. This secondary descent is the algorithmic signature that lets HHO locate gene combinations the other algorithms miss, and it provides a mechanistic explanation for the performance gap documented in Table 4. Fig. 3 complements this view by

showing the distribution of accuracy across the 30 runs, rather than only the mean trajectory.

Fig. 3 confirms that HHO's superiority is consistent rather than incidental, as shown by the narrow interquartile range of its accuracy distribution and the absence of low-accuracy outliers. HHO's median accuracy lies above the upper quartile of every other algorithm, implying that even the worst HHO run outperforms the best run of GA, ABC, DE, PSO, and WOA. The wide spreads in GA and ABC reflect their sensitivity to random initialization, a known weakness of classical optimizers in high-dimensional binary spaces that manifests clinically as unstable biomarker panels that cannot be reproduced across stochastic restarts. The compact distributions of AO and GWO, by contrast, indicate that these recent algorithms achieve reasonable stability even though their central tendency remains below that of HHO. Together, Figs. 2 and 3 establish that HHO's advantage is simultaneously a higher mean performance, a lower variance, and a qualitatively different convergence dynamic, which collectively justify the statistical analysis reported in the following subsection.

To rigorously establish that HHO's superiority is not due to random chance, we conducted a multi-faceted non-parametric statistical analysis. Because evaluating FS algorithms requires looking at both predictive power and subset size, the Wilcoxon signed-rank test at a significance level of p < 0.05 was applied not only to accuracy but across three distinct metrics: classification accuracy, number of selected features (where lower is better), and the composite fitness value. Table 5 reports the resulting p-values for the pairwise comparisons between HHO and each competitor.

Table 5 demonstrates that HHO is statistically significantly better than all ten competitors at the p < 0.05 level across every metric examined, with no exceptions. The p-values against the classical algorithms (GA, PSO, DE, ABC, WOA) are exceedingly small, on the order of 10^-4 to 10^-7, indicating dominance in both accuracy and feature reduction that is unlikely to be attributable to sampling variability under any reasonable significance threshold. Against recent high-novelty algorithms, the p-values are tighter, ranging from 0.004 to 0.031, confirming that while AO, GWO, and MPA are strong contemporary alternatives, HHO still maintains a statistically significant edge. Notably, the p-values for the selected-features metric remain highly significant even against the strongest competitors, such as GWO and AO, showing that HHO's minimal subset size is a consistent algorithmic behavior rather than a random occurrence. To consolidate this ranking into a single omnibus test, we applied the non-parametric Friedman test to the average accuracies of all eleven algorithms across the 30 independent runs, yielding the mean ranks reported in Table 6.

Table 6 reports that the Friedman test assigned HHO an absolute average rank of 1.00, placing it in the top tier well ahead of the second-place algorithm AO at rank 2.33, with an overall Friedman p-value below 0.001 that rejects the null hypothesis of equivalence among the eleven algorithms. The ranking also reveals a clear stratification: the recent algorithms (HHO, AO, GWO, MPA, SSA, SMA) occupy the top six positions, while the classical algorithms (WOA, PSO, DE, GA, ABC) occupy the bottom five. This stratification supports the broader claim that the field has made measurable progress over the past decade. Still, it also shows that the progress is heterogeneous, since HHO's lead over the next-best recent algorithm is comparable in magnitude to the lead of the recent tier over the classical tier. The combined evidence from the Wilcoxon and Friedman tests therefore establishes HHO's superiority as both pairwise-significant and globally dominant.

Predictive performance must be considered alongside computational cost when assessing the practical applicability of FS algorithms in clinical pipelines. Table 7 reports the average wall-clock runtime per run (in seconds), along with the average number of fitness evaluations and the per-iteration cost, for each of the eleven algorithms. We measured the runtime figures on the standardized computational node described in Section 4.4, including the full 10-fold cross-validation loop.

Table 7 indicates that all eleven algorithms execute approximately 3,000 fitness evaluations per run, corresponding to 30 agents multiplied by 100 iterations, which confirms that the comparison is conducted under a strictly equal computational budget. The per-iteration cost ranges from 1.31 seconds for PSO to 1.63 seconds for MPA, a spread of roughly 24 percent that reflects differences in the complexity of the internal update operators. HHO's per-iteration cost of 1.42 seconds places it in the middle of the range. It is markedly lower than MPA and SMA, which require additional matrix operations associated with their Lévy-flight and venous-network operators. Crucially, HHO delivers its superior predictive performance without incurring a proportional runtime penalty, which strengthens the case for its deployment in time-sensitive clinical pipelines. The classical algorithms PSO, DE, and GA are marginally faster per iteration. Still, their lower per-iteration cost is more than offset by their inferior predictive performance, so the cost-normalized advantage of HHO is even larger than the raw accuracy advantage documented in Table 4.

To validate the clinical utility of the FS process, we extracted the 15 most frequently selected genes by HHO across all 30 runs. We subjected them to GO and KEGG pathway enrichment analysis. KEGG pathway enrichment analysis using standard bioinformatics databases. Remarkably, HHO did not select random, statistically correlated noise; it consistently isolated genes with profound hematological and oncological relevance. Fig. 4 visualizes the selection frequency of the top 15 HHO-selected genes, providing a direct visual assessment of the algorithm's stability in identifying the same biomarkers across stochastic restarts.

Fig. 4 shows that HHO selects its top 15 genes with very high frequency, with the most stable marker CD33 appearing in 29 of 30 runs, indicating that the algorithm converges to a biologically coherent core rather than algorithm-specific artifacts. The contrast with GA and PSO is striking: the two classical algorithms show dark cells only for housekeeping genes such as ACTB and GAPDH, which are highly expressed across all tissues and carry no disease-specific signal, suggesting they overfit to expression variance rather than disease mechanisms. The absence of these housekeeping genes from the HHO selection is itself a meaningful finding, because it demonstrates that HHO's fitness pressure successfully discriminates between biologically informative markers and statistical noise. This qualitative difference in gene composition is ultimately more informative than the quantitative accuracy gap, because it determines whether the resulting biomarker panel can be mechanistically interpreted and translated into a clinical assay. Table 8 details the biological roles of the top five HHO-selected genes and their established associations with leukemia.

MYC 24 Master regulator of cell proliferation Amplified or dysregulated in aggressive leukemias IL8 23 Pro-inflammatory chemokine (CXCL8) Implicated in leukemia tumor microenvironment

Table 8 shows that each of the five most frequently selected HHO genes has a well-documented, mechanistically specific role in leukemia biology, providing strong evidence that the algorithm's selections are biologically meaningful rather than coincidental. CD33 is the definitive surface marker for AML and the target of an FDA-approved immunotherapy, which makes its consistent selection by HHO a strong validation of the algorithm's ability to recover clinically actionable biomarkers. TCF3 is a transcription factor whose chromosomal translocations are primary drivers of pre-B-cell ALL, and its selection by HHO indicates that the algorithm captures the lymphoid axis of the disease as well as the myeloid axis. ZYX, MYC, and IL8 extend the panel into cytoskeletal regulation, proliferative control, and microenvironmental signaling, respectively, spanning the principal biological axes along which ALL and AML diverge. KEGG pathway analysis of the full HHO-selected subset revealed significant enrichment in the Hematopoietic cell lineage pathway (adjusted p = 2.1 x 10^-4) and the Transcriptional misregulation in cancer pathway (adjusted p = 5.4 x 10^-3). In contrast, housekeeping genes and ribosomal proteins dominated the top 15 genes selected by GA and PSO. Fig. 5 visualizes this pathway enrichment as a network.

Fig. 5 confirms that the HHO-selected gene subset is not merely a collection of individually significant markers but a coherent biological network centered on two pathways directly implicated in leukemia pathophysiology. The Hematopoietic cell lineage pathway captures the developmental hierarchy from hematopoietic stem cells to committed myeloid and lymphoid progenitors, the axis along which ALL and AML diverge, and the Transcriptional misregulation in cancer pathway captures the dysregulated transcriptional programs that drive leukemic transformation. That HHO recovers both pathways simultaneously, without any explicit biological prior, indicates that its search dynamics naturally align with the modular structure of disease biology. By contrast, the corresponding network for the GA- and PSO-selected subsets is fragmented and dominated by metabolic and housekeeping nodes, reinforcing the interpretation that classical algorithms overfit to statistical variance rather than mechanistic signal. This biological validation elevates HHO from a purely computational optimizer to a credible tool for hypothesis-driven biomarker discovery.