To build this TRPM8 virtual screening (VS) benchmark, we compiled known TRPM8 modulators from three publicly available databases: PubChem, ChEMBL, and BindingDB (access date: June 2024; Table 1). All retrieved molecules were processed using a previously established data-curation workflow (Tran-Nguyen et al., 2023). Reported activity measurements were converted to their corresponding negative logarithmic values (pAct; Equation 1), and the median pAct value for each compound was used as its final activity metric. Compounds with pAct ≥ 6 (corresponding to activity ≤1 μM) were classified as actives. Only compounds classified as actives were retained for benchmarking. We assume that the reported TRPM8 activity values reflect ligand interactions at the well-characterized ligand-binding pocket within the voltage-sensing–like domain (VSLD), which mediates the interaction of most known TRPM8 modulators (Xu et al., 2020). This assumption is supported by available cryo-electron microscopy structures complexed with either the agonist icilin or the antagonist AMTB. p A c t = − log 10 IC 50 or  K d  or  EC 50 (1)

Agonist and antagonist annotations were assigned using a combination of functional assay descriptions and activity type (IC₅₀, K_d, and EC₅₀) information, depending on data availability. For compounds retrieved from ChEMBL and PubChem, functional assay descriptions were manually inspected to ensure that reported EC₅₀ and IC₅₀ values corresponded to receptor activation or inhibition, respectively. Compounds with ambiguous assay descriptions or unclear activity types (e.g., Kb or unspecified activity annotations) were excluded. For compounds retrieved from BindingDB, where functional assay descriptions were not available, classification was based on the reported activity types and their potencies. EC₅₀ values were used to identify agonists, while IC₅₀ values were used to identify antagonists. For compounds with multiple reported activity types, pEC₅₀ and pIC₅₀ values were compared, and compounds were classified according to the activity type associated with the stronger potency. For compounds with both pK_d and pEC₅₀ values, compounds with pK_d exceeding pEC₅₀ were classified as antagonists. The curated dataset was dominated by IC₅₀ measurements (∼2,300 compounds), a smaller number of EC₅₀ measurements (∼220 compounds) and two K_d values, with only minimal overlap between assay types.

We gathered and curated a total of 1,984 TRPM8-active compounds. To generate a challenging class-imbalanced dataset, we aimed at generating 50 property-matched decoys per active using DeepCoy (Imrie et al., 2021), a graph-based deep-learning approach that designs tailored decoys with similar physicochemical profiles but distinct topological structures. The resulting chemical space and property distributions are shown in Supplementary Figures S1A,B. As usual, not all the requested decoys were successfully generated. This resulted in an average of 48 decoys per active, leading to a final benchmark dataset comprising 1,984 TRPM8 actives and 94,650 decoys.

A total of 26 TRPM8 structures were available in the PDB (June 2024), all solved by cryo–electron microscopy. Of these, only one corresponds to the human protein (apo), while the remaining 25 are of mice (18) and avian (7) origin. For structure-based docking, we focused exclusively on ligand-bound target structures. Because mouse TRPM8 provides multiple ligand-bound structures and shares ∼94% sequence identity with the human TRPM8 channel (Yin et al., 2022), these structures were selected for docking studies. Among these, only those with the highest resolution (<3 Å) were retained.

Given that TRPM8 mediates both agonist and antagonist activities, we incorporated structural diversity in our docking experiments by using two distinct TRPM8 conformations: an agonist-bound structure TRPM8^7WRE (PDB ID: 7WRE; icilin-bound; 2.5Å) and an antagonist-bound structure TRPM8^9B6G (PDB ID: 9B6G; AMTB-bound; 2.81Å). This dual-structure strategy enables us to capture the conformational variability of the TRPM8 binding site and to evaluate which structural state provides the most reliable basis for VS.

We prepared the selected TRPM8 structures using the DockPrep module in UCSF Chimera v1.17.3. All non-standard residues were removed from the TRPM8 structures. The TRPM8^7WRE structure contains a Ca²⁺ ion located near the ligand-binding site within the calcium-binding pocket (Zhao et al., 2022). This ion was retained to preserve the agonist-bound conformation. Further, solvent molecules were removed, hydrogens were added and optimized for hydrogen bonding. Partial charges were assigned using the Gasteiger method, and any missing side chains were rebuilt using the Dunbrack 2010 rotamer library (Shapovalov and Dunbrack, 2011).

Each of the 96,634 molecules was prepared using the RDKit cheminformatics toolkit (version 2024.09.1). Salts were removed using the SaltRemover module, and hydrogens were added to the resulting chemical structures. One single 3D conformer per molecule was then generated using the ETKDGv3 method (Riniker and Landrum, 2015) employing the EmbedMolecule function. A fixed random seed (0xf00d) was applied to ensure reproducibility and the generated conformers were written in SD format.

rDock (Ruiz-Carmona et al., 2014) is an open-source docking program originally developed for VS against RNA targets but later extended to support protein–ligand docking. Its sampling protocol combines a genetic algorithm with Monte Carlo perturbations followed by local optimization. rDock employs an empirical scoring function (SF) with molecular-mechanics terms estimating the van der Waals, electrostatic, hydrogen-bonding and hydrophobic contributions to binding, with an optional desolvation term. The default scoring function (SF3) excludes the desolvation term and incorporates a repulsive polar penalty and works better with protein targets (Ruiz-Carmona et al., 2014). Accordingly, all docking scores in this work were computed using the default SF3 scoring function.

To analyze the molecules with rDock, we generated first an active-site parameter file specifying the binding site, search space, receptor flexibility, and SF. This file was processed using the rbcavity tool to build the receptor grid around the binding pocket. rDock also allows limited receptor flexibility during parameter-file generation, enabling rotation of terminal–OH and–NH₃ ⁺ groups on residues within a user-defined distance from the binding site.

To determine the optimal docking settings, redocking was performed for each protein structure, and the parameter configuration yielding the lowest RMSD between the bound and docked ligand were selected. Several parameters and combinations for docking box (6 Å, 10 Å and 15 Å) and binding site flexibility (3 Å, 5 Å and 6 Å) were tested. In the final parameter files, residues within 5 Å of the binding site were treated as flexible for TRPM8^7WRE, whereas a 3 Å flexibility radius was used for TRPM8^9B6G. For both structures, the docking cavity was defined as a 10 Å region surrounding the bound ligand.

We also dock the same molecules using Smina (Koes et al., 2013), a fork of AutoDock Vina (Trott and Olson, 2010). Smina employs a SF derived from AutoDock Vina, which is an empirical scoring model incorporating Gaussian steric terms, hydrophobic interactions, hydrogen bonding, and a rotatable-bond penalty. Ligand sampling is carried out using an iterated local search algorithm that combines stochastic perturbations with gradient-based local optimization to identify low-energy binding poses within the protein pocket. The degree of conformational search was controlled by the exhaustiveness parameter, where higher values perform more extensive sampling at the cost of increased computational time (Agarwal and Smith, 2023).

Smina docking parameters were optimized through redocking experiments analogous to those performed with rDock. Several values of exhaustiveness (1, 8, and 16) and search box dimensions (10 Å, 15 Å, 30 Å, and 35 Å) were evaluated (data not shown). Based on these tests, an exhaustiveness value of 8 was selected for TRPM8^7WRE whereas an exhaustiveness value of 1 was selected for TRPM8^9B6G combined with a 15 Å × 15 Å × 15 Å search box for both the protein structures, as these settings provided the best balance between pose recovery and computational efficiency. These parameters were therefore used for all subsequent large-scale VS calculations. Prior to docking the dataset molecules, they were converted to MOL2 format from SD format using OpenBabel v3.0.0. For each ligand, Smina was configured to generate and output a single best-scoring pose, selected from the full conformational search, and the predicted binding affinity of this pose was used for downstream performance evaluation.

Combining results from different docking programs to generate a consensus score or rank has been reported to improve VS performance (Chang et al., 2010; Charifson et al., 1999; Scardino et al., 2025). In this study, we used each docking tool to rank the benchmark molecules by decreasing predicted potency (i.e., lower ranks are assigned to molecules predicted to have more potent TRPM8 activity). To form the rDock-Smina consensus model, for each molecule, we averaged its two ranks, one from rDock and the other from Smina. The performance of this consensus protocol was then evaluated by examining the enrichment within the top 1% of the ranked list.

We also implemented a hierarchical protocol, in which molecules were filtered sequentially using the two docking tools. This approach is motivated by the observation that rDock is computationally faster, whereas Smina prioritizes more potent molecules (Supplementary Figure S2). First from the rDock run, we selected the top 10% of molecules using their rDock scores. We then retrieved the Smina scores of this rDock-filtered subset and further selected the top 10% using their Smina scores. This results in a subset with the top 1% of the molecules by this protocol, on which performance was evaluated. This protocol is practically as fast as rDock (the overhead is the time to docking 10% of the molecules with Smina).

The primary metric used for evaluating VS performance was the hit rate within the top 1% of the ranked molecules. Hit rate is defined as the proportion of true active compounds retrieved among the top-ranked candidates generated by the docking protocols. A higher hit rate indicates greater efficiency of a given protocol in identifying active hits (Equation 2). We also calculate this metric at the top 0.5%. Hit Rate  % = N u m b e r   o f   t r u e   a c t i v e s T o t a l   n u m b e r   o f   c o m p o u n d s   s c r e e n e d × 100 (2)

In addition to hit rate, we analyzed the predicted binding affinity scores across the full dataset to assess each protocol’s ability to distinguish true positives from false positives. The time to screen the molecules was taken as the CPU time required to complete that docking run.

The docking speeds of Smina and rDock were evaluated using a single-core, single-CPU setup. The ligand dataset was divided into six subsets according to molecular weight (MW), with each subset representing a distinct MW range. rDock consistently outperformed Smina in docking speed. As illustrated in Figure 1, rDock maintained near-constant docking times across all MW ranges, whereas Smina exhibited a strong dependence on ligand size, with docking time increasing progressively with MW.

Out of 96,634 ligands submitted for docking, rDock successfully docked 96,614 dataset molecules against both TRPM8^7WRE and TRPM8^9B6G. This corresponds to a docking success rate of 99.98% for each target, defined as the fraction of ligands for which a valid docking pose and score were generated. Whereas Smina docked 96,634 ligands against TRPM8^7WRE and 95,453 ligands against TRPM8^9B6G, corresponding to success rates of 100% and 98.78%, respectively. These results indicate good computational stability of both docking tools for large-scale VS. In terms of computational efficiency, rDock completed the full docking process in approximately 1 h, whereas Smina required approximately 4 h, making rDock about four times faster than Smina on this dataset.

rDock outperformed Smina by ranking a larger number of active molecules at the top position among their respective decoy sets (TRPM8^7WRE: rDock- 442; Smina- 317; and TRPM8^9B6G: rDock- 442; Smina- 294), reflecting superior early enrichment against both target structures (Figure 2). To further summarize overall ranking performance, we define here the target-wise weighted average rank (WAR) for both Smina and rDock using Equation 3, where a lower value indicates better overall ranking quality. WAR = Ʃ r a n k × n u m b e r   o f   a c t i v e   a t   t h a t   r a n k N u m b e r   o f   a c t i v e s (3)

For TRPM8^7WRE, Smina and rDock achieved WARs of 12.91 and 16.86, respectively, while for TRPM8^9B6G, the corresponding WARs were 11.62 and 16.61. These results indicate that, despite rDock exhibiting stronger first-rank enrichment, Smina achieves superior overall ranking quality across both targets. This is consistent with Figure 2, where Smina consistently ranks a large fraction of actives within approximately the top 2–15 positions relative to their associated decoys.

rDock demonstrated similar hit rates across both TRPM8 conformations, whereas Smina performed better against TRPM8^9B6G than against TRPM8^7WRE. At a hit rate of 1%, rDock achieved the highest performance across all protocols and structures, while Smina exhibited the lowest (Supplementary Figures S3A,B). At a hit rate of 0.5% (Figure 3), the consensus and hierarchical protocols showed the highest performance against TRPM8^9B6G (38% and 28%, respectively; Figure 3B), whereas the consensus protocol performed best against TRPM8^7WRE (30%, Figure 3A). These results indicate the benefit of using the integrated protocols to improve prioritization of true actives when selecting a small, high-priority subset of compounds for downstream experimental testing, compared to individual docking tools.

Analysis of the predicted docking score distributions revealed that rDock produced extreme range of scores, that is more extreme negative values, than Smina for TRPM8 structures (Supplementary Figures S3C,D). However, for both docking tools, the score distributions of true positives and false positives overlapped substantially, with very similar median values. This indicates that neither scoring function was able to reliably separate active molecules from decoys based solely on predicted docking scores, among the top 1% ranked compounds.

To further assess the added value of hierarchical and consensus protocols, we examined the overlap of true actives between methods at the top 0.5% of the ranked molecules for each method using UpSet plots (Figures 3C,D). rDock identified the largest number of unique or non-redundant active molecules (TRPM8^7WRE: 80, TRPM8^9B6G: 73), followed by consensus protocol (TRPM8^7WRE: 61, TRPM8^9B6G: 50) for both target structures, while hierarchical protocol showed considerable overlap with consensus protocol (TRPM8^7WRE: 41, TRPM8^9B6G: 82). Although Smina identified a small number of unique actives (TRPM8^7WRE: 14, TRPM8^9B6G: 49), when considering rDock and Smina alone, it added few actives beyond those already supported by either consensus or hierarchical protocols. This is reflected by the absence of an R–S only intersection in the UpSet plots. In contrast, both consensus and hierarchical screening recovered actives that were not identified by rDock and Smina alone. While consensus protocol therefore provides greater non-redundancy, it requires docking of the full compound library. We also provide UpSet plots for the less stringent 1% cutoff (Supplementary Figure S4).

From a practical screening perspective, these performances must be weighed against computational cost. Selecting the top 0.5% or 1% of compounds using rDock alone would require experimental testing of 483 or 966 molecules, respectively. However, hierarchical screening selects the same number of compounds as rDock at 0.5% and 1%, but includes additional true actives that are not prioritized by rDock alone. Importantly, this gain is achieved with limited additional computational expense, as the hierarchical protocol requires Smina docking of only 10% of the library (9,663 compounds), compared to full library docking for consensus protocol, thus providing an efficient compromise between enrichment performance and runtime.

Finally, Morgan fingerprint (ECFP4)-based pairwise Tanimoto similarity analysis between the true actives identified by Smina and rDock (Figures 3E,F) shows predominantly low-to-moderate similarity. This indicates that the two docking programs tend to retrieve chemically distinct compounds, consistent with limited scaffold overlap and the exploration of different regions of chemical space, as reported in previous studies (Guo et al., 2024). This suggests that compounds prioritized by the hierarchical protocol, which are supported by agreement between both docking tools, are likely to represent higher-confidence candidates for experimental validation.

Despite the use of both agonist-bound and the antagonist-bound TRPM8 structure, neither docking protocol retrieved any known agonists. Using TRPM8^7WRE, Smina retrieved 36 true actives in the top-ranked subset, all of which were antagonists. Similarly, Smina docking against TRPM8^9B6G retrieved 104 true actives, again consisting exclusively of antagonists.

Consistent results were obtained with rDock, which retrieved 266 and 258 true actives against TRPM8^7WRE and TRPM8^9B6G, respectively, all of which were antagonists. Notably, none of the 135 compounds classified as agonists in the dataset were recovered among the top-ranked compounds by either docking protocol or either receptor structure.