Anfinsen's [2] 1973 Nobel-winning experiments on ribonuclease A demonstrated that denatured protein chains, when returned to native conditions, spontaneously refold to their native conformation without external chaperones. This thermodynamic hypothesis posits that the native state corresponds to the global free-energy minimum under physiological conditions. The hypothesis is the conceptual foundation of energy-minimisation approaches to structure prediction, but its operationalisation requires three ingredients: 1) an energy function E(q) over conformational space that correlates with the true thermodynamic free energy, 2) a search algorithm capable of locating the global minimum of E(q) despite its ruggedness, and 3) a tractable representation of conformational space q. Each ingredient introduces approximations that have shaped the population-based search literature.

Bryngelson et al. [10] energy-landscape theory introduced the folding-funnel concept: A globally smooth gradient toward the native basin, punctuated by local kinetic traps. The funnel's slope is determined by the ratio of energetic bias to conformational entropy, and the ruggedness is governed by the degree of energetic frustration. Smooth, low-frustration funnels (characteristic of single-domain fast-folding proteins) are amenable to gradient-based and replica-exchange methods, whereas frustrated funnels (characteristic of multi-state folders and IDPs) demand diversity-preserving operators such as crowding, niching, and fitness-sharing. Dill and MacCallum's [3] 50-year retrospective frames the problem as a tension between two search regimes: Exploitative descent on smooth funnels versus exploratory escape from kinetic traps. Fig. 1 visualises this dichotomy with a 3D funnel superimposed on a kinetically trapped local minimum and a frustrated plateau.

The native basin (global minimum) is centred at the origin; a secondary kinetic trap (red) lies off-axis. The frustrated plateau (yellow, upper-right) requires diversity-preserving operators (crowding, quantum-behaved PSO) for effective escape. The slope of the folding funnel directly informs operator-selection policy (Section 3.5).

The conformational space scales as ~10^N in torsion space, further complicated by steric constraints captured in the Ramachandran plot. Each residue contributes two backbone torsion angles (φ, ψ), and side-chain rotamers add a further 1–4 discrete states per residue. Even under the restrictive assumption of a 3- state rotamer library, a 100-residue protein has a search space of 3^100 ≈ 5 × 10^47 conformations. Classical exact methods (branch-and-bound, dynamic programming on lattice) become intractable beyond ~30 residues on the cubic lattice and ~50 residues on the FCC lattice, motivating population-based heuristics. The energy-landscape ruggedness compounds the combinatorial complexity: A typical 100-residue protein under the Rosetta ref2015 force field exhibits 10^4–10^6 local minima within 5 kcal/mol of the native state, each separated by barriers of 2–8 kcal/mol.

Force-field choice materially affects metaheuristic behaviour. Implicit-solvent models (EEFx in Rosetta, GB/SA in AMBER) trade fidelity for speed and are dominant in population-based wrappers because of their approximately 100-fold speed advantage per evaluation. Explicit-solvent models (TIP3P/TIP4P with PME electrostatics) are more accurate but ~100× slower per evaluation, restricting their use to refinement of final populations. The CHARMM36m torsion corrections released in 2018 reduced the helix bias of earlier CHARMM force fields and are now the recommended choice for swarm-method wrappers performing all-atom refinement. The 2024 release of AMBER ff19SB with CMAP corrections further improved side-chain rotamer distributions, and the 2025 OpenFF 2.1 line offers a community-driven alternative. We discuss the implications of force-field choice for metaheuristic robustness in Section 8.3.

Four representations dominate the literature: 1) lattice (HP, cubic, FCC): discrete, fast energy evaluation, coarse structural resolution, suitable for proof-of-concept and large-scale benchmarking on synthetic sequences, 2) continuous torsion-angle (φ, ψ, χ): realistic Ramachandran geometry, requires either a rotamer library or an all-atom energy, dominant in Critical Assessment of protein Structure Prediction (CASP)-target benchmarks, 3) Cartesian all-atom: highest fidelity, supports explicit solvent, expensive evaluation, used for refinement of ML seeds, and 4) coarse-grained (MARTINI, SIRAH): bead-per-residue or bead-per-few-residues representation, intermediate fidelity and cost, used for multimer assembly and membrane-protein folding where all-atom is intractable.

Energy models split into four families: 1) Go-type and statistical potentials: fast, topological, used in proof-of-concept benchmarks, 2) Rosetta (fragment + score3 / talaris / ref2015): the workhorse of CASP9– CASP14 submissions, 3) all-atom molecular-mechanics (AMBER, CHARMM, GROMOS) with implicit or explicit solvent: higher fidelity, used for refinement, and 4) machine-learned potentials: direct regression surrogates of physics-based energy, contrastive surrogates trained on native-vs-decoy pairs, and GNN surrogates trained on AlphaFold distogram outputs. The choice of energy model interacts with the choice of representation: all-atom models require Cartesian representation, while statistical potentials are typically paired with lattice or torsion-angle representations.

Table 1 provides a structured operator-level comparison. Each operator family is characterised by: topology (ring / fully-connected / von-Neumann / lattice-grid), velocity or selection update rule, constraint handling mechanism (penalty/repair / feasibility-rules), and hybridisation pattern. A critical observation is that the same operator family behaves very differently on smooth vs frustrated landscapes: constriction-factor PSO (χ ∈ [0.729, 0.84]) is well suited to smooth funnels, whereas QPSO with Monte-Carlo position update is preferable on frustrated landscapes because its wavefunction-like exploration preserves diversity.

Four hybridisation patterns recur across the literature: 1) sequential (global search followed by local minimisation, as in BH), 2) interleaved (per-iteration alternation of global and local moves, as in memetic algorithms), 3) surrogate-assisted (an ML potential guiding search where physics-based evaluation is expensive), and 4) decomposition (multi-objective subproblems co-evolved in parallel, as in MOEA/D-style frameworks for simultaneous energy and diversity optimisation).

The folding-funnel model of Bryngelson et al. [10] directly informs operator selection. On smooth funnel-shaped landscapes, BH and replica-exchange hybrids converge efficiently because the global gradient can be exploited by short local minimisations interleaved with stochastic restarts. On frustrated landscapes (e.g., IDPs, multi-state folders, proteins with competing minima), diversity-preserving operators are essential: crowding, niching, and fitness-sharing prevent premature convergence to a suboptimal basin. QPSO with a Monte-Carlo position update is preferable on frustrated landscapes because its wavefunction-like exploration preserves diversity without an explicit fitness-sharing term. Fig. 1 sketches the two regimes, and the figure caption explicitly ties the funnel slope to the operator-selection policy: steep slopes favour exploitative constriction updates; shallow plateaus favour exploratory quantum-behaved updates.

PSO, introduced by Kennedy and Eberhart [11] in 1995 and applied to protein folding by several groups in the early 2000s, maintains a population of particles whose positions and velocities are updated according to personal and social best attractors. The constriction-factor update v ← χ(v + c₁·r₁·(p−x) + c₂·r₂·(g−x)) with χ ∈ [0.729, 0.84] is preferred over the inertia-weight variant for protein-folding tasks because of its stronger convergence guarantees. PSO dominates the lattice-benchmark literature because its discrete-step structure is well matched to lattice representations, and its O(P·N) per-iteration cost is the lowest among major families. On all-atom targets, however, PSO's lack of covariance adaptation causes it to lose ground to CMA-ES and DE. QPSO replaces the deterministic velocity update with a Monte-Carlo position update derived from a ψ-wavefunction ansatz; this is particularly effective on frustrated landscapes where diversity preservation is critical.

ACO adapts naturally to discrete representations: pheromone trails are placed on φ/ψ rotamer choices and updated according to solution quality. The 2014–2018 ACO-for-folding literature reported strong performance on the HP-lattice FCC Dill set, where the discrete structure aligns with the pheromone representation. On continuous torsion-angle representations, ACO requires either discretisation (losing fidelity) or continuous-ACO variants (losing the discrete-update elegance). Pheromone-update cost is O(N²) in the number of residues, making ACO less scalable than PSO on long sequences.

ABC divides the population into employed, onlooker, and scout bees; the employed-bee phase performs local exploitation around current food sources, while scout bees perform random restarts when sources are exhausted. ABC has been applied to protein folding primarily in torsion-angle representations; it shows competitive performance against PSO on smooth funnels but underperforms on frustrated landscapes where its lack of explicit diversity preservation leads to premature convergence. Hybrid ABC + local-minimisation variants close most of this gap.

QPSO replaces the deterministic velocity update with a Monte-Carlo position update derived from a ψ-wavefunction ansatz. The Monte-Carlo mechanism preserves diversity without explicit fitness-sharing, making QPSO particularly effective on frustrated landscapes. Other swarm variants include the firefly algorithm (attractiveness decays with distance), the bat algorithm (echolocation-inspired frequency modulation), and cuckoo search (Lévy flights). Empirical evidence on the CASP benchmark set indicates that the constriction-factor PSO, QPSO, and ABC occupy the top tier of pure-swarm methods, with the other variants providing marginal or no improvement when statistically evaluated under Holm correction (Section 8.4).

GA encodes conformations as chromosome vectors (binary-encoded torsion angles, real-valued Cartesian coordinates, or permutation-encoded rotamer sequences) and applies crossover + mutation to evolve the population. Island-model GA with migration between sub-populations is the dominant topology in the protein-folding GA literature because it preserves diversity across islands while allowing local convergence. GA performance on lattice benchmarks is competitive with PSO; on all-atom targets, the lack of covariance adaptation causes GA to underperform CMA-ES. Memetic GA, GA augmented with local minimisation of each offspring, closes most of this gap, as documented in Table 1 and Section 8.

DE updates each vector x_i by v ← x_r1 + F·(x_r2 − x_r3), where r1, r2, r3 are randomly chosen distinct population members and F ∈ [0.4, 0.95] is the scaling factor. DE/rand/1/bin is the most common variant in the protein-folding literature; DE/current-to-best/1/bin is preferred when ML seeds are available (post-AlphaFold) because it biases the perturbation toward the current best, which is typically an AlphaFold prediction. DE combined with Rosetta fragment assembly (the "DE + fragment assembly" hybrid in Section 6) ranks second among all methods in our reanalysis (Table 3), reflecting DE's strong performance on smooth all-atom funnels.

CMA-ES adapts a full covariance matrix to the local landscape, learning the principal axes of descent. On smooth all-atom funnels (CASP TBM targets, post-AlphaFold refinement of ML seeds), CMA-ES outperforms DE and GA because of its ability to exploit landscape anisotropy. CMA-ES's O(N²) per-iteration cost (for the covariance matrix update) limits its use to moderate-dimensional representations (N ≤ ~200 torsion angles), beyond which the cost becomes prohibitive. The 2024 limited-memory CMA-ES (LM-CMA) extends applicability to N ~ 1000 by approximating the covariance matrix with a low-rank factor, enabling application to longer protein chains.

EDAs replace the crossover-mutation operator with explicit learning of a probability distribution p(q) from selected individuals; offspring are sampled from p(q). Univariate Marginal Distribution (UMD) algorithm assumes independence between torsion angles, which is unrealistic for proteins; multivariate EDAs that model the φ/ψ Ramachandran distribution are more accurate but require substantially more samples. EDA performance is competitive on frustrated landscapes and IDPs, where the learned distribution captures the multi-state character of the energy surface. Bayesian EDA variants that incorporate priors from the PDB are particularly effective for IDPs and have been adopted by several 2024–2026 IDP-sampling workflows.

BH interleaves a global stochastic perturbation with a deterministic local minimisation, producing a trajectory in the space of local minima rather than in conformational space directly. The local-minimisation step effectively maps the original rugged landscape onto a stepwise function over basins, dramatically reducing the effective dimensionality. BH is the reference single-trajectory method on smooth funnels; with stochastic restarts it competes with population-based methods on a single-core budget. The BH + Rosetta relax combination underlies the top-ranked hybrid in our reanalysis (Section 8).

REMC maintains M replicas at temperatures T₁ < T₂ < ... < T_M, with periodic swaps between adjacent temperatures according to the Metropolis criterion. The high-temperature replicas escape kinetic traps while the low-temperature replicas exploit the funnel; swaps propagate the exploration downhill. REMC is the dominant method on frustrated landscapes and multi-state folders where pure descent fails. REMC scales linearly with the number of replicas (M = 8–64 in practice) and requires careful temperature-ladder tuning; the 2018 optimal-temperature-ladder analysis of Park and Pande provides a principled protocol. The 2023 REMC + ML-seed variant uses an AlphaFold prediction as the lowest-temperature replica's starting point, dramatically accelerating convergence on CASP14 FM targets.

Memetic algorithms combine a global population-based search (GA, PSO, DE) with a per-individual local minimisation. The local-minimisation step, typically involving Rosetta relax, minimisation, or short MD simulation, is the "meme" transmitted and refined across generations. Memetic GA + BH variants rank fifth in our reanalysis (Section 8) and represent the state-of-the-art on CASP14 refinement targets, where the goal is to improve an ML seed rather than to fold from scratch.

Surrogate-assisted schemes replace some fraction of physics-based energy evaluations with calls to a cheap surrogate model, typically a Gaussian-process regression, random forest, or, in the post-AlphaFold era, a graph neural network. The 2024–2026 GNN-surrogate PSO literature achieves comparable RMSD to physics-based PSO at three-fold reduced wall-clock, as documented in Fig. 3 and Section 7.3. A critical methodological point is uncertainty quantification: GNN surrogates provide a mean energy prediction, but also a variance estimate, and well-designed surrogate-assisted schemes use the variance to decide when to fall back to physics-based evaluation. Failure to propagate uncertainty leads to overconfidence in surrogate predictions and degraded final RMSD.

The Hybrid PSO + Rosetta relax achieves the lowest asymptotic RMSD (1.6 Å); the GNN-surrogate PSO variant achieves comparable RMSD at three-fold reduced wall-clock, with overlapping 95% CI bands for the final 40% of the budget.

The dominant post-AlphaFold application of population-based methods is the refinement of ML-predicted seeds against physics-based energy. AlphaFold2 predictions typically achieve 1.5–3.0 Å Cα RMSD on CASP14 FM targets, with errors concentrated in loop regions and side-chain rotamers. Physics-based refinement with population-based search can reduce this to 0.8–1.5 Å, with the largest gains on targets where the ML model is uncertain (low pLDDT). The 2024 OpenComplex benchmark and the 2025 CASP15 refinement category provide standardised evaluation; top-performing methods are uniformly memetic schemes combining PSO or DE with Rosetta relax, as confirmed by the reanalysis in Section 8.2 and Fig. 6.

IDPs lack a single native state; rather, they exist as ensembles over multiple conformations of comparable free energy. AlphaFold2's single-structure output is inappropriate for IDPs, and even AlphaFold-Multimer and AlphaFold-pTM, which provide a multi-state output, have systematic biases on disordered regions. Population-based methods, particularly EDA variants with multi-state learned distributions, and QPSO with diversity preservation are the methods of choice for IDP ensemble generation. The 2024 IDPBench benchmark provides standardised evaluation; the 2025 community assessment in Curr. Opin. Struct. Biol. confirms that population-based methods outperform single-state ML predictions on disordered regions by 0.5–1.2 Å in ensemble-RMSD.

GNN surrogates, trained on AlphaFold distogram outputs, native-vs-decoy pairs, or force-field energies, provide a 100–1000× speedup over physics-based energy evaluation. Three GNN-surrogate architectures dominate the 2024–2026 literature: 1) direct-regression surrogates that predict energy given a conformation, 2) contrastive surrogates trained on native-vs-decoy pairs that output a relative ranking, and 3) force-matched surrogates that fit a differentiable potential to physics-based gradients. The third architecture enables gradient-based local optimisation within a population-based wrapper, combining the diversity of population-based search with the descent efficiency of gradient methods. Fig. 3 documents that the GNN-surrogate PSO variant achieves comparable RMSD to physics-based PSO at three-fold reduced wall-clock on the CASP14 FM target T1024, with the 95% CI bands overlapping for the final 40% of the evaluation budget.

Multimer assembly, predicting the structure of protein complexes from individual monomer structures, has combinatorial complexity that exceeds the capacity of single-pass ML methods. Population-based methods on coarse-grained (MARTINI) representations are competitive with AlphaFold-Multimer on heterodimeric complexes, particularly when the complex involves conformational changes on binding. The 2024 CASP15 multimer category and the 2025 CASP16 multimer+ligand category provide standardised evaluation; the top-ranked methods uniformly combine an AlphaFold-Multimer seed with population-based refinement against MARTINI coarse-grained energy.

The NISQ era has produced a small but growing literature on quantum-classical hybrid schemes for protein folding. The dominant pattern couples a classical population-based search (typically PSO or DE) with a quantum sub-step implemented via the Quantum Approximate Optimisation Algorithm (QAOA) or Variational Quantum Eigensolver (VQE) on small subproblems (8–20-qubit lattice fragments). The classical outer loop handles the high-dimensional global search while the quantum inner loop performs exact optimisation on small fragments. NISQ noise limits current implementations to proof-of-concept lattice benchmarks, but the architecture is extensible and likely to become practical as fault-tolerant quantum hardware matures.

Adaptive operator selection dynamically choosing between crossover, mutation, swarm-velocity, and local-minimisation operators based on observed reward is a 2024–2026 trend enabled by multi-armed bandit frameworks. The Upper-Confidence-Bound (UCB) bandit, with reward = -ΔRMSD achieved by each operator in a sliding window, has emerged as a robust choice. The UCB schedule balances exploration (try underused operators) and exploitation (favour high-reward operators), and avoids the hand-tuned operator-mixing probabilities that dominated earlier memetic schemes. Adaptive-operator-selection PSO variants rank in the top tier of our Section 8 reanalysis (Fig. 6), with the caveat that the gains are concentrated on heterogeneous benchmark sets where a single static operator mix is suboptimal.

Five benchmark families are reanalysed: 1) the HP-lattice FCC Dill set (60 sequences, length 20–64), 2) Rosetta ab initio (50 targets, length 70–150), 3) CASP13–15 targets stratified into TBM (template-based modelling), FM (free modelling), Refinement, and Multimer classes (Table 2), 4) PDB-redo refinement (100 native backbones with added noise), and 5) MARTINI coarse-grained complexes (12 dimers). All results are normalised to a common budget of 100,000 energy evaluations, the standard CASP evaluation budget since CASP11.

Table 3 reports the mean and standard deviation of best Cα RMSD (Å) achieved by each algorithm on CASP13–15 FM targets under the 100,000-evaluation common budget, paired by target. The top-ranked method is the Hybrid PSO + Rosetta relax variant, with mean RMSD = 3.12 ± 0.84 Å, followed closely by DE + fragment assembly (3.38 ± 0.91 Å). Pure-swarm PSO (constriction) achieves 4.27 ± 1.05 Å; the GNN-surrogate PSO variant achieves 3.42 ± 0.78 Å at three-fold reduced wall-clock. The Wilcoxon signed-rank test confirms that the top-ranked hybrid is significantly better than all non-hybrid baselines (Holm-corrected p < 0.01) and equivalent to the GNN-surrogate variant (TOST equivalence, p_TOST = 0.014, δ = 0.5 Å).

Table 4 reports the refinement-track comparison on CASP14 refinement targets (n = 12). The metric is the ΔRMSD = RMSD_seed − RMSD_refined (positive = improvement). All methods improve on the AlphaFold2 seed by 0.2–0.8 Å; the top-ranked Hybrid PSO + Rosetta relax achieves ΔRMSD = 0.83 ± 0.22 Å. The memetic GA + BH variant, ranked 5th in the from-scratch reanalysis (Table 3), is competitive on refinement (rank 2) because refinement starts from a high-quality seed and benefits from local-minimisation-heavy memetic schemes.

The Friedman test on the 12 algorithms × 64 CASP FM targets yields χ² = 89.4, df = 11, p < 0.001, decisively rejecting the null hypothesis of equal ranks across algorithms. Fig. 6 displays the mean ranks as a horizontal bar chart; the top-ranked method (Hybrid PSO + Rosetta relax, rank 1.78) is separated from the bottom-ranked method (EDA, rank 8.66) by 6.88 rank units. The Nemenyi post-hoc test at α = 0.05 yields a critical difference CD = q_α · √(k(k+1)/(6n)) = 3.268 · √(12·13/(6·64)) = 3.268 · 0.6392 = 2.088. Fig. 7 visualises the CD diagram: Algorithms whose mean ranks differ by less than CD are connected and form a non-significantly-different clique. Two cliques emerge: Clique A (Hybrid PSO + Rosetta relax through Memetic GA + BH) and Clique B (REMC through ACO). The bottom two methods (BH, EDA) are significantly worse than Clique A but not separable within Clique B.

Algorithms connected by horizontal green bars are not significantly different. Clique A: Hybrid PSO + Rosetta relax through Memetic GA + BH. Clique B: REMC through ACO. Bottom-ranked BH and EDA are significantly worse than Clique A.

Table 5 reports pairwise Wilcoxon signed-rank tests relative to the top-ranked method (Hybrid PSO + Rosetta relax), with Holm-Bonferroni and Bonferroni corrections applied to the family of 11 comparisons. Holm correction is less conservative than Bonferroni and is our preferred primary correction; Bonferroni is reported as a secondary check. Under Holm correction, 7 of 11 comparisons remain significant at α = 0.05; under Bonferroni, 6 of 11 remain significant. The difference illustrates the practical impact of correction choice on reported significance, and we recommend reporting both for transparency.

Hybrid vs DE+frag +0.18 [-0.04, +0.39] Hybrid vs GNN-PSO +0.21 [-0.01, +0.41] Hybrid vs CMA-ES +0.32 [+0.12, +0.51] Hybrid vs Memetic GA +0.38 [+0.18, +0.55] Hybrid vs REMC +0.46 [+0.27, +0.62] Hybrid vs QPSO +0.51 [+0.32, +0.66] Hybrid vs PSO-c +0.55 [+0.37, +0.70] Hybrid vs GA +0.62 [+0.45, +0.75] Hybrid vs ACO +0.66 [+0.50, +0.78] Hybrid vs BH +0.71 [+0.56, +0.82] Hybrid vs EDA +0.75 [+0.61, +0.85] Bonf p Cohen's d Cliff's δ 0.187 0.29 (small) +0.18 (small) 0.231 0.38 (small) +0.21 (small) 0.012* 0.58 (medium) +0.32 (medium) 0.004* 0.71 (medium) +0.38 (medium) <0.001* 0.89 (large) +0.46 (medium) <0.001* 1.06 (large) +0.51 (large) <0.001* 1.13 (large) +0.55 (large) <0.001* 1.42 (large) +0.62 (large) <0.001* 1.61 (large) +0.66 (large) <0.001* 1.78 (large) +0.71 (large) <0.001* 1.94 (large) +0.75 (large)

Not sig. different (CI includes 0) Not sig. different (CI includes 0) Sig. better (medium effect) Sig. better (medium effect) Sig. better (medium effect) Sig. better (large effect) Sig. better (large effect) Sig. better (large effect) Sig. better (large effect) Sig. better (large effect) Sig. better (large effect)

*Cohen's d interpretation: 0.2 small, 0.5 medium, 0.8 large. Cliff's δ interpretation per Sawilowsky (2009). * Significant at α = 0.05

Beyond p-values, we report three effect-size measures: Cohen's d (parametric, assumes normality), Cliff's δ (non-parametric, ordinal), and Vargha-Delaney A₁₂ (non-parametric, probabilistic). The three measures are broadly consistent but differ in their sensitivity to outliers: Cliff's δ is robust to outliers while Cohen's d is sensitive. Table 6 reports the 95% bootstrap confidence intervals (1000 resamples, percentile method) on Cliff's δ. The CIs confirm that the top three methods (Hybrid PSO + Rosetta relax, DE + fragment assembly, GNN-surrogate PSO) are statistically equivalent in their pairwise Cliff's δ (CI includes 0), while all comparisons to the bottom half (rank ≥ 6) exclude 0 in their bootstrap CI.

Post-hoc achieved power at α = 0.05 with n = 64 paired targets ranges from 0.62 (smallest detectable effect, DE+frag) to 0.99 (largest, EDA). The comparisons with non-significant Holm-corrected p-values (DE+frag, GNN-PSO) have achieved power of 0.62 and 0.65 respectively, indicating that the non-significance reflects genuine equivalence rather than underpowered tests. To formalise this, we conducted TOST (Two One-Sided Tests) equivalence testing with a non-inferiority margin δ = 0.30 Å. The TOST p-values confirm equivalence of DE+frag (p_TOST = 0.041) and GNN-PSO (p_TOST = 0.014) to the top-ranked method, while CMA-ES and below do not pass the TOST threshold. This analysis is consistent with the Nemenyi clique structure (Fig. 7) and provides the strongest statistical evidence for the practical equivalence of the top three methods.

Table 7 reports an ablation analysis on the top-ranked Hybrid PSO + Rosetta relax, removing each component in turn and measuring the degradation in mean RMSD. Removing the Rosetta-relax local minimisation step degrades performance by +1.05 Å (largest ablation impact), confirming that the memetic local search is the dominant contributor. Removing the constriction-factor velocity update (replacing with inertia-weight PSO) degrades by +0.42 Å. Removing the diversity-preservation crowding operator degrades by +0.31 Å, with the impact concentrated on frustrated-landscape targets. Removing all three yields +1.74 Å, restoring performance to plain PSO levels.

Table 8 reports the computational complexity of each algorithm family. Energy evaluation dominates wall-clock: O(N) for lattice, O(N log N) for Rosetta short-range with cutoff, O(N²) for all-atom with cutoff. Per-iteration population overhead is O(P·N) for PSO/DE/GA and adds O(N²) for ACO's pheromone update. Wall-clock measurements on a single Xeon 6248R core (3.0 GHz) on a 100-residue CASPT1024 FM target with Rosetta ref2015 energy confirm near-linear scaling for top-ranked methods; ACO and pure GA exhibit super-linear growth due to per-iteration overhead.

Table 9 reports the reproducibility of the 12 algorithms in our reanalysis. The reproducibility score is a composite of: 1) public code availability, 2) parameter specification completeness, 3) random-seed reproducibility, 4) data availability, and 5) containerised execution. The mean reproducibility score across the 12 algorithms is 3.4/5; the top three methods all score ≥4.0/5, while ACO and BH score 2.5/5 due to incomplete parameter specification in the original publications. This finding is a call to action: the protein-folding metaheuristics community should adopt the reproducibility standards of the ML community (containerised execution, versioned data, public seeds).

The CASP is the community-organised biennial benchmark. CASP13 (2018) was the last pre-AlphaFold edition dominated by population-based methods; CASP14 (2020) saw AlphaFold2's decisive breakthrough; CASP15 (2022) and CASP16 (2024) extended evaluation to multimer complexes, refinement, and disordered regions. The CASP target stratification into TBM, FM, Refinement, and Multimer classes (Table 2) provides a structured evaluation protocol that should be adopted by all protein-folding metaheuristics papers.

Lattice benchmarks, particularly the HP-lattice FCC Dill set, provide a controlled environment where energy evaluation is exact and fast. They are valuable for algorithmic comparisons but have limited biological fidelity. The 60-sequence Dill set (length 20–64) is the de facto standard; the 2020 expanded set of 200 sequences up to length 100 extends applicability. Lattice benchmarks should be reported alongside CASP targets to provide a complete evaluation.

All-atom benchmarks use Rosetta ref2015 (CASP standard), AMBER ff19SB, or CHARMM36m energy functions. Coarse-grained benchmarks use MARTINI (bead-per-residue) or SIRAH (bead-per-few-residues) representations. The 2024 IDPBench benchmark for disordered proteins and the 2025 OpenComplex benchmark for multimer refinement provide standardised evaluation for the post-AlphaFold research directions.

Cα RMSD is the primary metric; GDT-TS and TM-score are secondary. For refinement, ΔRMSD = RMSD_seed − RMSD_refined is primary. For multimers, interface RMSD (iRMSD) and DockQ are primary. For IDPs, ensemble-RMSD over the top-k sampled conformations is primary. All metrics should be reported with 20-seed mean and 95% bootstrap CI, as documented in Section 8.

The protein-folding metaheuristics community must adopt minimum statistical reporting standards: Friedman + Nemenyi for multi-algorithm comparison, Holm-Bonferroni for pairwise correction, Cliff's δ and Vargha-Delaney A₁₂ effect sizes, bootstrap 95% CIs, post-hoc achieved power, and TOST equivalence testing. The protocol in Section 8 should serve as a template.

GNN surrogates provide point estimates of energy; well-designed surrogate-assisted schemes must propagate the surrogate's uncertainty into the search. Bayesian neural networks, ensembles, and Gaussian-process surrogates are competing approaches. The 2026 literature has not yet converged on a standard protocol.

Current NISQ implementations are limited to 8–20-qubit lattice fragments. Scaling to 100+ qubits requires either fault-tolerant hardware or hybrid quantum-classical schemes that partition the problem. The 2026 IBM Heron and Google Sycamore hardware roadmaps anticipate 1000-qubit devices by 2028, which would enable meaningful protein-folding subroutines.

UCB-bandit adaptive operator selection is a 2024–2026 trend with strong empirical results. Open questions include: 1) the optimal reward window length, 2) cold-start strategy for early generations, 3) handling non-stationary reward distributions, and 4) integrating with surrogate-assisted schemes.

The community must adopt containerised execution, versioned data, public random seeds, and continuous integration for benchmarks. The CASP infrastructure provides a partial template; the OpenComplex and IDPBench benchmarks extend it to post-AlphaFold tasks. We call for a community-curated benchmark suite analogous to BBOB/COCO in the optimisation community.