Structure-based virtual screening for TRPM8 modulators

Transient receptor potential melastatin 8 (TRPM8) is an emerging therapeutic target, yet the performance of available structural models for docking and optimal docking protocols for virtual screening (VS) remains unclear. Here, we benchmarked two available TRPM8 structural conformations (agonist-bound TRPM8 7WRE and antagonist-bound TRPM8 9B6G) using docking tools, Smina and rDock, using known TRPM8 inhibitors and property-matched decoys. rDock achieved the highest hit rates and outperformed Smina in ranking true actives at first-ranks than their corresponding decoys, whereas Smina showed a target-structure dependence on docking performance but delivered superior overall ranking quality across both target structures. Both docking tools displayed considerable overlap between active and decoy score distributions, indicating only moderate discriminatory power of docking scores alone. When prioritizing a small subset of top-ranked compounds, integrated screening approaches, particularly the consensus protocol, improved the recovery of true actives, while the hierarchical protocol achieved comparable performance at a substantially lower computational cost. Collectively, this work establishes a reproducible VS benchmark for TRPM8 and supports the use of different screening protocols to improve early hit identification.

In the early 2000’s, studies aimed at understanding the molecular basis of cold sensation revealed TRPM8 as a key cold-activated ion channel. This ground-breaking discovery was awarded the 2021 Nobel Prize for Physiology and Medicine to David Julius and Ardem Patapoutian. TRPM8 is a polymodal, nonselective cation channel activated by cold temperatures and cooling agents such as menthol, and it can be modulated by both agonists and antagonists. Since its discovery, TRPM8 has become one of the most extensively studied TRP channels (Wu et al., 2024). Subsequent research has expanded its functional repertoire beyond thermo-sensation to include roles in cold-induced and neuropathic pain Weyer and Lehto, 2017; Bianchini et al., 2021) multiple cancers (Ochoa et al., 2023), inflammation (Ramachandran et al., 2013) and a range of additional pathological conditions (Liu et al., 2020).

As of 28 May 2025, the U.S. Food and Drug Administration (FDA) approved the first-in-class drug targeting TRPM8 for the treatment of dry-eye disease (Zhou et al., 2025), nearly two decades after the channel’s initial discovery. This milestone underscores the need for more cost-effective and efficient strategies, such as VS, to accelerate the identification of early hit molecules for TRPM8. However, there is a notable lack of published studies evaluating optimal VS methodologies for this target.

Building on this need for more efficient drug lead discovery, the present study focuses on defining best practices for structure-based VS against TRPM8. Specifically, we aim to determine which available TRPM8 structural model(s) offer the most reliable basis for such VS, which molecular docking protocol(s) yield the highest predictive performance, and how different docking methods can be optimally combined to enhance hit identification. By systematically addressing these questions, this work also introduces a robust and reproducible VS benchmark tailored to TRPM8.

Preparing virtual screening benchmark for TRPM8

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 (p Act; Equation 1), and the median pAct value for each compound was used as its final activity metric. Compounds with p Act ≥ 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.

Agonist and antagonist annotations were assigned using a combination of functional assay descriptions and activity type (IC 50, K d, and EC 50) information, depending on data availability. For compounds retrieved from ChEMBL and PubChem, functional assay descriptions were manually inspected to ensure that reported EC 50 and IC 50 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 50 values were used to identify agonists, while IC 50 values were used to identify antagonists. For compounds with multiple reported activity types, pEC 50 and pIC 50 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 50 values, compounds with pK d exceeding pEC 50 were classified as antagonists. The curated dataset was dominated by IC 50 measurements (∼2,300 compounds), a smaller number of EC 50 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.

Selection of TRPM8 structures

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.

Docking preparation

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 2+ 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 protocol

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 3 + 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.

Smina protocol

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.

Consensus protocol

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.

Hierarchical protocol

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).

Virtual screening performance metrics

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%.

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.

rDock is up to 4x times faster than Smina

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 also outperforms Smina in ranking actives higher than their decoys

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.

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.

The hierarchical protocol provides the best efficiency-hit rate compromise

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.

Correctly predicted TRPM8 modulators are mostly antagonist

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.

This study evaluated the suitability of two TRPM8 structural conformations (agonist-bound: TRPM8 7WRE and antagonist-bound: TRPM8 9B6G) and two docking protocols, Smina and rDock, for VS. Across the two TRPM8 structures, rDock showed consistent hit rates, whereas Smina performed better against TRPM8 9B6G than against TRPM8 7WRE, indicating that TRPM8 9B6G is the more suitable template for antagonist-focused screening with Smina.

In addition to performance differences, the two docking tools displayed distinct computational characteristics. rDock showed superior computational efficiency, with faster docking and minimal dependence on ligand molecular weight, whereas Smina docking time increased substantially with ligand size. These differences have practical implications for large-scale screening, particularly when computational resources are limited. At a 1% hit rate, rDock provided the strongest overall performance across targets, indicating its suitability for broad, cost-effective prioritization. However, at the more stringent 0.5% cutoff (representative of scenarios where only a very small number of compounds can be experimentally tested) consensus and hierarchical protocols consistently outperformed individual docking tools, demonstrating improved early enrichment. Notably, although rDock demonstrated stronger early recognition by placing more actives at first-rank, Smina achieved better overall ranking quality across targets. However, both docking tools showed considerable overlap between the score distributions of actives and decoys, indicating limited ability to reliably distinguish true actives solely based on docking scores.

The observation that rDock and Smina sample distinct regions of chemical space further supports the use of consensus and hierarchical protocols, as agreement between these diverse methods provides additional confidence in compound prioritization. Although consensus protocol retrieves more actives at 0.5% cutoff, this strategy is computationally expensive. In contrast, hierarchical screening improves the likelihood of selecting true actives without substantially increasing computational cost, thereby offering a practical balance between screening performance and efficiency.

Despite the use of both an agonist-bound and an antagonist-bound TRPM8 structure, and a ligand dataset containing compounds with both functional annotations, docking consistently retrieved only compounds that are experimentally classified as antagonists when functional annotations were examined post hoc. Notably, even when docking was performed against the agonist-bound TRPM8 structure and limited receptor flexibility was introduced in rDock, none of the 135 known TRPM8 agonists present in the dataset were recovered among the top-ranked compounds. This behavior likely reflects the strong class imbalance in the dataset (1,849 antagonists vs. 135 agonists), as well as inherent limitations of standard docking approaches in capturing the molecular features required for TRPM8 agonism. These observations highlight that docking scores primarily reflect binding compatibility rather than functional outcome.

Future work will therefore focus on developing target-specific SFs, which are expected to provide improved accuracy and functional prediction (Tran-Nguyen et al., 2023). In addition, predictive modelling strategies will be explored to anticipate which predicted binders will be agonist and which will be antagonist. Complementary analyses, such as the evaluation of ligand-based features, residue-level interaction patterns, and molecular dynamics simulations, could further improve functional discrimination and provide deeper insight into agonist- and antagonist-specific receptor interactions.