1 Introduction↩︎

Generative models are widely used in molecular design for drug discovery, enabling the generation of large numbers of candidate compounds through efficient exploration and exploitation of chemical space [1]. However, the ability to generate chemically valid molecules alone is insufficient for practical drug development [2], [3]. To serve as meaningful starting points for medicinal chemistry, generated molecules must satisfy multiple criteria related to drug-likeness, developability, and biological relevance [4].

Standard evaluations of molecular generators emphasize intrinsic metrics such as validity, novelty, uniqueness, and distributional similarity [5]. These metrics are useful for assessing whether a model produces chemically plausible and nontrivial outputs; however, they are not sufficient on their own to establish practical value in early-stage drug discovery. For those applications, compounds are prioritized using multiple medicinal chemistry and drug discovery criteria beyond potency alone, including physicochemical properties [6], structural liabilities [7], synthetic tractability [3], [8], [9], and structure-based screening procedures [10], [11].

This gap between intrinsic generation quality and practical utility is especially important in hit identification. Early discovery campaigns do not operate on raw strings or abstract graph distributions. Instead, they use molecules that balance multiple properties [12]. Later analysis stages are more time- and resource-intensive and often constrained by limited downstream compute budgets. Therefore, these stages are applied selectively to a prescreened set of molecules rather than exhaustively on the full generated set [13], [14]. Molecules with implausible properties, reactive or unstable motifs, poor synthetic tractability, or incompatible binding poses consume resources without improving the quality of the process. As a result, evaluating generators only with broad distributional or structural metrics can substantially overestimate their usefulness for medicinal chemistry.

Several task-specific benchmarks have incorporated docking scores, property objectives, or synthetic accessibility proxies. However, these objectives are often evaluated individually or in simplified combinations, and they rarely reproduce the sequential attrition imposed by practical hit identification workflows. As a result, a model may appear successful by optimizing a property, pharmacophore, or docking proxy while still producing molecules that fail basic medicinal chemistry, retrosynthesis, or pose plausibility checks.

To address these limitations, we develop HEDGEHOG, a multi-stage benchmark designed to evaluate generated molecules under a computational triage workflow that more closely reflects hit identification practice. The benchmark applies a fixed sequence of filters covering preprocessing, descriptor screening, structural alerts, synthesis feasibility, docking and affinity estimation, and post-docking three-dimensional (3D) checks. Rather than relying on common generative metrics, HEDGEHOG measures whether model outputs remain plausible after running medicinal-chemistry-oriented criteria.

This work provides a unified protocol that can be applied to unconditional, ligand-based, and protein-based models, enabling stage-wise comparison under the same downstream constraints. We instantiate HEDGEHOG on molecules targeting the switch-II pocket of KRAS G12D, a therapeutically important and structurally challenging target [15], [16]. In addition to cumulative survival, we report conditional survival at each stage in order to identify where a model class loses viability as the screening procedure becomes more stringent. This framing makes it possible to separate strong performance on conventional intrinsic metrics from performance under a more demanding computational triage setting.

These results establish HEDGEHOG as a stress test for molecular generators, exposing failure modes that are not captured by conventional intrinsic metrics. The proposed benchmark reframes molecular generation evaluation as a problem of useful chemical space proposal rather than unconstrained molecule generation. By identifying where generated molecules fail during a realistic triage workflow, HEDGEHOG provides a more informative basis for developing molecular generative models that are better aligned with early drug discovery.

2 Related Work↩︎

Distribution-learning benchmarks. Existing benchmark suites for molecular generation largely assess two-dimensional (2D) molecular quality and how closely generated sets reproduce reference distributions. GuacaMol [5] provides standardized tasks for both distribution learning and goal-directed optimization, reporting metrics such as validity, uniqueness, novelty, KL divergence, and Fréchet ChemNet Distance [17]. It incorporates compound quality metrics using rule-based filters derived from SureChEMBL, Glaxo [18], PAINS [19], and in-house rule sets. MOSES [20] standardizes dataset splits and intrinsic distributional metrics. It reports the fraction of molecules passing medicinal chemistry filters such as MCF and PAINS, and evaluates SA score, although this score is a heuristic proxy for synthesis feasibility. These benchmarks do not include comprehensive synthetic feasibility estimation, docking and binding pose assessment, or 3D structure-dependent tests.

Structure-based benchmarks. Recent benchmarks move closer to practical structure-based evaluation. DOCKSTRING [21] packages docking into a reproducible benchmark with a large target panel and tasks such as virtual screening and de novo docking score optimization. DOCKSTRING also reports a small set of molecular descriptors, including HBDs, HBAs, rotatable bond count, and logP. GenBench3D [22] evaluates 3D molecular generators inside binding pockets, introduces geometry-aware Validity3D criteria for conformational plausibility, and reports docking-based scores from Vina [23], Glide [24], and Gold PLP [25]. It computes a narrow panel of molecular descriptors, including molecular weight, logP, and ring proportion, and uses SA score as a proxy for synthesizability. Durian [26] further expands structure-based 3D evaluation by incorporating protein–ligand complexes with experimental affinity annotations. It uses a broader panel of physicochemical and geometric metrics, and evaluates docking-based affinity with QuickVina2 [27], Surflex [28], and GNINA [29]. However, these benchmarks do not natively integrate synthesizability estimation, structural alert filtering, or a representative descriptor panel.

Synthesis-aware benchmarks. Synthesis-aware evaluation has also received increasing attention. SDDBench [30] argues that standard synthetic-accessibility scores are insufficient and instead evaluates molecules using a retrosynthesis-based criterion together with search success rate. TARTARUS [31] introduces practical objectives based on physical simulation and explicit computational budgets. Its drug design tasks combine docking objectives with structural constraints and medicinal chemistry filters, while also incorporating metrics such as QED, TPSA, and SA score. These works move evaluation closer to downstream utility, but they are not designed for efficient reduction to actionable chemical space.

General-purpose evaluation framework. Frameworks such as MolScore [32] and TDC [33] provide flexible infrastructure for composing objectives, datasets, and evaluation modules. They are valuable for building and testing molecular design workflows, and in some cases contain many of the components that a practical pipeline would use. Our aim is complementary. Rather than proposing another general-purpose scoring framework, we define a highly configurable benchmarking protocol with a fixed staged structure and stage-wise reporting, motivated by hit identification triage. HEDGEHOG evaluates whether generated molecules continue to survive under a sequence of increasingly restrictive and practically motivated filters. Table 1 summarizes the relationship between prior work and HEDGEHOG.

4pt

HEDGEHOG benchmark protocol. HEDGEHOG evaluates model outputs by whether compounds pass a staged filtration cascade that reflects how medicinal chemistry teams triage molecules before hit identification. It analyzes how different generator classes fail under fixed downstream constraints. Beyond standard intrinsic metrics, HEDGEHOG computes \(21\) physicochemical descriptors, and applies structural alerts and graph sanity checks using approximately \(2,460\) SMARTS patterns and \(10\) rule sets. It assesses synthesizability using four criteria together with explicit synthesis route finding in AiZynthFinder [34], and evaluates docking using smina [35], GNINA [29], and Matcha [36] alongside binding affinity prediction using Boltz-2 [37]. Finally, our benchmark applies 3D molecular filters, including pose quality, conformer deviation, similarity, and interaction-based assessments. Appendix Table 3 contains software versions for each tool.

3 Methods↩︎

HEDGEHOG (Fig. 1) comprises six stages: molecular preprocessing, physicochemical filtering, structural filtering, synthesis feasibility assessment, docking and binding affinity estimation, and 3D post-docking analysis. The staged design is motivated not only by medicinal chemistry practice but also by computational efficiency: retrosynthetic route finding, docking, and binding affinity estimation are substantially more expensive than early descriptor- and structure-based checks. HEDGEHOG therefore uses a coarse-to-fine filtration cascade that removes implausible molecules before sending a much smaller candidate set to the slowest downstream stages. The overall thresholds and rule sets are target-agnostic. In this study, the benchmark is instantiated for KRAS G12D at the switch-II pocket (PDB ID: pdb_00007ew9 [38]). Full configuration details per stage and thresholds are provided in Appendix Section 6.

3.1 Stage 1: Molecular preprocessing↩︎

All generated molecules are first passed through a common preprocessing procedure implemented with RDKit [39] and DataMol [40]. This stage converts generated strings into molecular objects, removes salts and solvents, retains the largest relevant organic fragment, disconnects metals, standardizes representation through normalization and sanitization, and converts back to SMILES strings. Molecules containing atom types outside the supported element set, radicals, isotopes, invalid valence states, or unresolved multi-fragment artifacts are removed at this stage. This preprocessing serves two purposes. First, it standardizes outputs from heterogeneous generators before evaluation. Second, it reduces the risk that downstream calculations are affected by malformed or chemically inconsistent inputs (Appendix Section 6.1).

3.2 Stage 2: Physicochemical descriptors↩︎

After preprocessing, HEDGEHOG computes a panel of 21 physicochemical descriptors and filters molecules using configurable threshold ranges (Appendix Table 4). The panel includes size- and composition-related descriptors, lipophilicity and polarity descriptors, hydrogen bond features, ring statistics, fraction of sp³-hybridized carbon atoms, and additional medicinal chemistry checks. We do not adopt a single canonical rule set like Lipinski [41] or Veber [6], because none of them covers the full descriptor panel used in our benchmark, and published thresholds vary across target classes, screening settings, and medicinal chemistry objectives. Instead, we use broad descriptor bounds that combine literature-reported thresholds with manually curated checks for descriptors where no widely accepted cutoff exists. These thresholds are chosen to remove extreme outliers and clearly implausible molecules while preserving broad chemical diversity.

3.3 Stage 3: Structural filters↩︎

Stage 3 removes molecules with undesirable substructures or poor medicinal chemistry profiles. This stage combines public structural alert libraries with additional rule-based graph checks to remove chemotypes that are excluded in small molecule hit identification workflows, including unstable, highly reactive, toxic, and assay-interfering motifs, as well as structures with inappropriate graph-level features. The alert layer covers medicinal chemistry and assay-interference substructures such as Dundee [42], BMS [43], Glaxo [18], PAINS [19], Lilly Medchem Rules [7], NIBR [44], and related alert sets (Appendix Table 5). The rule-based layer includes graph-level checks for ring infractions, protecting groups, halogen- and stereochemistry-related constraints, and other structural features that are typically excluded during medicinal chemistry review. Because some alert categories can be context-dependent in broader drug discovery settings, the structural filters stage is fully configurable, both at the level of rule set selection and at the level of individual SMARTS patterns (Supplementary File 2).

3.4 Stage 4: Synthesis feasibility↩︎

We evaluate whether molecules are likely to be synthetically tractable by combining heuristic synthesis-related scores with explicit retrosynthesis search. In the default configuration, HEDGEHOG computes SA score [45], RA score [46], and SYBA score [47], and then runs AiZynthFinder [34] to test whether at least one synthetic route can be found under the chosen retrosynthesis setup. Because retrosynthetic search is computationally expensive, it is applied only after faster physicochemical and structural filters have substantially reduced the candidate pool. The default thresholds for synthesis scores are SA \(\leq 4.5\), SYBA \(\geq 0\), and RA \(\geq 0.5\). We use this combination because heuristic scores alone are insufficient. They are fast and informative; however, they do not provide an explicit route. At the same time, route finding success depends on the template library, stock set, search policy, and compute budget. An example of a synthetic route provided by AiZynthFinder is shown in Appendix Fig. 5.

3.5 Stage 5: Docking and binding affinity estimation↩︎

Molecules that pass the earlier stages are evaluated in a target-aware manner by docking them into the KRAS G12D switch-II pocket. HEDGEHOG supports several docking engines, including smina [35], GNINA [29], and Matcha [36], a recent neural docking method based on multi-stage Riemannian flow matching. These tools were selected based on their strong performance in independent evaluations on diverse docking datasets [48]. In the default configuration, HEDGEHOG uses all three docking tools together and estimates binding affinity with Boltz-2 [37]. Docking and affinity estimation are among the slowest stages in the pipeline, so they are performed only on molecules that survive the earlier low-cost filters. Because docking scores and affinity estimates are imperfect proxies for true binding, we extend this stage with multiple models, rather than relying on any single proxy. Survival through Stage 5 is a convergent computational support for favorable predicted activity and plausible binding.

Molecules that reach docking have already passed validity, descriptor, structural, and synthesis feasibility checks. These early stages are essential because docking scores alone can be exploited by unrealistic molecules (Appendix Fig. 6). A detailed explanation and the corresponding thresholds are provided in Appendix Section 6.5.

3.6 Stage 6: Three-dimensional filtration↩︎

The final stage evaluates docked poses using a sequence of post-docking 3D checks. These checks are intended to remove molecules whose nominal docking success is undermined by poor poses or implausible geometry. That is why we include search-box containment, pose quality assessment, protein-ligand interaction checks, and a conformer-deviation analysis. For the KRAS G12D benchmark, we impose a target-specific interaction requirement involving Asp12. The interaction is illustrated in Appendix Fig. 7. HEDGEHOG separates target-agnostic 3D checks from target-specific interaction rules so that future benchmark instances can define different target-aware criteria without changing the broader framework. The interaction stage can enforce required or forbidden residue contacts through ProLIF-based analysis [49].

4 Results and discussion↩︎

We evaluate three classes of molecular generators, including unconditional, ligand-based, and protein-based models. For each model, we sample \(N_{gen} = 10,000\) distinct SMILES strings. Invalid strings were not included in the \(10,000\) SMILES input set. During this step, uniqueness was enforced only by exact string matching. Therefore, different SMILES strings encoding the same molecular graph could still be retained. However, during Stage 1 (molecular preprocessing), redundant representations of generated molecules may be reduced through SMILES sanitization and standardization of structures, applying structure-level deduplication, and removal of chemically inconsistent molecules. For model \(m\) and stage \(s\), survival with respect to initial number of molecules was defined as \(N_{m,s}/N_{m,0}\), where \(N_{m,0}=10{,}000\). Survival with respect to the previous stage was defined as \(N_{m,s}/N_{m,s-1}\). For model level summaries, we report the mean and standard deviation of the number of retained molecules across individual generators within each model class. All experiments were conducted using an NVIDIA A100 40 GB GPU. Runtime details are provided in Appendix Section 6.7.

4.1 HEDGEHOG benchmark design and evaluated generators↩︎

HEDGEHOG is a hierarchical filtration benchmark designed to evaluate molecular generators under a staged computational workflow. Generated molecules are first evaluated with inexpensive, broadly applicable medicinal chemistry filters and only then passed to more expensive target-aware calculations. The complete workflow is shown in Fig. 1, and full implementation details and thresholds are provided in Methods and Appendix Section 6.

HEDGEHOG applies six sequential filters to the generated molecules. First, SMILES strings are standardized and chemically invalid structures, unsupported atoms, salts, solvents, and disconnected fragments are removed. Second, molecules are screened by 21 physicochemical descriptors that capture composition, lipophilicity, polarity, and molecular complexity [50]. Third, undesirable structural motifs are removed using public alert sets such as PAINS [19], Lilly Medchem Rules [7] and NIBR [44]. Fourth, synthesis feasibility is assessed with SA [45], RA [46] and SYBA [47] scores and explicit route search in AiZynthFinder [34]. Fifth, retained molecules are docked into the target protein with smina [35], GNINA [29] and Matcha [36], and evaluated by Boltz-2 [37] affinity prediction. Finally, docked poses are checked for 3D plausibility, including pocket containment, pose quality, conformer deviation and target-specific interaction.

We evaluated \(23\) molecular generators grouped into three classes of unconditional, ligand-based, and protein-based models. For each model, we generated \(10,000\) valid and distinct SMILES, excluding unparsable strings, and with uniqueness enforced at the raw string level. This yielded \(80,000\) molecules from unconditional models, \(70,000\) molecules from ligand-based models, and \(80,000\) molecules from protein-based models, for a total of \(230,000\) generated molecules. Evaluation details are provided in Methods section 3.

4.2 Survival under the HEDGEHOG workflow↩︎

We first quantified how many generated molecules survive the complete filtration cascade and at which stages attrition occurs. Table ¿tbl:tab:class-pass-rates? summarizes outputs at each stage for each model class, and reports both survival relative to the initial generated set and survival relative to the previous stage. This distinction is important because the lowest final survival is not always caused by the final target-aware filters. The per model statistics provide a complementary view of these class-level trends by reporting the mean and standard deviation of the number of retained molecules across individual generators within each model class. Since each generator contributes \(10,000\) valid and distinct SMILES strings, these values quantify both the typical survival of a single generator and the dispersion of performance within the class. Fig. 2 reports cumulative survivors after each stage.

Early-stage filtering. Preprocessing already shows differences in raw generation quality. Conditional models produce molecules that pass standardization at much higher rates than unconditional models, with average preprocessing pass rates \(97.56\%\) and \(75.51\%\), respectively. At preprocessing, unconditional models retain \(7551 \pm 3476\) molecules per generator, whereas ligand-based and protein-based models retain \(9837 \pm 196\) and \(9675 \pm 743\), respectively. Thus, conditional models improve early standardization not only in aggregate, but also with greater consistency across generators. As shown in Fig. [fig:un95survival95plot], several unconditional generators, such as E(3)DM and TGM-DLM, lose a substantial fraction of molecules before any chemical filtering is applied, highlighting poor raw generation quality.

Descriptor filtering is the first major bottleneck in the benchmark. High retention rates on individual filters do not guarantee simultaneous satisfaction, which is substantially more difficult. Therefore, the main challenge is multi-property balance rather than isolated property failure. According to Table ¿tbl:tab:class-pass-rates?, unconditional, ligand-based, and protein-based model classes are reduced to smaller sets of \(24.93\%\), \(28.54\%\), and \(24.27\%\) respectively. This convergence at the class level masks substantial variation between individual generators. The corresponding per model values are \(2493 \pm 2062\), \(2854 \pm 2109\), and \(2427 \pm 1423\) molecules for unconditional, ligand-based, and protein-based models, respectively. The large standard deviations indicate that descriptor compliance is highly model-dependent, even within the same generation setting. Ligand-based models perform slightly better, but the variation is large, for example, GCPG and GENTRL retain 5,849 and 5,451 molecules respectively, while PGMG drops to 97 molecules. Structural filters remove an additional fraction of molecules. Unconditional, ligand-based, and protein-based model classes fall to \(5.82\%\), \(5.90\%\), and \(3.62\%\) of their initial molecules, respectively. This suggests that many molecules that satisfy a broad panel of physicochemical descriptors still contain prohibited substructures. Without a structural evaluation stage, models may appear competitive in later non-structural filtering stages while systematically enriching undesirable chemotypes. Since the descriptor and structural filtering stages are target-agnostic and rely on simple, broadly applicable criteria, the fact that only \(5.11\%\) of all initial molecules remain for synthesis evaluation suggests that many generated compounds fail general standards of chemical plausibility and developability. This implies that the main weakness of most models is not only in later proxy-based optimization, but already in their ability to generate molecules that satisfy common baseline requirements for drug-like chemical space. The staged design of HEDGEHOG also serves a practical purpose as a coarse-to-fine computational cascade. By removing most candidates before expensive downstream analyses, these early stages reduce the number of molecules entering synthesis feasibility and docking evaluation, thereby lowering overall computational time and cost.

Synthesis evaluation and docking-based filtering. Explicit retrosynthetic planning is substantially more selective than heuristic synthesizability scores alone. After synthesis feasibility filtering, unconditional, ligand-based, and protein-based model classes retained \(3.47\%\), \(2.12\%\), and \(1.65\%\) of their initial molecules, respectively. This shows that satisfying score-based proxies alone can overestimate practical synthetic accessibility. Conditional pass rates show that synthesis filtering is not equally severe for all classes. Among molecules that passed structural filtering, unconditional, protein-based, and ligand-based models retained \(59.72\%\), \(45.44\%\), and \(35.89\%\), respectively. Thus, ligand-based or protein-based generation did not necessarily improve robustness under broader synthetic tractability constraints. The remaining molecules were then evaluated by docking and binding affinity filters. Relative to the initial set, unconditional models retained the highest fraction of molecules, compared with ligand-based and protein-based models, retaining \(1.80\%\), \(1.55\%\), and \(0.96\%\), respectively. As shown in Supplementary File 1, docking with a single tool appears to result in tool-specific false positive pass rates. The bottleneck at this stage is not strong performance in any single docking or affinity model, but agreement across models with different inductive biases. Molecules that satisfy score-based docking criteria may still fail geometric or conformational plausibility checks. The final 3D filters assess whether each retained docked pose is geometrically plausible, conformationally accessible, and consistent with the target-specific interaction hypothesis.

Model-class performance overview. Unconditional models retain the highest fraction of molecules through the docking and binding-affinity stage relative to their initial set (\(1.80\%\)), compared with ligand-based models (\(1.55\%\)) and protein-based models (\(0.96\%\)). This result should be interpreted after the preceding target-agnostic filters. Under HEDGEHOG, where docking is evaluated only after these chemically implausible molecules are removed, explicit target conditioning does not lead to higher survival at the docking and binding affinity estimation stage.

Ligand-based models show the highest conditional survival from the synthesis feasibility stage to the docking and binding affinity stage, retaining \(73.10\%\) of synthesis stage survival molecules. However, unconditional models retain the largest absolute number of molecules at this stage because they enter it with a larger post-synthesis pool. Therefore, the lower end-to-end survival of ligand-based models (Fig. [fig:lb95survival95plot]) is the combined effect of earlier losses.

Among model classes, protein-based generators have the lowest cumulative survival through docking and affinity estimation (Fig. [fig:pb95survival95plot]), yet the highest conditional survival at the 3D stage, retaining \(63.15\%\) of docking survivors. This suggests that, once protein-conditioned molecules pass the earlier filters and docking affinity criteria, their retained poses are more likely to satisfy the final geometric and interaction-based constraints. However, because fewer protein-based molecules survive the earlier target-agnostic stages, this local advantage does not translate into the highest end-to-end survival.

The last stage reduces the benchmark to a small set of molecules that satisfy all HEDGEHOG filters. Unconditional, ligand-based, and protein-based model classes retain \(609\), \(396\), and \(485\) molecules in total, corresponding to \(0.76\%\), \(0.57\%\), and \(0.61\%\), respectively, suggesting that target conditioning alone does not guarantee higher end-to-end survival. Protein-based models retain \(61 \pm 110\), where the standard deviation is substantially larger than the mean. This means that the protein-based class total is driven disproportionately by the strongest Dragonfly model rather than by uniformly strong performance across the class.

HEDGEHOG selects for a narrow combination of valid generation, medicinal chemistry quality, synthetic accessibility, binding tool agreement, and physically plausible 3D poses. This explains why models with high early-stage retention do not necessarily produce the best final candidates. The differences between ligand-based REINVENT4 configurations show that the broad out-of-the-box prior outperforms both the similarity-biased and transfer-learned versions, suggesting that stronger target bias can overconstrain generation and reduce robustness under broad medicinal chemistry filters.

To provide a compact cross-class summary, Table 2 reports the best performing generators within each model class, ranked by the number of molecules that pass the full HEDGEHOG pipeline. Dragonfly is the best performing model in this KRAS G12D benchmark instance, retaining nearly twice as many final molecules as the second-best REINVENT4 (V) model with 182 molecules. The result is consistent with the design of Dragonfly, which combines interactome-based learning with molecular generative priors and can incorporate ligand templates, 3D binding site information, and physicochemical property constraints. The result suggests that protein conditioning is most effective when it is combined with strong chemical priors that preserve scaffold coherence, structural plausibility, and drug-like molecular features. The comparison between Dragonfly and Dragonfly (b) illustrates this point. Dragonfly (b), which uses target-ligand physicochemical descriptor generation, passes a descriptor stage similar to the out-of-the-box Dragonfly model, retaining \(9,990\) molecules and \(9,967\), respectively. However, it falls to \(87\) molecules after structural alerts evaluation, while the out-of-the-box Dragonfly retains \(2,488\) molecules. This indicates that matching ligand descriptor profiles can improve early physicochemical compatibility, but does not by itself prevent structural liabilities.

Across \(230,000\) initial molecules, only \(1,490\) molecules survive all stages. The stage-wise results show that benchmark success is highly local, because model rankings vary across different filtering stages. The diagnostic value of HEDGEHOG is not limited to ranking unconditional, ligand-based, and protein-based models. The stage-wise failures also separate architectural and training objective effects. Scaffold- and motif-aware graph models tend to preserve chemically coherent structures through early filters, ligand-based models can improve target-specific relevance while losing robustness under synthesis and structural alert constraints, and protein-based 3D generators can show stronger pose-level survival only after severe target-agnostic attrition. Thus, the main failure mode is not a single missing property but the inability to maintain chemical realism, synthesizability, and target compatibility simultaneously. Results show that HEDGEHOG is a stringent and highly discriminative benchmark that reduces the candidate pool under a realistic multi-stage computational triage workflow and clearly separates performance across model classes.

5 Conclusions↩︎

Molecular generator evaluation is misleading when based on a limited set of scores and metrics. Our results suggest that generator evaluation should report end-to-end survival across medicinal chemistry, structure-based, and synthesis filters. We show that evaluated generators are still poorly aligned with stringent end-to-end computational triage. Under HEDGEHOG, the large majority of generated molecules fail before the proxy stages, and success at one isolated stage does not guarantee survival through the full pipeline. Under this KRAS G12D benchmark instance, Dragonfly has the highest end-to-end survival with \(345\) final molecules, followed by REINVENT4 (V) with \(182\), and REINVENT4 with \(163\) final molecules. The benchmark therefore exposes a central limitation of the current paradigm, that many models can generate molecules that look promising under partial or local criteria, yet very few produce outputs that remain acceptable under a broader sequence of medicinal-chemistry-oriented constraints.

A limitation of our evaluation is that we do not estimate variance across independent generation seeds for all models. Some baseline implementations did not expose or consistently document seed control, and for fairness we therefore use a common evaluation protocol based on \(10,000\) generated SMILES strings per model. This benchmark instance is centered on the KRAS G12D switch-II pocket and includes a target-specific Asp12 interaction requirement. Model rankings may therefore change for other targets, pocket classes, or interaction hypotheses.

While our pipeline integrates physicochemical, structural, synthesis, docking, and post-docking 3D evaluation stages, it does not yet explicitly model long-term ADMET properties, pharmacokinetics, or clinical-success predictors. Future work should extend the benchmark with uncertainty-aware generative modeling, multi-objective optimization under experimentally grounded constraints, broader target coverage, and prospective validation. By shifting focus from isolated benchmark metrics toward survival in actionable chemical space with molecules that satisfy basic medicinal chemistry criteria, HEDGEHOG provides a stricter and more practically motivated framework for evaluating molecular generators.

Data Availability statement↩︎

The HEDGEHOG source code is available at https://github.com/LigandPro/hedgehog.

Appendix↩︎

6 Detailed molecular workflow of the HEDGEHOG pipeline↩︎

The code is available at https://github.com/LigandPro/hedgehog. Unless otherwise stated, all experiments use the core benchmark dependencies listed in Table 3.

Supplementary File 1 reports the number of generated molecules that satisfy the corresponding threshold for each HEDGEHOG filtration stage, split by model class: unconditional, ligand-based, and protein-based. The bold rows report the number of molecules that satisfy all filter thresholds simultaneously at each stage.

At each stage, a molecule is evaluated at stage \(s+1\) only if it passed all criteria at stage \(s\). All thresholds were fixed before computing model-level pass rates and were applied identically to all generated molecules and model classes.

6.1 Molecular preprocessing definitions↩︎

Stage 1 standardizes all generated SMILES strings so that downstream calculations are performed on a consistent and comparable set of molecules. Different model architectures, training datasets and objectives, and sampling procedures can produce outputs in different formats. For that reason, HEDGEHOG first parses each SMILES with DataMol and RDKit, repairs valence issues when possible, and sanitizes the resulting molecule. Then explicit hydrogens are added to each molecule and checked for atom-type compatibility with a predefined restricted set of \(10\) elements: H, C, N, O, F, P, S, Cl, Br, I. Although some approved drugs contain metals, boron, silicon, or other elements, this restricted atom set is appropriate for a general small-molecule generation benchmark because it avoids model-dependent failures caused by uncommon chemistry.

Generated structures, especially from diffusion-based models, may contain disconnected fragments, salts, or solvent molecules. Therefore, HEDGEHOG retains only the largest molecular fragment and standardizes it by disconnecting metals, normalizing functional groups, and reionizing the molecule. Salts and solvents are removed because they also may be artifacts of crystallization, and we evaluate the potential of a generator to produce a drug-like molecule, rather than a full crystallographic assembly. Finally, stereochemical annotations are removed and the molecule is converted to a canonical non-isomeric SMILES representation. This makes deduplication deterministic because standardized molecular strings can then be compared directly regardless of artifacts. Figure [fig:details_preproc_heatmap] summarizes per filter pass rates.

6.2 Physicochemical descriptors definitions and thresholds↩︎

HEDGEHOG computes \(21\) physicochemical descriptors using RDKit, including one custom molecular complexity term, MCE-18 [50]. All descriptors are two-dimensional (2D) and do not depend on three-dimensional (3D) conformations.

No single literature rule set defines thresholds for all physicochemical descriptors used in the HEDGEHOG panel. Many medicinal chemistry filters were designed for specific screening pipelines and cannot be transferred directly to every generative model setting. We therefore assigned thresholds in three ways, summarized in Table 4. When a descriptor and cutoff matched a published filter directly, we used the reported threshold without modification. For descriptors that were related to published rules but not defined in exactly the same way, we used the rule as a guide for the corresponding HEDGEHOG descriptor. Finally, for descriptors that were not covered by broad drug-likeness filters, thresholds were chosen empirically after examining molecules.

Figure 4 shows how generated molecules populate each descriptor range, whereas Figure [fig:details_desc_heatmap] summarizes how many molecules from each model satisfy each descriptor threshold.

6.3 Structural filters definitions and thresholds↩︎

Structural filters remove molecules with structural motifs that are likely to be reactive, unstable, synthetically undesirable, or associated with well-known medicinal chemistry liabilities. Figure [fig:details_sf_heatmap] highlights pass rates of molecules during structural filters evaluation, and Table 5 summarizes the structural filter configuration used in the HEDGEHOG pipeline, resulting in the set of \(2,458\) underlying SMARTS patterns provided in Supplementary File 2.

6.4 Synthesis feasibility definitions and thresholds↩︎

We evaluate synthesis feasibility using four complementary metrics. The Synthetic Accessibility score (SA score) [45] is a fragment- and complexity-based heuristic that ranges molecules from 1, for molecules expected to be easy to synthesize, to 10, indicating molecules that are difficult to synthesize. We apply an SA score cutoff of 4.5. The Retrosynthetic Accessibility score (RA score) [46] is a machine-learned classifier that predicts whether an automated retrosynthesis planner is likely to identify a route to the target. We apply the natural decision cutoff of 0.5 for the predicted probability to find a synthesis route. SYBA [47] is a fragment-based Bernoulli Naive Bayes classifier. Positive values indicate molecules expected to be easier to synthesize, and we therefore apply a threshold of 0.

AiZynthFinder is a machine-learned tool for retrosynthesis route planning. In this work, we used a ZINC-derived stock of building blocks from the ZINC database [61], with a maximum search depth of 5, and a time limit of 300 seconds per molecule. If AiZynthFinder finds a route within these constraints, the search is considered successful. An example route proposed under this setup is shown in Figure 5. Failure to find a route does not imply that the molecule is intrinsically non-synthesizable, because the result depends on the selected stock, policies, depth, and time budget. Model-wise retention after the individual criteria and the final Stage 4 pass is summarized in Figure [fig:details_retro_heatmap].

6.5 Docking and binding affinity definitions and thresholds↩︎

Stage 5 is the target-aware structure-based evaluation for KRAS G12D at the switch-II pocket. HEDGEHOG supports multiple docking approaches, including smina, GNINA, Matcha and binding affinity estimation tool Boltz-2.

For molecular docking with smina, GNINA, and Matcha, we use a search box of side length \(10~\AA\) centered at \([-1.258, -12.520, 11.390]\), with exhaustiveness set to eight, nine output modes, and an energy range of \(3 kcal/mol\). Docking poses were retained if their predicted docking score was lower than -\(6.5 kcal/mol\). This threshold was chosen to remove clearly unfavorable poses while retaining molecules with potentially meaningful target interactions in the binding pocket. For Boltz-2, we used the continuous affinity output defined as \(\log_{10}(IC_{50}/1 \mu M)\). Lower values correspond to stronger predicted binding. We applied a cutoff of \(\log_{10}(100) = 2\), corresponding to a predicted \(IC_{50} \leq 100 \mu M\). Figure [fig:details_dock_heatmap] shows per filter independent pass rates with an overall pass rate of all criteria simultaneously.

The tools used in this stage make different assumptions about protein-ligand binding. Smina is a classical empirical docking engine, whereas GNINA augments the same general docking framework with convolutional neural network scoring. Matcha is a recent docking pipeline based on multi-stage flow matching with explicit physical validity filtering and GNINA-based final ranking. Boltz-2 is a biomolecular foundation model that jointly predicts complex structure and binding affinity. We therefore do not interpret a favorable score from any single tool as sufficient evidence of binding. Instead, we prioritize molecules that remain favorable across multiple tools, because agreement across these distinct modeling paradigms is less likely to reflect a tool-specific false positive result and provides a more conservative basis for subsequent pose-level validation.

Some models appear to optimize docking scores at the expense of chemical plausibility, resulting in favorable docking scores but chemically poor molecules. To illustrate this, we sampled three representative molecules. These molecules passed Stage 1 only and were docked into the KRAS G12D switch-II pocket under the same conditions as in the main pipeline. We report average docking score as the unweighted arithmetic mean of the best docking scores returned by smina, GNINA, and Matcha. Figure 6 shows that these molecules satisfy the docking scoring function while violating basic medicinal chemistry or synthesizability constraints. This highlights that a model may appear strong under a docking-only metric by generating high-scoring but chemically implausible candidates.

6.6 Three-dimensional filtration definitions and thresholds↩︎

At this stage, the goal is no longer to assess whether a molecule can in principle be placed in the pocket, but whether the retained docked pose is geometrically plausible, conformationally accessible, and consistent with the target-specific interaction hypothesis. Figure [fig:details_threeD_heatmap] summarizes per filter pass rates with an overall pass rate of all criteria simultaneously.

The final post-docking stage applies a set of complementary 3D pose filters. First, docked poses are required to remain inside the predefined docking box, ensuring that the final ligand geometry stays within the intended binding region. In the KRAS G12D setting, this region was defined around the crystallographic 05C ligand. The filter checks final ligand coordinates atom by atom and rejects any pose with atoms outside the box. Second, pose quality is evaluated using a PoseCheck Fast [52] geometry check to remove poses with severe intermolecular clashes, excessive volumetric overlap with the protein, large displacement from the pocket, or internal steric inconsistencies. The mean speed of \(0.05 ms\) per pose comes from keeping the test purely geometric and batchable. Heavy-atom coordinates are evaluated with vectorized PyTorch distance calculations. Third, conformer deviation is evaluated with a symmetry-corrected RMSD backend to reject poses that deviate too strongly from accessible ligand conformations. Each docked geometry was compared with a molecule-specific RDKit ETKDGv3 conformer ensemble. We used the minimum symmetry-corrected RMSD to this ensemble and rejected poses above the \(3.0~\AA\) threshold. Finally, protein–ligand interaction fingerprints computed with ProLIF are used to assess target-specific interaction requirements. For the KRAS G12D setting, interactions with Asp12 are biologically motivated because salt-bridge formation with the mutant Asp12 residue has been reported as a route to KRAS G12D selectivity [15]. Figure 7 illustrates this interaction. This requirement is not presented as a universal rule for all targets, but as part of the example studied here and is motivated by the target biology and binding site context.

6.7 Per-stage runtime↩︎

The HEDGEHOG workflow was implemented as a staged, parallelized screening cascade. Early, low-cost filters removed unsuitable molecules before computationally intensive retrosynthesis, docking, and affinity-prediction steps were applied. This design substantially reduced the number of candidates entering the slowest stages and enabled the complete benchmark to be run within practical wall-clock times. Per-stage wall-clock times are reported in Table 6.

7 Models overview↩︎

We categorize molecular generators by generative strategy and architecture because both impose distinct inductive biases, including validity and grammar errors for strings, geometry handling for 3D models, pocket alignment for pocket-based models [62], [63]. This taxonomy is summarized in Table 7.

3pt

Genetic algorithm (GA). GA is a heuristic optimizer that evolves molecules via crossover and mutation operations. MolFinder [64] applies conformational space annealing to SMILES [81] and requires no generative model pretraining.

Variational autoencoder (VAE). VAE models learn a latent distribution over chemical space with an encoder–decoder pair optimized through the evidence lower bound [82]. JT-VAE [66], HierGraphVAE [65], and MoLeR [67] operate on graphs with scaffold-aware decoders. GENTRL [68] is a VAE with a structured latent prior and reinforcement-learning fine-tuning.

Autoregressive models. String-based autoregressive models factorize sequence likelihood as \(\prod_i P(t_i \mid t_{<i})\). MolGPT [69] is a decoder-only Transformer. REINVENT4 [72] fine-tunes a SMILES prior into an agent via policy gradient with recurrent or transformer-based prior.

For structure-based design, 3D autoregressive models condition on a pocket \(P\) and learn the conditional likelihood of a molecule \(M\) as \(p_{\theta(M | P)} = \prod_{t = 1}^T p_{\theta(z_t | z_{< t}, P)}\), where each step \(z_t\) adds atoms, bonds, or coordinates. Pocket2Mol [73] and ResGen [74] are autoregressive 3D generators that build molecules inside the pocket. Dragonfly [4] is an interactome-based graph-to-sequence model that can incorporate ligand templates or 3D binding site information by encoding pocket geometry with SE(3) or E(3)-equivariant encoders.

Pharmacophore-based models use condition \(c\) as a set of 3D interaction features and geometry introducing latent \(z\) which, via \(p\left(x | c\right) = \int p_{\theta(x | c, z)} p\left(z\right) d z\), models the many-to-many relationship between pharmacophores and ligands. PGMG [71] represents a pharmacophore as a fully connected graph. GCPG [70] is a transformer encoder–decoder whose hidden state is modulated by pharmacophore embeddings.

Diffusion models. Diffusion models learn to approximate the reverse process of a predefined forward noising process [83]. E(3)DM [75] is an E(3)-equivariant model that jointly denoises atom coordinates and types. DiffSBDD [77] is an SE(3)-equivariant 3D-conditional model. TargetDiff [79] conditions the diffusion process on a protein binding site. TGM-DLM [76] is a text-guided diffusion language model over SMILES token embeddings. In our benchmark, it is used in a non-target-specific prompting setting and is therefore grouped with unconditional baselines. ProtoBind-Diff [78] is a structure-free diffusion language model that takes a protein amino-acid sequence as a condition.

Flow matching. Flow-matching models learn continuous-time velocity fields that transport a base distribution to the data distribution. DrugFlow [80] is a pocket-conditioned ligand generation model with flow-based sampling.

8 KRAS G12D benchmark instance↩︎

The KRAS G12D receptor was prepared from PDB 7EW9 (PDB ID: pdb_00007ew9 [38]). The crystallographic inhibitor TH-Z816, GDP, Mg\(^{2+}\), and crystallographic water molecules were removed before docking. Residue numbering was kept consistent with the PDB file, with Asp12 used as the target-specific interaction residue. The prepared receptor was used for smina, GNINA, Matcha, and Boltz-2 calculations and Stage 6 post-docking checks.

We used KRAS G12D as the benchmark target for HEDGEHOG. KRAS G12D is the most prevalent KRAS oncogenic mutation [15], and for many years was considered undruggable. Unlike KRAS G12C, it lacks the Cys12 required for covalent binding, but at the same time the negatively charged Asp12 can be targeted via salt-bridge interaction, as demonstrated in the development of TH-Z816 [84]. Given its challenging nature and clinical significance, we considered KRAS G12D a relevant and mechanistically informative benchmark target for the HEDGEHOG evaluation pipeline [15].

For unconditional models, KRAS G12D is used as the target in docking, affinity estimation and post-docking 3D checks. For protein-based models, it also provides the KRAS G12D structure used as target-conditioning input. For ligand-based models, we used a set of KRAS G12D inhibitors provided by [85] as conditioning data.

The initial set contains 645 molecules. We assessed validity using RDKit and identified 13 molecules with invalid structures due to valency issues. These molecules were removed, leaving 632 unique and valid ligands.

9 REINVENT4 configurations↩︎

We evaluated REINVENT4 in one de novo configuration and three ligand-conditioned configurations. REINVENT4 was used for this analysis because the same framework provides public priors, de novo SMILES generation, molecule-based generation, and transfer-learning workflows. This allows prior choice to be varied while keeping the downstream HEDGEHOG evaluation fixed, isolating how prior choice and transfer-learning within models affect downstream utility. For sampling, we used beam search with temperature set to \(1.0\) and number of SMILES per input SMILES set to \(500\). All REINVENT4 configurations were evaluated on one NVIDIA A100 40GB GPU. Prior checkpoints are listed in [86] and are open-source.

The unconditional baseline was denoted as REINVENT4 (checkpoint, reinvent.prior) in sampling mode. No KRAS ligand SMILES were supplied to this run.

REINVENT4 (V) (Vanilla, checkpoint: reinvent.prior). Prior was trained on ChEMBL data [87]. We input \(632\) KRAS G12D inhibitors to the model and obtained \(49,119\) valid SMILES strings in \(8\) hours. The REINVENT4 (P) (provided prior, checkpoint: mol2mol_medium_similarity.prior) was trained with an explicit medium Tanimoto similarity objective. Because the public Mol2Mol prior was trained on PubChem-derived molecular pairs [88], we removed KRAS seed ligands with exact PubChem matches before sampling to prevent data leakage. This filtering reduced the seed set to 583 molecules. The model outputs \(126,001\) valid SMILES strings in \(8\) hours. For REINVENT4 (TL) (transfer-learning), we fine-tuned a REINVENT4 prior (reinvent.prior) on \(632\) known KRAS G12D inhibitors with the Mol2Mol training pipeline. We defined a similarity threshold of \(0.7\), and fine-tuned for \(1,000\) epochs for \(72\) hours. For further evaluation, we selected the checkpoint with the lowest loss on the validation set.

References↩︎