The molecular docking problem can be formally defined as a continuous optimization problem over a high-dimensional search space determined by the degrees of freedom of the ligand relative to the receptor. A candidate docking solution (pose) is encoded as a vector of decision variables encompassing: three translational parameters (x, y, z) specifying the position of the ligand center of mass within the binding site; four rotational parameters (q1, q2, q3, q4) representing a unit quaternion that defines the orientation of the ligand; and n torsion angles (φ1, φ2, …, φn) corresponding to the rotatable bonds of the ligand. The total dimensionality of the search space is therefore D = 7 + n, where n varies from zero (rigid ligand) to over twenty for highly flexible drug-like molecules [10], [30].

The quality of a given pose is assessed by a scoring function that estimates the binding free energy of the ligand–receptor complex. In the AutoDock 4 framework, the total binding free energy ΔGbind is decomposed as follows:

ΔGbind = Einter + Eintra + Etor, (1)

where Einter is the intermolecular energy between ligand and receptor, Eintra is the intramolecular energy of the ligand, and Etor is the torsional free energy penalty proportional to the number of rotatable bonds. The intermolecular energy is further decomposed into component terms:

Einter = Evdw + Ehbond + Eelec + Edesolv. (2)

Here, Evdw represents van der Waals interactions modeled by a Lennard-Jones 12–6 potential; Ehbond captures hydrogen bonding interactions using a directional 12–10 potential; Eelec accounts for electrostatic interactions via a screened Coulomb potential; and Edesolv estimates the desolvation penalty upon ligand binding using a volume-based pairwise approach [10]. Each term is parameterized using empirical coefficients calibrated against a training set of experimentally determined protein–ligand binding affinities.

AutoDock Vina [11] employs an alternative empirical scoring function based on a combination of Gaussian steric terms (attractive and repulsive), hydrophobic interactions, hydrogen bonding, and a torsional penalty for the number of rotatable bonds. The Vina scoring function emphasizes computational speed by using pre-computed grid maps that discretize the binding site into a regular lattice of energy values, enabling rapid energy evaluation during the search [11].

The energy landscape defined by these scoring functions is highly multimodal, characterized by numerous local minima of varying depth separated by energy barriers. This ruggedness arises from the complex, non-linear interactions between ligand and receptor atoms and is exacerbated by the high dimensionality of the search space for flexible ligands. Consequently, deterministic optimization methods, including gradient descent and Newton’s method, are generally inadequate for molecular docking, as they are prone to entrapment in local minima and provide no guarantee of finding the global optimum [15], [16].

Metaheuristic algorithms are stochastic optimization strategies designed to explore complex search spaces efficiently by balancing global exploration with local exploitation. Their population-based nature, ability to escape local optima, and minimal reliance on gradient information make them well-suited for molecular docking problems [17], [31].

The LGA implemented in AutoDock [10] is perhaps the most extensively used metaheuristic for molecular docking. LGA combines a standard genetic algorithm employing selection, crossover, and mutation operators with a Solis–Wets local search procedure. The “Lamarckian” aspect refers to the inheritance of locally optimized phenotypic improvements, whereby the results of local search are written back into the genotype for subsequent evolutionary operations. This hybridization significantly enhances the algorithm’s ability to refine promising docking poses while maintaining population diversity through genetic operators [10]. PSO, inspired by the collective behavior of bird flocks and fish schools [32], has been adapted for molecular docking in several studies. Namasivayam and Günther [22] developed PSO@AutoDock, a PSO-based docking engine integrated into the AutoDock framework, demonstrating competitive performance with LGA on a diverse set of protein–ligand complexes. More recently, Li et al. [23] proposed RDPSOVina, incorporating random drift PSO into AutoDock Vina to improve sampling efficiency in high-dimensional docking search spaces.

DE, introduced by Storn and Price [24], has been applied to docking by several groups, exploiting the algorithm’s effective perturbation mechanisms for continuous optimization. ACO, inspired by the foraging behavior of ants [25], has also been explored for docking, particularly in the context of fragment-based approaches. Additionally, nature-inspired algorithms including the Grey Wolf Optimizer (GWO) [33], Whale Optimization Algorithm (WOA) [34], Harris Hawks Optimization (HHO) [35], and Gravitational Search Algorithm (GSA) [36] have been applied to various molecular optimization tasks, reflecting the growing diversity of metaheuristic paradigms in computational chemistry.

The METADOCK 2 platform developed by Imbernón et al. [27] represents a significant advance in high-throughput docking, leveraging parallel computing architectures to execute metaheuristic searches across thousands of protein–ligand pairs simultaneously. Chen et al. [5] presented an algorithm selection framework for docking, analyzing the performance of multiple metaheuristics on Angiotensin-Converting Enzyme (ACE) targets and providing guidance for algorithm–problem matching.

Recent years have witnessed increasing interest in deep learning-based approaches to molecular docking, including methods that predict binding poses using graph neural networks and geometric deep learning. However, these approaches remain limited by the availability and quality of training data, and physics-based scoring functions evaluated through metaheuristic search continue to offer superior generalization in many settings, particularly for novel chemotypes and targets [5], [37].

The multi-objective formulation of molecular docking seeks to optimize two or more conflicting objectives simultaneously, yielding a set of Pareto optimal solutions rather than a single optimum. In the bi-objective formulation adopted in this study, the two objectives are:

Minimize f1(x) = Einter(x) (intermolecular energy), (3) Minimize f2(x) = Eintra(x) (intramolecular energy), (4)

where x = [x, y, z, q1, q2, q3, q4, φ1, …, φn] is the decision vector representing a ligand conformation. A solution xA is said to Pareto dominate another solution xB (denoted xA ≤ xB) if and only if f(xA) ≤ f(xB) for all i ∈ {1, 2} and f(xA) < f(xB) for at least one j ∈ {1, 2}. The Pareto optimal set consists of all non-dominated solutions, and its image in objective space forms the Pareto front [37], [38].

Several well-established multi-objective evolutionary algorithms have been applied to molecular docking. NSGA-II, proposed by Deb et al. [39], employs fast non-dominated sorting and crowding distance-based selection to produce diverse, well-converged Pareto front approximations. SMPSO (Speed-constrained Multi-objective PSO) by Nebro et al. [40] adapts PSO for multi-objective optimization with velocity clamping and polynomial mutation. GDE3 by Kukkonen and Lampinen [41] extends DE to multi-objective problems using Pareto dominance and crowding distance. MOEA/D by Zhang and Li [42] decomposes the multi-objective problem into a set of scalar subproblems using weight vectors. SMS-EMOA (S-Metric Selection Evolutionary Multi-Objective Algorithm) by Beume et al. [43] uses the Hypervolume (HV) indicator for selection pressure. In their seminal work, Garcia Godoy et al. [20] compared these five algorithms on a benchmark suite of eleven protein–ligand complexes, evaluating performance using HV, spread, and epsilon indicators alongside RMSD analysis of the predicted binding poses. Their study demonstrated that multi-objective approaches could identify high-quality docking conformations competitive with or superior to the single-objective LGA, while providing a richer set of alternative binding modes for medicinal chemistry analysis. A subsequent study by the same group [21] applied multi-objective docking to EGFR kinase mutants relevant to non-small cell lung cancer, demonstrating the translational value of Pareto-based approaches for understanding drug resistance mechanisms.

The WAO was introduced by Daliri et al. [1] as a population-based metaheuristic inspired by the chemical and physical properties of water molecules. Unlike many nature-inspired algorithms that draw metaphors from biological behavior, WAO is grounded in the physicochemical phenomena governing water at the molecular level, specifically hydrogen bonding, evaporation, and molecular motion.

In WAO, each solution candidate is represented as a water molecule in the search space. The algorithm proceeds through three principal phases in each iteration:

Hydrogen bonding phase: Water molecules form intermolecular hydrogen bonds based on proximity and fitness quality. Molecules with superior fitness attract neighboring molecules, analogous to the way hydrogen bonds create transient clusters in liquid water. This phase facilitates information sharing and drives the population toward promising regions, serving as the primary exploitation mechanism. Mathematically, the position of molecule i is updated as:

xi(t+1) = xi(t) + β × (xj − xi(t)) + α × r × (xbest − xi(t)), (5)

where xj is the bonding partner, xbest is the global best solution, β is the bonding coefficient controlling the strength of attraction, α is a scaling factor, and r is a uniformly distributed random number in [0, 1] [1]. Evaporation phase: Molecules with poor fitness are subject to probabilistic evaporation, whereby they are reinitialized to random positions within the search space. This mechanism is analogous to the phase transition of water from liquid to gas and serves as the primary diversity maintenance (exploration) mechanism, preventing the population from collapsing into a single basin of attraction. The evaporation probability is inversely related to the fitness of the molecule: poorly performing molecules are more likely to evaporate and be replaced [1].

Particle motion phase: All molecules undergo position updates governed by velocity-based dynamics, analogous to the kinetic motion of water molecules. Each molecule’s velocity is influenced by its own inertia, attraction to its personal best position, and attraction to the global best position found by the swarm:

vi(t+1) = w × vi(t) + c1 × r1 × (xlocal − xi(t)) + c2 × r2 × (xglobal − xi(t)), (6) xi(t+1) = xi(t) + vi(t+1), (7)

where w is an inertia weight, c1 and c2 are acceleration coefficients, and r1, r2 are random numbers in [0, 1] [1].

The synergistic interplay among these three phases endows WAO with a robust balance between exploration and exploitation. Daliri et al. [1] validated WAO extensively on CEC benchmark functions and engineering design problems, demonstrating competitive performance against established algorithms including PSO, GA, GWO, WOA, HHO, Sine Cosine Algorithm (SCA), and others. The algorithm showed particular strength on multimodal and high-dimensional problems, which is encouraging in the context of molecular docking, where both characteristics are prevalent.

It is worth noting the relationship between WAO and other water-inspired algorithms. The Intelligent Water Drops (IWD) algorithm proposed by Shah Hosseini [44] models the dynamics of water flowing in rivers and forming paths, and is designed primarily for combinatorial and graph-based problems. In contrast, WAO operates in continuous search spaces and draws its inspiration from molecular-level phenomena (bonding, evaporation, kinetic motion) rather than macroscopic hydrodynamics, making it conceptually and mechanistically distinct from IWD.

The transition from a single-objective to a multi-objective optimization framework necessitates fundamental modifications to the WAO algorithm’s selection, ranking, and archive mechanisms. In the single-objective WAO, the scalar fitness value directly determines bonding partners, evaporation candidates, and the global best solution. In the proposed MO-WAO, these decisions must be guided by Pareto dominance relationships and diversity preservation criteria, as no single solution can be declared universally optimal when multiple conflicting objectives are present.

The proposed MO-WAO incorporates the following key multi-objective components:

Non-dominated sorting: at each iteration, the population is partitioned into successive non-dominated fronts F1, F2, …, Fk using a fast non-dominated sorting procedure analogous to that of NSGA-II [39]. Front F1 contains all non-dominated solutions in the population; front F2 contains all solutions dominated only by members of F1; and so forth. The front index of a molecule determines its rank for selection and evaporation purposes.

Crowding distance computation: within each non-dominated front, a crowding distance metric is computed for each solution, measuring the average side length of the cuboid defined by its nearest neighbors in objective space [39]. Solutions with larger crowding distances are located in less crowded regions of the front and are therefore preferred for selection, promoting diversity along the Pareto front approximation.

External archive: an external archive A of maximum size Amax is maintained to store the best non-dominated solutions encountered throughout the search. At each iteration, newly generated non-dominated solutions are added to the archive, and dominated solutions are removed. When the archive exceeds Amax, solutions with the smallest crowding distances are pruned to maintain a well-distributed representation of the Pareto front. Modified hydrogen bonding: in the multi-objective context, bonding partners are selected from the external archive rather than based on scalar fitness. A roulette wheel selection mechanism is employed, with selection probability proportional to crowding distance: archive members in less crowded regions of the front are more likely to serve as bonding attractors. This selection strategy encourages the population to explore underrepresented portions of the Pareto front. The position update for the hydrogen bonding phase becomes:

xi(t+1) = xi(t) + β × (xarchive − xi(t)) + α × r × (xleader − xi(t)), (8)

where xarchive is the selected archive member and xleader is another randomly chosen archive member.

Modified evaporation: the evaporation probability for molecule i is determined by its domination count (the number of solutions that dominate it) rather than by scalar fitness. Molecules that are dominated by many others in the population have a higher probability of evaporation and reinitialization, directing computational resources away from poor regions of the objective space:

pevap(i) = ndom(i) / N, (9)

where ndom(i) is the number of solutions dominating molecule i, and N is the population size.

Modified particle motion: the velocity update in the motion phase replaces the single global best with an archive leader selected via roulette wheel on crowding distance:

vi(t+1) = w × vi(t) + c1 × r1 × (xlocal − xi(t)) + c2 × r2 × (xarchive − xi(t)), (10)

where xlocal is the personal best position based on Pareto dominance (updated when the new position dominates the stored personal best), and xarchive is a leader selected from the external archive [40], [41].

Each water molecule in the MO-WAO population represents a candidate ligand conformation, encoded as a real-valued vector:

x = [x, y, z, q1, q2, q3, q4, φ1, φ2, …, φn]. (11)

The translational variables (x, y, z) specify the position of the ligand center of mass within the binding site and are bounded by the dimensions of the docking grid box, which is centered on the active site and typically extends 15–30 Å in each dimension. The rotational quaternion (q1, q2, q3, q4) defines the orientation of the ligand and is normalized to unit length (||q|| = 1) after each position update to maintain a valid rotation representation. The torsion angles (φ1, …, φn) are bounded within [−π, π] and define the internal conformation of the ligand through rotation about its n rotatable bonds. The total solution dimensionality is D = 7 + n, which ranges from 7 for rigid ligands to over 25 for highly flexible molecules [10], [20].

The two objective functions optimized by MO-WAO are defined as follows:

Objective 1 (intermolecular energy (f1)). This objective quantifies the strength of the interaction between the ligand and the receptor. It is computed as the sum of pairwise interaction energies between all ligand atoms and all receptor atoms:

f1(x) = Einter = ∑i∈L ∑j∈R [Evdw(rij) + Ehbond(rij) + Eelec(rij) + Edesolv(rij)], (12)

where L and R denote the sets of ligand and receptor atoms, respectively, and rij is the interatomic distance. More negative values indicate stronger, more favorable binding interactions.

Objective 2 (intramolecular energy (f2)). This objective captures the internal conformational strain of the ligand arising from intramolecular contacts and torsional strain:

f2(x) = Eintra = ∑i,j∈L, i<j [Evdw(rij) + Ehbond(rij) + Eelec(rij)], (13)

where the summation is over all non-bonded pairs of ligand atoms separated by more than three covalent bonds; lower intramolecular energy indicates a less strained, more relaxed ligand conformation.

Both objectives are to be minimized. Energy evaluation is performed using the AutoDock 4.2 force field, with pre-computed grid maps enabling rapid interpolation of receptor-dependent energy contributions at arbitrary ligand atom positions [10], [30].

The complete MO-WAO procedure for molecular docking is summarized in the following pseudocode: Algorithm 1. MO-WAO for Molecular Docking.

Input: Protein structure, Ligand structure, Grid parameters, Population size N, Max iterations T, Archive size A_max. Output: Pareto front approximation (Archive A) 1. Initialize population P of N water molecules with random docking conformations within search bounds 2. Evaluate f₁(x) and f₂(x) for each molecule in P 3. Initialize external archive A with non-dominated solutions from P 4. FOR t = 1 TO T DO: 4a. Perform fast non-dominated sorting on P → fronts F₁, F₂, ..., F_k 4b. Compute crowding distance for each molecule within each front ── HYDROGEN BONDING PHASE ── 4c. FOR each molecule i in P: • Select bonding partner x_archive from A via roulette wheel (probability ∝ crowding distance) • Select leader x_leader randomly from A • Update position: x_i(t+1) = x_i(t) + β·(x_archive − x_i(t)) + α·rand·(x_leader − x_i(t)) ── EVAPORATION PHASE ── 4d. FOR each molecule i in P: • Compute domination count n_dom(i) • Compute p_evap = n_dom(i) / N • IF rand < p_evap THEN: Reinitialize x_i randomly within search bounds ── PARTICLE MOTION PHASE ── 4e. FOR each molecule i in P: • Select archive leader x_archive via roulette wheel on crowding distance • Update velocity: v_i(t+1) = w·v_i(t) + c₁·rand·(x_local − x_i(t)) + c₂·rand·(x_archive − x_i(t)) • Update position: x_i(t+1) = x_i(t) + v_i(t+1) ── BOUNDARY HANDLING ── 4f. FOR each molecule i in P: • Clamp translation variables to grid box • Normalize quaternion to unit length • Wrap torsion angles to [−π, π] ── EVALUATION AND ARCHIVE UPDATE ── 4g. Evaluate f₁(x) and f₂(x) for all updated molecules 4h. Update personal bests using Pareto dominance 4i. Add non-dominated solutions from P to A 4j. Remove dominated solutions from A 4k. IF |A| > A_max THEN:

Remove members with smallest crowding distance until |A| = A_max 5. RETURN Archive A as Pareto front approximation

The computational complexity of MO-WAO per iteration is analyzed as follows. The fast non-dominated sorting procedure requires O(M × N2) operations, where M = 2 is the number of objectives, and N is the population size [39]. Crowding distance computation requires O(M × N log N) for sorting along each objective. The hydrogen bonding, evaporation, and motion phases each require O(N) position updates with O(D) operations per update, where D is the solution dimensionality. Energy evaluation dominates the per-molecule cost at O(E) per evaluation, where E depends on the number of ligand and receptor atoms and the grid interpolation complexity. The overall time complexity for T iterations is:

O(T × (M × N2 + N × D + N × E). (14)

This complexity is comparable to that of NSGA-II, which also has O(M × N2) non-dominated sorting cost per generation. The space complexity is O((N + Amax) × D) for storing the population and archive.

The parameter configuration of MO-WAO is summarized in Table 1. Parameters were selected based on preliminary sensitivity analysis and guidelines from the original WAO study [1].

The population size of 100 and the maximum function evaluation budget of 25,000 are consistent with the default settings of AutoDock and the configurations used in the comparative studies of Garcia Godoy et al. [20], [21], ensuring fair comparisons across algorithms. The inertia weight is linearly decreased from 0.9 to 0.4 over the course of the run to shift the balance from exploration to exploitation as the search progresses, a strategy well-established in PSO literature [32]. All competitor algorithms are configured with identical population size, evaluation budget, and archive size to enable equitable comparison.

The performance of MO-WAO is evaluated on a set of eleven protein–ligand complexes drawn from the Protein Data Bank (PDB) [45], consistent with the benchmark suite employed by Garcia Godoy et al. [20]. These complexes span diverse protein families and ligand flexibility levels, providing a representative cross-section of molecular docking challenges encountered in early-stage drug design. Table 2 summarizes the benchmark instances.

The benchmark instances range in complexity from the relatively rigid biotin–streptavidin complex (1STP, 5 rotatable bonds, D = 12) to the highly flexible SB203238–HIV-1 protease complex (1HBV, 19 rotatable bonds, D = 26). Several HIV-1 protease complexes are included, reflecting the therapeutic importance of this target in antiretroviral drug development and enabling systematic assessment of algorithm performance across different inhibitor scaffolds bound to the same protein. The acetylcholinesterase complex (1D4K) and the CDK2 complex (1G5S) provide additional target diversity relevant to Alzheimer’s disease and cancer therapeutics, respectively.

The proposed MO-WAO is compared against five well-established multi-objective metaheuristics and one single-objective reference method:

I. MO-WAO (proposed): MO-WAO as described in Section 3.

II. NSGA-II [39]: NSGA-II employs fast non-dominated sorting, crowding distance-based selection, and elitist preservation through a combined parent–offspring population. NSGA-II is widely regarded as the gold-standard evolutionary multi-objective optimizer and serves as the primary baseline for comparison.

III. SMPSO [40]: SMPSO applies velocity clamping and polynomial mutation to a PSO framework with Pareto dominance-based archiving. SMPSO has demonstrated strong performance on continuous multi-objective problems.

IV. GDE3 [41]: GDE3 extends DE to multi-objective optimization using Pareto dominance for survivor selection and crowding distance for tie-breaking. GDE3 benefits from the effective perturbation operators of DE.

V. MOEA/D [42]: MOEA/D decomposes the multi-objective problem into a set of scalar subproblems defined by uniformly distributed weight vectors and solves them collaboratively through neighborhood-based reproduction. The Tchebycheff decomposition approach is used.

VI. SMS-EMOA [43]: the S-Metric Selection Evolutionary Multi-Objective Algorithm uses the HV indicator as a selection criterion, replacing the solution whose removal results in the smallest HV loss. SMS-EMOA provides strong convergence pressure but is computationally more expensive per generation.

VII. LGA [10]: the LGA from AutoDock 4.2 serves as a single-objective reference, optimizing the total binding free energy (ΔGbind = Einter + Eintra + Etor). LGA provides a baseline for assessing whether multi-objective approaches offer advantages over the established single-objective standard.

All multi-objective algorithms are configured with a population size of 100, a maximum of 25,000 function evaluations, an archive size of 100, and 30 independent runs per benchmark instance. These settings ensure fair, controlled comparisons and are consistent with the experimental protocols of Garcia Godoy et al. [20].

The quality of Pareto front approximations is assessed using four widely used performance indicators [46], [47]:

HV: the HV indicator measures the volume of objective space that is dominated by the Pareto front approximation and bounded by a reference point r. Higher HV values indicate superior approximations that are both well-converged toward the true Pareto front and well-distributed across the objective space. The reference point is set to the worst objective values observed across all algorithms for each benchmark instance. HV is the only unary indicator that is both Pareto-compliant and sensitive to both convergence and diversity [47].

Inverted Generational Distance (IGD): IGD computes the average Euclidean distance from each point in the reference Pareto front to the nearest point in the approximation. Lower IGD values indicate closer proximity to the reference front. IGD captures both convergence and diversity in a single scalar metric and requires a reference front, which is constructed from the union of non-dominated solutions found by all algorithms across all runs [46].

Spread (Δ): the spread indicator evaluates the extent and uniformity of distribution of solutions along the Pareto front approximation. It considers both the Euclidean distances between consecutive solutions and the distances to the extreme points. Lower Δ values indicate more uniform distribution and better coverage of the objective space [39].

Epsilon (ε): the additive epsilon indicator measures the minimum translation factor by which every point in the reference front can be ε-dominated by the approximation. Lower ε values indicate better convergence of the approximation toward the reference front [47].

RMSD: for each benchmark complex with a known crystallographic binding pose, the RMSD between the predicted and experimental ligand heavy-atom positions is computed. An RMSD below 2.0 Å is conventionally considered a successful docking prediction [10], [11]. RMSD analysis bridges the gap between algorithmic performance indicators and the practical objective of reproducing experimentally observed binding modes.

MO-WAO is implemented in Java within the jMetal 5.0 multi-objective optimization framework [48], which provides standardized implementations of NSGA-II, SMPSO, GDE3, MOEA/D, and SMS-EMOA, ensuring consistent algorithm implementations and fair comparisons. The jMetal framework offers modular components for encoding, operators, archiving, and quality indicator computation. Energy evaluation is performed using the AutoDock 4.2 scoring engine [10]. Protein and ligand structures are prepared using AutoDockTools 1.5.6 [30]: polar hydrogens are added, Gasteiger charges are computed and assigned, and non-polar hydrogens are merged. Ligand rotatable bonds are identified and set using the torsion tree mechanism. Grid maps for each atom type are pre-computed using AutoGrid 4.2 with a spacing of 0.375 Å, centered on the experimentally determined ligand binding site with sufficient margin to accommodate alternative poses [10], [30].

Statistical analysis is conducted at two levels. First, the Wilcoxon signed-rank test [49] with a significance level of α = 0.05 is applied for pairwise comparison of MO-WAO against each competitor algorithm on each quality indicator. Second, the Friedman test [49] is used to obtain an overall ranking of all algorithms across the benchmark suite, providing a non-parametric assessment of relative algorithm performance. All statistical tests are performed using 30 independent runs per algorithm per instance. Experiments are conducted on a workstation equipped with an Intel Core i7 processor, 32 GB RAM, running Ubuntu 20.04 LTS. The computational time for each algorithm run is recorded to assess practical efficiency.

The HV indicator provides a comprehensive measure of both convergence and diversity of the Pareto front approximation. Table 3 presents the median HV values (and interquartile ranges) obtained by each algorithm across 30 independent runs on each benchmark instance. Higher HV values indicate superior performance. The HV analysis is expected to reveal performance differences across complexes of varying flexibility. For rigid to moderately flexible ligands (e.g., 1STP, 1G5S, 1D4K with ≤ 6 rotatable bonds), all algorithms are anticipated to achieve comparable HV values, as the lower-dimensional search spaces are more amenable to thorough exploration within the allocated evaluation budget. For highly flexible ligands (e.g., 1HBV with 19 rotatable bonds, 1BDR and 1HTF with 14 rotatable bonds), differences among algorithms should become more pronounced, as the effectiveness of the exploration–exploitation balance becomes critical in high-dimensional spaces. The hydrogen bonding phase of MO-WAO, which provides strong directed exploitation through archive-guided attraction, combined with the diversity maintenance afforded by the evaporation phase, is hypothesized to offer advantages on these challenging instances. Wilcoxon signed-rank tests should be applied to determine the statistical significance of observed differences at the α = 0.05 level.

The IGD provides a complementary assessment of convergence and diversity. Table 4 presents the median IGD values for all algorithms. Lower IGD values indicate better performance.

The IGD analysis is expected to complement the HV results by providing a reference front-based assessment. Algorithms that achieve low IGD values produce Pareto fronts that closely approximate the reference front constructed from the union of all non-dominated solutions across all algorithms and runs. Consistent rankings between HV and IGD would strengthen confidence in the relative performance assessment. Discrepancies between HV and IGD rankings may arise from differences in reference point sensitivity (HV) versus reference front construction (IGD) and should be discussed in the context of each benchmark instance.

The spread indicator evaluates the distribution uniformity of solutions along the Pareto front, independent of convergence quality. Table 5 presents the median spread values. Lower values indicate better distribution uniformity.

Good spread performance is particularly important in the drug design context, as a well-distributed Pareto front provides medicinal chemists with a comprehensive spectrum of binding mode alternatives, from conformations with maximally favorable intermolecular interactions (potentially strained) to relaxed conformations with moderate binding. The crowding distance-based archive management and leader selection strategies in MO-WAO are specifically designed to promote distribution uniformity, and their effectiveness should be evaluated through the spread analysis.

The epsilon indicator provides a convergence-focused assessment, measuring the minimum additive factor needed to translate the reference front to be dominated by the approximation. Table 6 reports the median epsilon values. Lower values indicate better convergence.

The epsilon indicator isolates convergence quality from diversity considerations, providing a focused measure of how close the approximation comes to the reference front in the worst case. Strong epsilon performance by MO-WAO would indicate that the hydrogen bonding phase’s exploitation pressure, combined with the particle motion phase’s directed search toward archive leaders, effectively drives solutions toward the Pareto front. Algorithms with strong convergence (low epsilon) but poor diversity (high spread) may indicate premature convergence issues.

Visual inspection of the Pareto fronts produced by each algorithm provides qualitative insight into the shape, extent, and distribution of solutions in the objective space. For each selected benchmark complex, the Pareto front approximation is plotted in the (f1, f2) = (intermolecular energy, intramolecular energy) plane. Representative complexes for visualization include:

I. 1HVR (HIV-1 protease, XK263): A moderately flexible complex representing a therapeutically important target–ligand pair.

II. 1HPV (HIV-1 protease, KNI-272): Another HIV-1 protease complex with a different inhibitor scaffold, enabling within-target comparison.

III. 1HTF (HIV-1 protease, GR126045): A flexible ligand (14 rotatable bonds) presenting greater search difficulty. IV. 1STP (Streptavidin, biotin): A rigid ligand (5 rotatable bonds) representing the easiest benchmark instance. [Figures to be generated from experimental runs. Each figure should display the Pareto fronts of all six multi-objective algorithms overlaid on the same axes, with distinct markers and colors for each algorithm. The x-axis represents intermolecular energy (kcal/mol) and the y-axis represents intramolecular energy (kcal/mol). Inset panels or supplementary figures showing the convergence of the front over function evaluations may be included.] The Pareto front visualizations are expected to reveal characteristic differences among the algorithms. Well-performing algorithms should produce fronts that extend across a wide range of the intermolecular– intramolecular energy trade-off space, with solutions uniformly distributed along the front. Algorithms with strong convergence but poor diversity may produce fronts that cluster in a narrow region, while algorithms with good diversity but weak convergence may produce widely spread but suboptimal fronts. The shape of the front itself provides insight into the nature of the trade-off for each complex: concave fronts suggest complementary objectives, convex fronts indicate competing objectives, and discontinuous fronts may reveal separate binding modes.

To assess the practical utility of the multi-objective docking results, the RMSD between the best predicted binding pose and the crystallographic reference pose is computed for each algorithm and each benchmark complex. Table 7 presents the best RMSD values achieved across 30 independent runs.

For the multi-objective algorithms, the best RMSD is selected from the Pareto front solution that minimizes the distance to the crystallographic pose. This analysis demonstrates that the multi-objective formulation does not sacrifice binding pose accuracy for the sake of providing diverse trade-off solutions. In many cases, Pareto fronts contain solutions with RMSD values competitive with or superior to the single-objective LGA result, as the multi-objective search may discover conformations that a single-objective approach misses due to its focus on minimizing the scalar sum of energies [20], [21]. A detailed analysis of specific protein–ligand interactions hydrogen bonds, hydrophobic contacts, and π–π stacking for the best-RMSD poses from MO-WAO should be compared against the crystallographic reference to validate the physical plausibility of the predicted binding modes.

The convergence behavior of MO-WAO is analyzed through the evolution of the HV indicator as a function of the number of function evaluations consumed. Convergence curves for representative complexes (e.g., 1HVR, 1HBV, 1STP) illustrate the search dynamics of each algorithm over time.

[Convergence curve figures to be generated from experimental runs. Each figure should display HV vs. function evaluations (0–25,000) for all algorithms on a single instance.] The convergence analysis is expected to reveal the characteristic search dynamics of MO-WAO. During the early iterations, the evaporation phase should promote broad exploration of the docking landscape, preventing premature convergence to a narrow region of the objective space. As the search progresses and the inertia weight decreases, the hydrogen bonding phase should increasingly dominate, providing focused exploitation of the most promising Pareto front regions. The particle motion phase contributes to both exploration and exploitation throughout the search, with its balance shifting as the inertia weight decreases. This three-phase dynamic is expected to produce a convergence curve that initially rises rapidly (broad exploration discovers many non-dominated solutions) and then transitions to a gradual, sustained improvement (focused exploitation refines the front), avoiding the premature stagnation observed in algorithms with insufficient diversity maintenance.

The statistical significance of performance differences between MO-WAO and each competitor is assessed using the Wilcoxon signed-rank test at the α = 0.05 significance level. Table 8 summarizes the pairwise comparison results for the HV indicator.

as a drug target for antiretroviral therapy.

sacrifice some binding contacts.

vs. GDE3 vs. MOEA/D vs. SMS-EMOA [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [p-value, symbol] [X/Y/Z] [X/Y/Z] [X/Y/Z]

Spread Rank Epsilon Rank Average Rank [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK] [RANK]

The Friedman test is applied to obtain an overall ranking of the algorithms across all eleven benchmark instances for each quality indicator. Table 9 summarizes the Friedman rankings.

The statistical analysis should reveal the instances and indicators where MO-WAO achieves statistically significant advantages or disadvantages relative to each competitor. A comprehensive discussion should address: 1) the number of instances where MO-WAO significantly outperforms, ties with, or is outperformed by each algorithm, 2) whether performance advantages are correlated with ligand flexibility (number of rotatable bonds) or other problem characteristics, 3) the overall Friedman ranking and its consistency across indicators, and 4) the practical significance of observed statistical differences in the context of drug design

To illustrate the practical utility of the MO-WAO approach in a drug design context, a detailed case study is conducted on a therapeutically relevant target. The HIV-1 protease complex 1HVR, featuring the inhibitor XK263, is selected due to the well-characterized binding mode and the clinical importance of HIV-1 protease

The Pareto front obtained by MO-WAO for 1HVR is analyzed in terms of its drug design implications. Solutions along the front represent different compromises between intermolecular binding strength (f1) and intramolecular conformational strain (f2). At one extreme, conformations with the most favorable intermolecular energy achieve strong receptor interactions but may impose significant strain on the ligand. At the other extreme, conformations with minimal intramolecular strain adopt relaxed geometries that may From a medicinal chemistry perspective, the Pareto front provides valuable information for lead optimization decision-making. Solutions near the “knee” of the Pareto front, the region of maximum marginal trade-off between the two objectives, are of particular interest, as they represent conformations that achieve a near-optimal balance between binding affinity and conformational feasibility. These knee-region solutions often correspond to binding modes that are both energetically favorable and geometrically accessible, making them the most promising candidates for further optimization through Structure–Activity Relationship (SAR) studies.

Furthermore, the diversity of solutions along the Pareto front may reveal alternative binding modes that differ in their hydrogen bonding networks, hydrophobic contacts, or orientation within the binding site. These alternative modes can suggest novel chemical modifications to improve selectivity, reduce off-target effects, or overcome drug resistance mutations. In the context of HIV-1 protease, where drug resistance is a major clinical challenge, the availability of multiple binding mode hypotheses from a single Pareto-based docking run is particularly valuable for designing inhibitors that remain effective against resistant viral variants. [Detailed interaction analysis, including hydrogen bond donor–acceptor distances, hydrophobic contact surfaces, and comparison with crystallographic data, to be completed with experimental results.]