Towards high-resolution modeling of small molecule-ion channel interactions

Rosetta, ligand docking, ion channel, computational modeling, ligand – protein 15 interactions, computer-aided drug discovery. 16 17

18 Ion channels are critical drug discovery targets for a wider range of pathologies, such as epilepsy, 19 chronic pain, itch, autoimmunity, and cardiac arrhythmias. To develop effective and safe 20 therapeutics, it is necessary to design small molecules with high potency and selectivity for specific 21 ion channel subtypes. There has been increasing implementation of structure-guided drug design for 22 the development of small molecules targeting ion channels. We evaluated the performance of two 23 Rosetta ligand docking methods, RosettaLigand and GALigandDock, on structures of known ligand - 24 cation channel complexes. Ligands were docked to voltage-gated sodium (NaV), voltage-gated 25 calcium (CaV), and transient receptor potential vanilloid (TRPV) channel families. For each test case, 26 RosettaLigand and GALigandDock methods were able to frequently sample a ligand binding pose 27 within 1-2 Å root mean square deviation (RMSD) relative to the native ligand – channel structure. 28 However, RosettaLigand and GALigandDock scoring functions cannot consistently identify native 29 drug binding poses as top-scoring models. Our study reveals that the proper scoring criteria for 30 RosettaLigand and GALigandDock modeling of ligand - ion channel complexes should be assessed 31 on a case-by-case basis using sufficient ligand and receptor interface sampling, knowledge about 32 state specific interactions of the ion channel and inherent receptor site flexibility that could influence 33 ligand binding. 34 35

36 The voltage-gated cation channel families consist of pore-forming transmembrane proteins that 37 selectively conduct ions across lipid bilayers, and mediate physiological processes such as signal 38 transduction, gene expression, synaptic transmission, and the activation and proliferation of cells in 39 the immune system (Catterall, 1995; Hille, 2001; Catterall, 2011; Feske, Wulff, and Skolnik, 2015; 40 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Nanou and Catterall, 2018). Cation channels function in a finely regulated manner across spatial and 41 temporal domains to complete these cellular functions. Current drug discovery efforts aim to 42 modulate channel activity by targeting specific channel domains. For instance, a voltage-gated 43 sodium channel's therapeutically relevant structural domains are the selectivity filter (otherwise 44 known as the outer pore vestibule), the central pore cavity (otherwise known as the inner pore 45 vestibule), and the voltage-sensing domain (Nguyen and Yarov-Yarovoy, 2022). There have been 46 considerable academic and industry efforts to identify therapeutically relevant small molecules that 47 selectively target ion channels (Bagal et al., 2013; Zhang et al., 2020). However, developing effective 48 and safe therapeutics targeting ion channels has been challenging (Wulff et al., 2019). 49 To address these challenges, drug discovery pipelines are trending towards incorporation of 50 virtual drug screening and computer-aided drug design processes, for their ability to minimize drug 51 development time and cost (Maia et al, 2020). Among these processes, molecular docking has 52 demonstrated its usefulness for structure-based drug discovery. Molecular docking involves 53 predicting conformations of the small molecule’s orientation with the protein (known as the pose) 54 and scoring the poses to rank the likely protein-ligand interaction (Meng et al, 2011). 55 Among the numerous molecular docking software packages, Rosetta is a protein modeling 56 software and design suite with two established small molecule docking methods: RosettaLigand 57 (Meiler and Baker, 2006; Davis and Baker, 2009; Smith and Meiler, 2020) and GALigandDock (Park 58 et al, 2021). RosettaLigand uses a Monte Carlo minimization procedure using the Rosetta energy 59 function (Alford et al, 2017) to dock a pre-generated set of ligand conformers, while allowing side 60 chain flexibility within a protein-receptor site. GALigandDock utilizes a different approach with two 61 distinct features. First, the scoring function, RosetteGenFF, is a new generalized energy function 62 tailored for small molecules. RosettaGenFF was trained from the Cambridge Structural Database 63 (Groom et al, 2016), which at the time contained 1,386 small molecule crystal lattice arrangements to 64 create a balanced force field that discriminates true lattice packing arrangements of the ligand from 65 decoy (alternative lattice packing and conformational) arrangements. During docking, an orientation-66 dependent water-bridging energy term is incorporated within RosettaGenFF to further discriminate 67 the protein-ligand orientation (Pavlovicz, Park, and DiMaio, 2020). Second, GALigandDock samples 68 conformational space using a genetic algorithm. The ligand rigid body degrees of freedom and 69 rotatable torsions are encoded as ‘genes’ to generate new ligand inputs for successive docking 70 iterations. This allows efficient sampling of the protein-ligand energetic landscape when paired with 71 the RosettaGenFF score function, canonical Monte Carlo optimization, and quasi-Newtonian 72 minimization procedures within the Rosetta framework (Park et al, 2021). 73 While Rosetta protein-ligand docking methods are able to perform well with soluble protein – 74 ligand benchmarks, the application of these methods to membrane-embedded ion channels has not 75 been explored. Since there is a need to better assess and screen small molecules targeting different 76 ion channel domains, we selected a diverse set of ten known cation channel-ligand structures for 77 evaluation. From this set of high-resolution ion channel – small-molecule complexes, we assessed the 78 accuracy of the RosettaLigand and GALigandDock methods in sampling ligand poses near the native 79 structural coordinates and predicting the closest matching pose by energy ranking. 80 Our case studies include four voltage-gated sodium (NaV) channel structures, five voltage-gated 81 calcium (CaV) channel structures, and one transient receptor potential vanilloid (TRPV) channel 82 structure. The ion channel ligand binding sites include the voltage-sensing domain, the selectivity 83 filter, and the central pore cavity. Our results demonstrate that RosettaLigand and GALigandDock 84

can frequently sample ligand binding poses within 1-2 Å root-mean-square deviation 85 (RMSD) from the native channel – ligand structure. However, the ability to identify a near-native 86 pose from energy ranking remains a challenge. When considering factors like the targeted ion 87 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint channel domain, the ligand library features, and the sampling of ligand and receptor site 88 conformations, our work demonstrates that high-resolution structures paired with RosettaLigand or 89 GALigandDock can support drug discovery pipelines on a case-by-case basis. 90 91

92 Ligand generation 93 Ligands were extracted as Structure Data Files (.sdf) from PubChem (Kim et al, 2023). Using the 94 Avogadro software (Hanwell et al, 2012), each ligand structure underwent bond-correction, 95 protonation at pH 7.4, and energy minimization using the Merck Molecular Force Field (Halgren, 96 1996a; b; c; Halgren and Nachbar, 1999; Halgren, 1999a; b). The resulting models were saved as 97 Tripos Mol2 (.mol2) files. The protonation and bond order of saxitoxin and tetrodotoxin was matched 98 to experimentally reported work (Hinman and DuBois, 2003; Thomas-Tran and Du Bois, 2016). Both 99 experimentally resolved structures of verapamil docked to rabbit CaV1.1 were tested (Zhao et al, 100 2019). 101 Next, using Amber Tools’ Antechamber protocol, the partial atomic charge, atom, and bond type 102 assignments for each ligand were AM1-BCC corrected (Case et al, 2021; Salomon-Ferrer, Case, and 103 Walker, 2012; Appendix S1). AM1-BCC correction is commonly used in Rosetta-based ligand 104 docking protocols (Smith and Meiler, 2020; Park et al, 2021) and has demonstrated a similar 105 performance correlation with other RosettaLigand input preparation protocols (Smith and Meiler, 106 2020). 107 The AM1-BCC corrected ligands were used for subsequent steps specific to each method. For 108 RosettaLigand, an in-house script (Appendix S2) using the OpenEye Omega toolkit (Hawkins et al, 109 2010) was used to generate the conformer library, followed by Rosetta to generate the associated 110 ligand parameters file. For GALigandDock, the input conformer was generated using the 111 RosettaGenFF crystal structure prediction protocol (Park et al, 2021; Appendix S3), taking the 112 lowest energy packing arrangement as input. 113 Ion channel preparation 114 Ion channel structures were downloaded from the Protein Data Bank (Berman et al, 2000). Prior 115 to RosettaLigand docking, structures were relaxed with backbone constraints using the RosettaRelax 116 protocol (Nivón, Moretti, and Baker, 2013). This protocol allows the repacking of protein sidechains 117 and minimization of the structure into the Rosetta score function for comparison between poses. The 118 lowest energy pose from 100 relaxed poses was used for docking. 119 RosettaLigand docking 120 RosettaLigand docking was performed using previously described RosettaScripts protocols 121 (Davis and Baker, 2009; DeLuca, Khar, and Meiler, 2015; Appendix S4-S7). Briefly, the initial 122 placement of all ligands into their respective ion channels was performed by superimposing the initial 123 ligand poses onto the experimental ligand coordinates. RosettaLigand uses grid-based sampling to 124 score the ion channel – ligand interface, whereby all ligand atoms and ion channel atoms a set 125 distance from the ligand are scored with the defined scoring function. The scoring grid width was 126 calculated uniquely for each ligand to ensure all ligand atoms were bound to the scoring grid. As 127 used previously, the scoring grid width was calculated as the maximum conformer atom-atom 128 distance plus twice the box size value used in the Transform mover (Moretti, Bender, Allison, and 129 Meiler, 2016; Appendix S4). 130 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint RosettaLigand has a low-resolution and high-resolution sampling phase. For an unbiased 131 sampling of the domain, an initial transformation during the low-resolution sampling phase is 132 performed on the initial ligand position using Rosetta’s Transform mover with a box size of 7-8 Å. 133 The low-resolution Transform step is a grid-based Monte Carlo simulation, where the ligand is 134 translated up to 0.2 Å and rotated up to 20 degrees per iteration, for a total of 500 iterations. The pose 135 with the lowest score is then used as the starting pose for high-resolution docking. For high-136 resolution docking using Rosetta’s HighResDocker mover, six docking cycles of rotamer 137 sampling are performed with repacking of sidechains every third iteration. Lastly, the ion channel – 138 ligand complex is minimized using Rosetta’s FinalMinimizer mover, and the interface scores 139 are reported using the InterfaceScoreCalculator. 140 For all RosettaLigand docking runs, the ligand_soft_rep and hard_rep scoring functions 141 were reweighted based on previous work assessing Rosetta score functions with the Comparative 142 Assessment of Score Function 2016 (CASF-2016) dataset (Su et al., 2019; Smith and Meiler, 2020). 143 Specifically, reweights of Coulombic electrostatic potential (fa_elec), Lennard-Jones repulsive 144 energy between atoms in the same residue (fa_intra_rep), sidechain-backbone hydrogen bond 145 energy (hbond_bb_sc), sidechain-sidechain hydrogen bond energy (hbond_sc), and 146 Ramachandran preferences (rama) were applied for the soft-repulsive and hard-repulsive docking 147 phases (Appendix S6). 148 For each docking run, either 20,000 poses or 100,000 poses were generated to assess if there are 149 statistically significant differences in nearest native pose RMSD. In RosettaLigand, a ligand interface 150 is defined either by a representative ligand atom (a “neighbor atom”, defined as the geometric center 151 of mass by default), or all ligand atoms, relative to all ion channel Cβ atoms within a specified radius 152 from the ligand (commonly default to 6 or 7 Å). For RosettaLigand, two mutually exclusive ligand 153 area interface modes are available for scoring the pose interface: the ligand neighbor atom cutoffs 154 mode (add_nbr_radius=”True”), and the all ligand atom cutoffs mode 155 (all_atom_mode=”True”). In previous work, both modes were used (Moretti, Bender, Allison, 156 and Meiler, 2016; Smith and Meiler, 2020), thus, both modes were evaluated for any differences in 157 performance. Four individual RosettaLigand docking sets were performed for each PDB structure, by 158 combining different pose totals and ligand area interface modes, for performance comparison. 159 GALigandDock docking 160 As described in the original study, docking was performed using the RosettaScripts’ 161 GALigandDock mover in DockFlex mode (Park et al, 2021; Appendix S8-S10). Replicate runs of 162 GALigandDock were performed in parallel for each structure evaluated. Each run consisted of 20 163 generations with a pool of 100 poses, where each generation updates the pool by total energy. By 164 default, GALigandDock outputs the top 20 structures from the final generation, however, for this 165 study, the entire pool of poses were used. For each docking run, 20,000 poses (1,000 runs) or 166 100,000 poses (5,000 runs) were output with a padding value of 2 Å, 4 Å, or 7 Å to test for statistical 167 differences in the lowest RMSD pose (RMSDMin). In total, six individual GALigandDock docking 168 sets were performed for each PDB structure, by combining different pose totals and padding sizes, 169 for performance comparison. 170 PNear calculation 171 PNear is a quantitative metric to evaluate an energy funnel’s quality from a population of poses, 172 with the lowest-scoring pose as the converged minima, or native state. PNear was calculated as 173 previously defined (Bhardwaj, Mulligan, Bahl, et al., 2016), and applied as previously described for 174 ligand docking (Smith and Meiler, 2020), with the native state being the ligand structure coordinates. 175 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint PNear ranges from 0 (protein will not converge to native state) to 1 (protein will always converge to 176 the native state). 177 𝑃!"#$ = ∑ 𝑒%$&'(! " )"! *+, 𝑒 % -! .#/ ∑ 𝑒 % -$ .#/! 0+, 185 The numerator is the Boltzmann probability of an individual pose being near the native state, 178 governed by the “native-ness” parameter (λ), and the thermal energy, the product of the Boltzmann 179 constant and absolute temperature (kBT). The denominator is the partition function of a canonical 180 ensemble. For small-molecule docking, the root mean square deviation (RMSD) of the pose’s ligand 181 coordinates relative to the structure’s ligand coordinates when the protein pose is superimposed to the 182 entire protein structure, or the receptor site is used. The energy scoring metric (Ei) is the interface 183 energy: the energy solely composed of protein-ligand interactions. 184 A previous Rosetta small-molecule docking study (Park et al., 2021) categorized native-like 186 models at 1 Å or 2 Å RMSD without calculating PNear, while another study evaluating RosettaLigand 187 performance calculated PNear using λ =1.5 Å, and kBT =0.62 (Meiler and Smith, 2020). Therefore, we 188 calculated PNear with native-like models defined by λ=2.0 Å and kBT =0.62 Rosetta energy units 189 (REU), however we calculated PNear using all previously reported values for the parameters λ and kBT 190 (Tables S4.1-4.10). A specific PNear cutoff indicative for drug discovery pipelines has not been 191 established, hence we will refer to PNear ≥ 0.5 as a ‘first pass’ cutoff for this study when evaluating 192 energy funnel convergence to the ligand structure coordinates. 193 Statistical tests 194 Tests for normality, heteroscedasticity, and Pearson’s correlation between covariates were 195 performed in Python using NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), and Pingouin 196 (Vallat, 2018) with a significance level (α) of 0.05. Population data consisted of the RMSDMin for 197 each docking set. Not all RMSDMin data for each docking set fit a normal distribution. Since 198 RMSDMin data should skew towards zero Å, a lognormal base 10 transformation was applied to all 199 RMSDMin population data when comparing across docking sets. Shapiro Wilk tests for normality and 200 Levene heteroscedasticity tests were performed to ensure transformed data were normal and of equal 201 variance, respectively. 202 Covariates for both methods were the number of rotatable bonds of the ligand, the number of 203 ligand heavy atoms, and the resolution of the structure. Statistical tests for RosettaLigand also 204 included the number of conformers generated, transformed with a lognormal base 10, as a covariate. 205 The ligand molecular weight was initially included as a covariate but was discarded after identifying 206 strong Pearson’s correlation with the number of ligand heavy atoms (Table S3, Figure S1.). 207 To assess RMSDMin across all docking sets for a given method, a repeated measures factorial 208 analysis of variance (ANOVA) with covariates were performed using IBM SPSS (version 29). For 209 RosettaLigand, the two factors were sample size and ligand area interface mode, while for 210 GALigandDock the two factors were sample size and padding value. Mauchly’s test for sphericity 211 was performed for levels of padding (p = 0.63). 212 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint

213 Table 1. Summary of ligand – ion channel structures docked, domain, and covariates included in the study. 214 PDB ID Description Ion channel modulation / therapeutic class Domain Total ligand rotatable bonds Total ligand heavy atoms Ligand molecular weight (Da) Structure resolution (Å) Total RosettaLigand Conformers 5EK0 human NaV1.7- VSD4-NaVAb / GX-936 Inhibitor / NA Voltage- sensing domain 8 39 575.6 3.53 2179 6J8G human NaV1.7 / Saxitoxin Blocker / NA Selectivity filter 3 21 299.3 3.20 28 6J8I human NaV1.7 / Tetrodotoxin Blocker / NA Selectivity filter 1 22 319.3 3.20 2 6JP5 rabbit CaV1.1 / Nifedipine Blocker / vasodilator, antihypertensive, antianginal CPC† 5 25 346.3 2.90 22 6JP8 rabbit CaV1.1 / (S)-(−)-Bay K 8644 Activator / NA CPC 3 25 356.3 2.70 3 6JPA rabbit CaV1.1 / (S)-Verapamil Blocker / antiarrhythmic, vasodilator CPC 13 33 454.6 2.60 785 6JPB rabbit CaV1.1 / Diltiazem Blocker / antihypertensive, vasodilator CPC 7 29 414.5 2.90 45 6KZP human CaV3.1 / Z944 Blocker / NA CPC 6 26 383.9 3.10 1992 6U88 rat TRPV2 / Cannabidiol Activator / anticonvulsants CPC 6 23 314.5 3.20 71 6UZ0 rat NaV1.5 / Flecainide Blocker / antiarrhythmic class 1c CPC 7 28 414.3 3.24 1169 † CPC = Central pore cavity.215 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint 216 In this study, we evaluated ligands that are suggested to modulate ion channel activity for a given 217 clinical effect, such as antiarrhythmic, anticonvulsant, and antihypertensive (Table 1). Our study also 218 includes canonical ion channel blockers, such as tetrodotoxin and saxitoxin, small molecules whose 219 structures have been resolved for drug discovery campaigns, such as GX-936, Z944, and Bay K 220 8544, and drugs currently approved for therapy, such as nifedipine, diltiazem, and flecainide (Figure 221 1).222 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint 223 Figure 1. Two-dimensional structures of ligands docked in this study. The depicted 224 stereochemistry reflects that resolved from structures. The protonation and bond order of saxitoxin 225 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 9 and tetrodotoxin from prior reported work was used (Hinman and DuBois, 2003; Thomas-Tran and 226 Du Bois, 2016). 227 228 Performance criteria for Rosetta ligand docking 229 A previous study has emphasized that two criteria must be satisfied for accurate modeling of 230 protein-ligand interactions (Kaufmann and Meiler, 2012). First, the method must produce native-like 231 poses through sufficient sampling. For small molecule docking, a native-like ligand pose has a root 232 mean square deviation (RMSD) below a context-dependent predetermined value. Previously reported 233 RMSD ranges when assessing a native-like pose are 1.0 Å or 2.0 Å, with 2.0 Å being a common 234 cutoff for drug discovery pipelines (Maia et al., 2020; Park et al., 2021). When calculating RMSD, 235 the protein pose is superimposed to the entire protein structure or the receptor site, and then the 236 ligand coordinates of the pose are evaluated relative to the ligand coordinates of the structure. 237 Second, the pose population must produce an energy funnel converging to the lowest-energy, 238 native-like pose (Meiler and Baker, 2006). Rather than the lowest-energy pose, we use empirical 239 structural evidence in this study to define the native pose as the ligand’s structural coordinates. An 240 energy funnel is qualitatively evaluated using the interface energy between ligand and protein with 241 respect to ligand RMSD (Meiler and Baker, 2006, Smith and Meiler, 2020). A quantitative metric of 242 funnel convergence, PNear, has been utilized in Rosetta protein design (Guarav, Mulligan, Bahl, et al., 243 2016; Mulligan et al., 2021), and has been adopted for small molecule docking (Smith and Meiler, 244 2020). Due to the large number of poses, funnel quality is often assessed with a population subset 245 consisting of the best-scoring poses (Combs et al., 2010; Lemmon, Kaufmann, and Meiler, 2012; 246 Shim et al., 2019). 247 With these criteria in mind, we evaluated the performance of RosettaLigand and GALigandDock 248 for the following ion channel receptor sites: the voltage-sensing domain, the pore-forming domain, 249 and the central pore cavity (Table 1). We evaluated RosettaLigand using combinations of sample 250 size (20,000 poses vs 100,000 poses) and ligand area interface mode (ligand neighbor atom vs all 251 ligand atom). We evaluated GALigandDock using combinations of sample size (20,000 poses vs 252 100,000 poses) and padding of the sampling grid (2 Å, 4 Å, and 7 Å). Increasing the padding of the 253 sampling grid enables additional rotamer sampling around the receptor and increases translational 254 sampling of the ligand around the receptor site.255 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Table 2. The minimum RMSD (Å) pose for each ligand-channel docked that is detailed in the 256 primary text. 257 PDB ID Description RosettaLigand, all ligand atom interface cutoff, 100,000 poses GALigandDock, padding 7 Å, 100,000 poses 5EK0 human NaV1.7-VSD4- NaVAb / GX-936 0.91 0.65 6J8G human NaV1.7 / Saxitoxin 0.70 0.76 6J8I human NaV1.7 / Tetrodotoxin 0.54 0.81 6JP5 rabbit CaV1.1 / Nifedipine 0.77 1.3 6JP8 rabbit CaV1.1 / (S)-(−)- Bay K 8644 0.93 0.97 6JPA-1† rabbit CaV1.1 / (S)-Verapamil 1.4 2.0 6JPA-2 rabbit CaV1.1 / (S)-Verapamil 2.2 2.1 6JPB rabbit CaV1.1 / Diltiazem 2.5 1.1 6KZP human CaV3.1 / Z944 0.73 0.84 6U88 rat TRPV2 / Cannabidiol 1.0 0.90 6UZ0 rat NaV1.5 / Flecainide 1.2 0.94 Average 1.2 ± 0.6 1.1 ± 0.5 Median 0.93 0.94 † Verapamil was resolved with two binding orientations to the rabbit CaV1.1 central pore cavity in 258 6JPA and are named 6JPA-1 and 6JPA-2 in this study respectively. 259 260 261 Covariates potentially influencing RMSD 262 We controlled for the following covariates that are dependent upon the ligand or structure used 263 for docking. We speculated that the number of ligand rotatable bonds and the number of heavy atoms 264 could influence RMSDMin by increasing the amount of internal sampling required during the docking 265 run. We also speculated that a poorer, higher reported structural resolution could result in an 266 increased overall RMSDMin. The total number of ligand heavy atoms was also used as a covariate for 267 RMSDMin, which is strongly correlated to the molecular weight of each ligand (Figure S1). Lastly, 268 for RosettaLigand, the total number of conformers provided as input could affect the RMSDMin since 269 more conformers could inherently require more sampling. 270 271 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 11 RosettaLigand and GALigandDock generate native-like ion channel-ligand poses 272 This study evaluated docking sets using population data consisting of the lowest RMSD pose 273 (RMSDMin) from each channel – ligand test case. We provide four docking sets for RosettaLigand 274 and six docking sets for GALigandDock, using combinations of the sample size with either the 275 RosettaLigand-specific ligand area interface calculation, or GALigandDock-specific padding value. 276 For brevity, we will discuss the RosettaLigand docking set using a sample size of 100,000 poses and 277 the all-ligand atoms interface mode, and the GALigandDock docking set using a sample size of 278 100,000 poses and padding value of 7 Å. These specific docking sets provide the greatest breadth of 279 receptor sampling and scoring among our docking sets by appropriately encompassing the entirety of 280 each ion channel domain tested (Table 2). The results for other docking sets are provided in the 281 Supplemental Material (Tables S1.1-1.2). 282 For RosettaLigand, the average RMSDMin across the docking set was 1.2 ± 0.6 Å (Table 2). 283 When comparing RosettaLigand docking sets and using repeated measures ANOVA with covariates, 284 an interaction from both factors (sample size and ligand area interface mode) did not result in 285 rejection of the null hypothesis of equivalent means from common logarithm transformed RMSDMin 286 data (p=0.71). Likewise, individual interactions from sample size (p=0.38) and ligand area interface 287 mode (p=0.72) did not result in rejection of the null hypothesis of equivalent means. Each 288 RosettaLigand docking set produced comparable RMSDMin averages and standard deviations within a 289 few sub-Angstroms (Table S1.1). Therefore, our results suggest that using either ligand area 290 interface mode with 20,000 or 100,000 total poses could generate similar near-native models. 291 For GALigandDock, the average RMSDMin across the docking set was 1.1 ± 0.5 Å (Table 2), 292 while all docking sets produced similar RMSDMin averages and standard deviations within a few sub-293 Angstroms (Table S1.2). Again, using repeated measures ANOVA with covariates, an interaction 294 from both factors (sample size and padding size) failed to result in a rejection of the null hypothesis 295 of equivalent means (p=0.91). Likewise, individual interactions from sample size (p=0.59) and 296 padding size (p=0.34) did not result in rejection of the null hypothesis of equivalent means. While 297 there was no statistical advantage to using a specific padding value to achieve a lowered average 298 RMSDMin, we note that ion channel structures within the central pore cavity had greater sampling 299 with a padding value of 7 Å, since a padding value of 5 Å did not encompass the entire pore. 300 Therefore, the appropriate padding size is context dependent and should be adjusted when using 301 GALigandDock to provide sufficient sampling grid space.302 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Table 3. PNear of RMSDMin and interface energy for each ligand-channel docked that is detailed in the 303 primary text, calculated with kBT=0.62 and λ=2.0. 304 RosettaLigand, all ligand atom interface cutoff, 100,000 poses GALigandDock, padding 7 Å, 100,000 poses PDB ID Description Full population Top 10% total_score Full population Top 10% total_score 5EK0 human NaV1.7- VSD4-NaVAb / GX-936 0.22 0.26 0.83 0.83 6J8G human NaV1.7 / Saxitoxin 1.2•10-2 9.3•10-3 1.6•10-2 1.6•10-2 6J8I human NaV1.7 / Tetrodotoxin 0.12 8.8•10-2 0.27 0.27 6JP5 rabbit CaV1.1 / Nifedipine 0.21 0.28 0.24 0.29 6JP8 rabbit CaV1.1 / (S)-(−)-Bay K 8644 0.62 0.66 0.60 0.61 6JPA- 1 rabbit CaV1.1 / (S)-Verapamil 4.7•10-3 4.3•10-3 2.4•10-3 2.8•10-3 6JPA- 2 rabbit CaV1.1 / (S)-Verapamil 5.0•10-4 1.2•10-4 1.8•10-5 1.9•10-5 6JPB rabbit CaV1.1 / Diltiazem 6.8•10-4 6.4•10-4 2.8•10-6 7.8•10-6 6KZP human CaV3.1 / Z944 0.50 0.62 1.2•10-4 1.1•10-4 6U88 rat TRPV2 / Cannabidiol 8.3•10-2 4.3•10-2 5.0•10-3 3.8•10-3 6UZ0 rat NaV1.5 / Flecainide 0.14 0.12 5.9•10-6 2.9•10-5 305 306 RosettaLigand and GALigandDock energy funnels for ion channel-ligand docking 307 We evaluated whether the entire population of generated poses and the top 10% of the lowest 308 total energy-scoring poses would achieve PNear values indicative of an energy funnel converging onto 309 the ligand structural coordinates (Tables S4.1-4.10). Since interface energy describes how a ligand is 310 interacting with the ion channel, PNear was calculated with the interface energy rather than the total 311 energy. A specific PNear cutoff indicative for drug discovery pipelines has not been established. 312 Hence, we will refer to PNear ≥ 0.5 as a ‘first pass’ cutoff for this study when evaluating energy funnel 313 convergence to the ligand structure coordinates. For brevity, we will discuss PNear with a “native-314 ness” parameter (λ) of 2.0, and a thermal energy parameter (kBT) of 0.62. Other PNear values 315 matching parameter values reported in other work are provided in the Supplemental Material 316 (Bhardwaj, Mulligan, Bahl, et al., 2016; Smith and Meiler, 2020; Tables S4.1-4.10). 317 The PNear of the full population and the top 10% of total energy-scoring poses were similar within 318 methods, with the bulk of PNear values being in the thousandths or lower (Table 3). However, 319 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 13 RosettaLigand and GALigandDock were able to identify energy funnels for unique cases. For 320 RosettaLigand, the human CaV3.1- Z944 complex (PDB: 6KZP) achieved PNear ≥ 0.5, while for 321 GALigandDock, human NaV1.7-VSD4-NaVAb - GX-936 (PDB: 5EK0) achieved PNear ≥ 0.5. For 322 both methods, rabbit CaV1.1-Bay K8644 (PDB: 6JP8) achieved PNear ≥ 0.5. 323 Following the standard of reporting a percentage of poses (Combs et al., 2010; Lemmon, 324 Kaufmann, and Meiler, 2012; Shim et al., 2019), specific PNear values reported herein refer to the top 325 10% of total energy-scoring poses (Table 3). 326 Notably, there are visual distinctions in interface energy funnel plots between RosettaLigand and 327 GALigandDock (Figures S2-S12). Generally, RosettaLigand can sample more poses up to 2 Å 328 RMSD compared to GALigandDock, but with a less distinguishable energy funnel since some poses 329 score greater than zero, indicating an unfavorable energy score. Since GALigandDock uses the 330 lowest total energy poses as input for future conformer generations, poses with interface energy 331 greater than zero are infrequent. 332 Lastly, when using the lowest interface energy pose from the full docking population as ranking 333 criteria, the RMSDMin increases to a range between 1.2 Å and 10.8 Å (mean 4.5 ± 2.8 Å) with 334 RosettaLigand (Table S2.1), and a range between 0.83 Å and 14.5 Å (mean 6.4 ± 5.0 Å) with 335 GALigandDock (Table S2.5). This suggests that extracting the lowest energy pose is not a reliable 336 indicator of a near-native pose and does not necessarily reflect the native binding mode for ion 337 channel-ligand docking. 338 Ligand Docking into the voltage-sensing domain 339 The only ion channel-ligand structure evaluated for ligand docking into the voltage-sensing 340 domain (VSD) was GX-936 in complex with the VSD4 of the hNaV1.7-NaVAb chimera (PDB: 341 5EK0). The hNaV1.7 channel has been validated as a drug target for pain signaling, and aryl 342 sulfonamides have been reported as selective inhibitors. Specifically, GX-936 exhibits selectivity for 343 hNaV1.7 compared to other hNaV isoforms (McCormack et al., 2013; Ahuja et al., 2015; Nguyen and 344 Yarov-Yarovoy, 2022). 345 After docking GX-936 to the hNaV1.7 VSD4, the RMSDMin poses were 0.91 Å using 346 RosettaLigand and 0.65 Å using GALigandDock (Figure 2, Table 2). The largest deviations from 347 the native structure were observed at the peripheral pyrazole ring and the ethyl azetidine (Figure 2). 348 Notably, GX-936 has eight rotatable bonds, making it the second most flexible ligand used in this 349 study (Table 1; Figure 1). However, the sampled ligand poses were consistently below 2 Å RMSD 350 for each docking run (Table S1.1-1.2). RosettaLigand was unable to identify the RMSDMin pose as 351 the lowest interaction energy pose (Table S2.1), scoring the lowest energy pose 2.7 Å from the 352 structure coordinates. This pose has the same binding orientation, but the pyrazole ring clashes with 353 the structurally resolved lipid, while the ethyl azetidine withing the VSD is orientated up towards the 354 extracellular space, rather than down. Indeed, the lowest ten interface energy poses all possess these 355 features (Figure S13.1). GALigandDock scored the lowest energy pose 0.83 Å from the structure 356 coordinates (Table S2.5). Compared to RosettaLigand, the lowest ten interface energy poses 357 converge well with the RMSDMin pose, with only one of the ten poses clashing with the lipid (Figure 358 S13.2). Further, GALigandDock was able to discriminate with good confidence an energy funnel 359 (PNear = 0.83; Table 3), whereas RosettaLigand was unable to distinguish a clear energy funnel (PNear 360 = 0.26; Table 3). 361 362 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 14 363 Figure 2. GX-936 docking to hNa V1.7-NaVAb voltage sensing domain 4 (PDB 5EK0). (A) 364 RosettaLigand (RMSDMin = 0.91 Å, PNear = 0.26) and (B) GALigandDock (RMSDMin = 0.65 Å, PNear 365 = 0.83). Left: transmembrane view . Middle: extracellular view . Right: ligand RMSD vs interface 366 energy distribution of the top 10% of total energy poses. Grey dots: 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 368 GALigandDock in dark pink, native structure coordinates in dark grey and phospholipid head group 369 resolved from the structure in light grey. Non-carbon atoms match the Corey-Pauling-Koltun coloring 370 scheme. Black star indicates RMSDMin pose from the entire population. 371 372 Ligand Docking into the Pore-Forming Domain 373 Tetrodotoxin (TTX) and saxitoxin (STX) are both guanidinium-based small molecules derived 374 from puffer fish and shellfish, and act as selective blockers of sodium channels by binding to the 375 pore-forming domain (Hille, 2018; Shen et al., 2019). Being potent pore blockers for only some NaV 376 channel isoforms, hNaV channel isoforms are classified in physiology as TTX-resistant (hNaV1.5, 377 hNaV1.8-1.9) and TTX-sensitive (hNaV1.1-1.4, hNaV1.6-1.7) (Stevens, Peigneur, and Tytgat, 2011). 378 Like aryl sulfonamides for VSDs, the discovery of TTX and STX binding to the pore-forming 379 domain of NaV channels has spurred the design of similar blockers with greater selectivity for a 380 specific hNaV isoform, usually hNaV1.7 for pain therapy (Hagen et al., 2017; Pajouhesh et al., 2020). 381 We chose STX and TTX as test cases to evaluate ligand docking into the pore-forming domain of 382 hNaV1.7 (PDB: 6J8G and 6J8I, respectively; Shen et al., 2019). The protonation and bond order of 383 STX and TTX from a previously reported work was used (Thomas-Tran and Du Bois, 2016). Based 384 on their cage-like structures, STX has only 3 rotatable bonds and TTX has only 1 rotatable bond: 385 making them the most rigid ligands in this study (Table 1; Figure 1). In both cases, RosettaLigand 386 docking resulted in RMSDMin values of 0.70 Å for STX and 0.54 Å for TTX, while GALigandDock 387 reported RMSDMin values of 0.76 Å for STX and 0.81 Å for TTX (Table 2; Figures 3 and 4). The 388 PNear values suggested little to no energy funnel convergence with STX (RosettaLigand PNear = 389 9.3•10-3 and GALigandDock PNear = 1.6•10-2), and TTX (RosettaLigand PNear = 7.8•10-2 and 390 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 15 GALigandDock PNear = 0.27) (Table 3). This lack of convergence is due to other energy minima 391 occurring within 3-6 Å RMSD (Figures 3 and 4). 392 For STX, the lowest ten interface energy poses with RosettaLigand converged at 4.4 Å RMSD, 393 with the carbamoyl group positioned deeper into the selectivity filter, while GALigandDock rendered 394 a 3.8-6.4 Å RMSD, with the center of STX further away from the selectivity filter in various 395 rotations (Figure S14.1-2). For TTX, the lowest ten interface energy poses with RosettaLigand 396 contained poses between 2.9-4.2 Å RMSD; five of ten poses converged to the structural coordinates 397 with the amino group pointing away from the selectivity filter, while the other five poses converged 398 to the structural coordinates with the amino group pointing perpendicular to the selectivity filter path, 399 (Figure S15.1). With GALigandDock the poses were between 2.2-5.9 Å RMSD, with five of ten 400 poses having the TTX amino group pointing towards the pore and one set of P1/P2 helices, three of 401 ten poses with the TTX amino group pointing toward the selectivity filter path, and one pose with the 402 TTX amino group pointing toward the extracellular space. (Figure S15.2). 403 404 Figure 3. Saxitoxin docking to hNaV1.7 selectivity filter (PDB 6J8G). (A) RosettaLigand 405 (RMSDMin = 0.70 Å, PNear = 9.3•10-3) and (B) GALigandDock (RMSDMin = 0.76 Å, PNear = 1.6•10-2). 406 Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 407 distribution of the top 10% of total energy poses. Grey dots: 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 409 in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-410 Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 411 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 16 412 Figure 4. Tetrodotoxin docking to hNaV1.7 selectivity filter (PDB 6J8I). (A) (RMSDMin = 0.54 Å, 413 PNear = 7.8•10-2) and (B) GALigandDock (RMSDMin = 0.81 Å, PNear = 0.27). Left: transmembrane view. 414 Middle: extracellular view. Right: ligand RMSD vs interface energy distribution of the top 10% of total 415 energy poses. Grey dots: 4 Å. Carbon atoms 416 of RosettaLigand molecule are shown in cyan, GALigandDock in dark pink, and native structure 417 coordinates in dark grey. Non-carbon atoms match the Corey-Pauling-Koltun coloring scheme. Black 418 star indicates RMSDMin pose from the entire population. 419 420 Ligand docking into the central pore cavity 421 Seven test cases were used to evaluate ligand docking into the central pore cavity involving the 422 following channels: four cases from rabbit CaV1.1, one from human CaV3.1, one from rat NaV1.5, 423 and one from TRPV2. The central pore cavity is a broad target, with small molecules reported to 424 traverse through the gate or fenestration sites to reach their binding site and modulate channel activity 425 via pore blockade or allosteric mechanism (Hille, 1977; Hockerman et al., 1997; Jiang et al., 2020). 426 Drugs targeting this region can act as vasodilators (dihydropyridines), antiarrhythmics 427 (benzothiazepines, phenylalkylamines, flecainide), antiepileptics (Z944, cannabidiol), or local 428 anesthetics (flecainide) (Pumroy et al., 2019; Zhao et al., 2019a; Zhao et al., 2019b; Jiang et al., 429 2020). 430 We docked nifedipine (dihydropyridine channel blocker), Bay K 8644 (dihydropyridine channel 431 activator), diltiazem (benzothiazepine channel blocker), and two conformations of verapamil 432 (phenylalkylamine channel blocker) into rabbit CaV1.1. The number of rotatable bonds was five, 433 three, seven, and thirteen, respectively (Table 1; Figure 1). 434 Nifedipine docking resulted in a RosettaLigand RMSDMin of 0.77 Å and a GALigandDock 435 RMSDMin of 1.3 Å (Figure 5, Table 2). The calculated PNear suggested little energy funnel 436 convergence (RosettaLigand PNear = 0.28 and GALigandDock PNear = 0.29) (Table 3). This lack of 437 energy convergence is exemplified by the interaction energy distributions containing multiple low 438 energy minima 1-2 Å and 4-6 Å RMSD away from the ligand’s structural coordinates (Figure 5). 439 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 17 The lowest ten interface energy poses with RosettaLigand contained poses between 1.1-4.4 Å 440 RMSD, with most poses being 1.7 Å or less; eight of ten poses converged to the structural 441 coordinates with slight variation in rotamers, while two of the ten poses flipped the position of the 442 dihydropyridine and carbon ring relative to the structural coordinates (Figure S16.1). With 443 GALigandDock, the lowest ten interface energy poses were at 1.3 Å RMSD, in a similar position and 444 orientation to the low RMSD conformations from RosettaLigand (Figure S16.2). 445 446 447 Figure 5. Nifedipine docking to rabbit CaV1.1 central pore cavity (PDB 6JP5). (A) 448 RosettaLigand (RMSDMin = 0.77 Å, PNear = 0.28) and (B) GALigandDock (RMSDMin = 1.3 Å, PNear = 449 0.29). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs Interface 450 energy distribution of the top 10% of total energy poses. Grey dots: 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 452 GALigandDock in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match 453 the Corey-Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire 454 population. 455 456 Bay K 8644 docking resulted in a RosettaLigand RMSDMin of 0.93 Å and GALigandDock 457 RMSDMin of 0.97 Å (Figure 6, Table 2). The calculated RosettaLigand PNear was 0.66 and the 458 GALigandDock PNear was 0.61 (Table 3). The PNear values, paired with the interaction energy 459 distribution data, indicate well-converged energy funnels. Further, the lowest interaction energy 460 poses with RosettaLigand was within 0.3 Å of the RMSDMin pose (Table S2.1) and for 461 GALigandDock was within 0.2 Å of the RMSDMin pose (Table S2.5). With RosettaLigand, the 462 lowest ten interface energy poses were between 1.2-1.3 Å RMSD, with all ten poses converged to the 463 structural coordinates with slight variation in rotamers (Figure S17.1). With GALigandDock, the 464 lowest ten interface energy poses were at 1.1 Å or 4.6 Å RMSD. Eight of ten poses had 1.1 Å RMSD 465 with slight deviation in position to structural coordinates, while two of ten poses had 4.6 Å RMSD 466 with the dihydropyridine in the correct position, but the carbon ring flipped 180 degrees such that the 467 trifluoromethyl group was pointed in the opposite direction (Figure S17.2). 468 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 18 469 Figure 6. Bay K 8644 docking to rabbit CaV1.1 central pore cavity (PDB 6JP8). (A) 470 RosettaLigand (RMSDMin = 0.93 Å, PNear = 0.66) and (B) GALigandDock (RMSDMin = 0.97 Å, PNear 471 = 0.61). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface 472 energy distribution of the top 10% of total energy poses. Grey dots: 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 474 GALigandDock in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match 475 the Corey-Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire 476 population. 477 Diltiazem docking resulted in a RosettaLigand RMSDMin of 2.5 Å and a GALigandDock 478 RMSDMin of 1.1 Å (Figure 7, Table 2). The calculated PNear suggested no energy funnel convergence 479 for either method (RosettaLigand PNear = 6.4•10-4 and GALigandDock PNear= 7.8•10-6) (Table 3). 480 This lack of energy funnel convergence is exemplified by the interaction energy distributions 481 containing local minima 5-14 Å RMSD away from the ligand’s structural coordinates (Figure 7). 482 With RosettaLigand, the lowest ten interface energy poses were between 5.2-7.6 Å RMSD, with all 483 ten poses at a similar channel depth at the pore center, but with no convergence in local position or 484 rotamer placement. (Figure S20.1). With GALigandDock, the lowest ten interface energy poses were 485 between 10.5-14.5 Å RMSD. Rather than being positioned in the pore center, all ten poses were in a 486 similar channel depth to the structural coordinates but positioned in the fenestration with no 487 convergence in local position or rotamer placement (Figure S20.2). 488 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 19 489 Figure 7. Diltiazem docking to rabbit CaV1.1 central pore cavity (PDB 6JPB). (A) RosettaLigand 490 (RMSDMin = 2.5 Å, PNear = 6.4•10-4) and (B) GALigandDock (RMSDMin = 1.1 Å, PNear = 7.8•10-6). 491 Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 492 distribution of the top 10% of total energy poses. Yellow dots: 2-4 Å, green dots: > 4 Å. Carbon 493 atoms of RosettaLigand molecule are shown in cyan, GALigandDock in dark pink, and native 494 structure coordinates in dark grey. Non-carbon atoms match the Corey-Pauling-Koltun coloring 495 scheme. Black star indicates RMSDMin pose from the entire population. 496 497 Verapamil was previously resolved in a complex with rabbit CaV1.1 in two binding modes with 498 flipped orientations (Zhao et al., 2019a). We evaluated if RosettaLigand and GALigandDock 499 sampling favored a particular orientation. For the first orientation, docking resulted in a 500 RosettaLigand RMSDMin of 1.4 Å and a GALigandDock RMSDMin of 2.0 Å (Figure 8, Table 2). The 501 calculated PNear suggested no energy funnel convergence (RosettaLigand PNear = 4.3•10-3, 502 GALigandDock PNear = 2.8•10-3) (Table 3). This nonexistent energy funnel is evident by local 503 interaction energy minima beyond 4 Å RMSD from the ligand structure coordinates (Figure 8). With 504 RosettaLigand, the lowest ten interface energy poses were between 4.1-9.7 Å RMSD, with all ten 505 poses at a similar channel depth, and some poses converging in local position and rotamer placement 506 at 5.7 Å and 9.2 Å RMSD (Figure S18.1). With GALigandDock, the lowest ten interface energy 507 poses were between 3.8-11.9 Å RMSD. All ten poses were at different channel depths, positioned in 508 the pore center or at the pore center – fenestration interface, and did not converge in local position or 509 rotamer placement (Figure S18.2). 510 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 20 511 Figure 8. Verapamil docking in its first orientation to rabbit CaV1.1 central pore cavity (PDB 512 6JPA). (A) RosettaLigand (RMSDMin = 1.4 Å, PNear = 4.3•10-3) and (B) GALigandDock (RMSDMin = 513 2.0 Å, PNear = 2.8•10-3). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD 514 vs interface energy distribution of the top 10% of total energy poses. Blue dots: 1-2 Å, yellow dots: 515 2-4 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, 516 GALigandDock in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match 517 the Corey-Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire 518 population. 519 520 For the second orientation, docking resulted in a RosettaLigand RMSDMin of 2.2 Å and a 521 GALigandDock RMSDMin of 2.1 Å (Figure 9, Table 2). The calculated PNear values suggest no 522 energy funnel convergence (RosettaLigand PNear = 1.2•10-4, GALigandDock PNear = 1.9•10-5) (Table 523 3). Similar to the first orientation, the interaction energy distribution for the second orientation yields 524 local energy minima greater than 4 Å RMSD from the ligand structure coordinates (Figure 9). With 525 RosettaLigand, the lowest ten interface energy poses were between 8.5-10.8 Å RMSD, with all ten 526 poses at a similar channel depth, and some poses converging in local position, similar to the docking 527 set for the first structural orientation of verapamil (Figure S19.1). With GALigandDock, the lowest 528 ten interface energy poses were between 6.5-10.8 Å RMSD. All ten poses were at different channel 529 depths, positioned in the pore center or at the pore center – fenestration interface, and did not 530 converge in local position or rotamer placement, similar to the docking set for the first structural 531 orientation of verapamil (Figure S19.2). 532 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 21 533 Figure 9. Verapamil docking in its second orientation to rabbit CaV1.1 central pore cavity 534 (PDB 6JPA). (A) RosettaLigand (RMSDMin = 2.2 Å, PNear = 1.2•10-4) and (B) GALigandDock 535 (RMSDMin = 2.1 Å, PNear = 1.9•10-5). Left: transmembrane view. Middle: extracellular view. Right: 536 ligand RMSD vs interface energy distribution of the top 10% of total energy poses. Yellow dots: 2-4 537 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 538 in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-539 Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 540 541 The small molecule Z944 (channel blocker) has 6 rotatable bonds and is resolved bound to 542 human CaV3.1. Docking with RosettaLigand and GALigandDock both resulted in a RMSDMin of 0.73 543 Å and 0.84 Å, respectively (Figure 10, Table 2). Interestingly, the calculated PNear suggests good 544 energy funnel convergence for RosettaLigand (PNear = 0.62), but no energy funnel convergence for 545 GALigandDock (PNear = 1.1•10-4) (Table 3). For RosettaLigand, the interaction energy produced an 546 energy funnel between 1-2 Å, while for GALigandDock produced an energy funnel near 10-15 Å 547 (Figure 10). With RosettaLigand, the lowest ten interface energy poses were between 0.90-1.7 Å 548 RMSD, with all ten poses converged to the structural coordinates with slight variation in rotamers. 549 (Figure S21.1). With GALigandDock, the lowest ten interface energy poses were between 10.6-12.4 550 Å RMSD, with four sets of unique poses identified. The first set, containing three of the ten poses, is 551 in a slightly lower channel depth and in a similar orientation to the structural coordinates, but rotated 552 approximately 180 degrees about the pore center. The second set is one pose and is in a similar 553 location as the first set but is rotated such that the phenyl ring and tertiary butylamine positions are 554 flipped. The third group, containing four poses, is similar to the second set but rotated approximately 555 180 degrees about the pore center. The last group, containing two poses, is in a similar channel depth 556 and location to the structural coordinates, but the phenyl ring and the tertiary butylamine positions 557 are flipped (Figure S21.2). 558 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 22 559 Figure 10. Z944 docking to human CaV3.1 central pore cavity (6KZP). (A) RosettaLigand 560 (RMSDMin = 0.73 Å, PNear = 0.62) and (B) GALigandDock (RMSDMin = 0.84 Å, PNear = 1.1•10-4). 561 Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 562 distribution of the top 10% of total energy poses. Grey dots: 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 564 in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-565 Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 566 567 Flecainide (channel blocker) has 7 rotatable bonds and is resolved bound to rat NaV1.5. 568 RosettaLigand docking resulted in a RMSDMin of 1.2 Å, while GALigandDock resulted in a 569 RMSDMin of 0.94 Å (Figure 11, Table 2). The calculated PNear values suggest no energy funnel 570 convergence with either method (RosettaLigand PNear = 0.12, GALigandDock PNear = 2.9•10-5) 571 (Table 3). The interaction energy distribution from RosettaLigand produced an energy minimum 572 between 2-4 Å and 7-10 Å, while from GALigandDock produced energy minima near 10-12 Å 573 (Figure 11). With RosettaLigand, the lowest ten interface energy poses were between 2.1-4.0 Å 574 RMSD, with all ten poses in the same channel depth as the structural coordinates, but with a tilt such 575 that the trifluoromethyl groups are level in channel depth, rather than the trifluoromethyl group 576 facing the pore being lower in depth. In one pose, the piperidin-2-yl-methylamine was positioned 577 lower into the channel than the rest of the ligand and the structural coordinates, while the rest of the 578 ligand pose was in a similar position to other poses (Figure S23.1). With GALigandDock, the lowest 579 ten interface energy poses were between 9.7-11.9 Å RMSD, with all ten poses positioned at the 580 fenestration in the same channel depth as the structural coordinates, but oriented outside of the pore 581 with one of the trifluoromethyl groups pointing towards the pore center but positioned at the exterior 582 of the channel fenestration (Figure S23.2). 583 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 23 584 Figure 11. Flecainide docking to rat NaV1.5 central pore cavity (PDB 6UZ0). (A) RosettaLigand 585 (RMSDMin = 1.2 Å, PNear = 0.12) and (B) GALigandDock (RMSDMin = 0.94 Å, PNear = 2.9•10-5). Left: 586 transmembrane view. Middle: extracellular view. Right: ligand RMSD vs interface energy 587 distribution of the top 10% of total energy poses. Blue dots: 1-2 Å, yellow dots: 2-4 Å, green dots: > 588 4 Å. Carbon atoms of RosettaLigand molecule is indicated in cyan, GALigandDock in dark pink, and 589 native structure coordinates in dark grey. Non-carbon atoms match the Corey-Pauling-Koltun 590 coloring scheme. Black star indicates RMSDMin pose from the entire population. 591 592 The transient receptor potential vanilloid (TRPV) subfamily of cation channels broadly play roles 593 in thermo-sensation and thermoregulation in response to heat (>53 °C), notably in noxious heat 594 sensing (Nilius, 2007; Gees et al., 2012). TRPV2 channels are widely expressed, with identification 595 in dorsal root ganglion neurons, as well as brain, heart, and smooth muscle tissue, among others 596 (Gees et al., 2012). The TRPV2 channel has been implicated in thermal pain sensing, muscular 597 dystrophy, and cardiomyopathy, among other diseases (Nilius, 2007; Gees et al., 2012). Structurally, 598 TRPV2 channels contain six transmembrane spanning domains and commonly assemble as a 599 homomer with four identical subunits (Zheng and Trudeau, 2015). They resemble canonical voltage-600 gated ion channels with a pore domain and a weak voltage sensor domain, while their gating is 601 regulated by heat and a diverse set of agonists such as 2-aminoethoxydiphenyl borate and cannabidiol 602 (Gees et al., 2012; Zheng and Trudeau, 2015; Pumroy et al., 2019). 603 Cannabidiol (channel activator) has 6 rotatable bonds and is resolved bound to rat TRPV2. 604 Docking with RosettaLigand resulted in a RMSDMin of 1.0 Å, while GALigandDock resulted in a 605 RMSDMin of 0.90 Å (Figure 12, Table 2). The calculated PNear values suggest no energy funnel 606 convergence for both methods (RosettaLigand PNear = 4.3•10-2, GALigandDock PNear = 3.8•10-3) 607 (Table 3). For RosettaLigand, the interaction energy distribution did not demonstrate a clear 608 minimum, with minima seen between 2-8 Å (Figure 12). For GALigandDock, there is a clear 609 minimum between 6-7 Å (Figure 12). With RosettaLigand, the lowest ten interface energy poses 610 were between 2.8-7.0 Å RMSD, with two groups of conformations in the same channel depth and 611 interfaced with the S5 and S6 helices of adjacent TRPV2 monomers like the structural coordinates. In 612 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 24 one group, cannabidiol is positioned parallel to the S6 helical segment, with the cyclohexene at the 613 lowest channel depth. In the second group, the poses are positioned in the same orientation as the 614 structural coordinates, with the pentyl chain facing away from the pore center, but the overall 615 position of cannabidiol is slightly elevated relative to the structural coordinates (Figure S22.1). With 616 GALigandDock, the lowest ten interface energy poses were either 6.4 Å or 6.7 Å RMSD, with all ten 617 poses positioned at the fenestration in the same channel depth and interfacing with the S5 and S6 618 helices of adjacent TRPV2 monomers like the structural coordinates, but with the ligand rotated 619 around the cyclohexene plane such that the pentyl chain is oriented parallel to the S6 helix, rather 620 than pointing outward (Figure S22.2). 621 622 Figure 12. Cannabidiol docking to rat TRPV2 central pore cavity (PDB 6U88). (A) 623 RosettaLigand (RMSDMin = 1.0 Å, PNear = 4.3•10-2) and (B) GALigandDock (RMSDMin = 0.90 Å, 624 PNear = 3.8•10-3). Left: transmembrane view. Middle: extracellular view. Right: ligand RMSD vs 625 interface energy distribution of the top 10% of total energy poses. Blue dots: 1-2 Å, yellow dots: 2-4 626 Å, green dots: > 4 Å. Carbon atoms of RosettaLigand molecule are shown in cyan, GALigandDock 627 in dark pink, and native structure coordinates in dark grey. Non-carbon atoms match the Corey-628 Pauling-Koltun coloring scheme. Black star indicates RMSDMin pose from the entire population. 629 630

631 Previous studies have underscored the importance of ligand docking methods for generating ion 632 channel structure-based hypotheses (Yang et al., 2013; Shim et al., 2019). Furthermore, ligand 633 docking methods when combined with high-resolution structures can aid in rational drug design 634 (Wang et al., 2007; Liu et al., 2018; Wulff et al., 2019; Maia et al., 2020). Notably, RosettaLigand 635 has been extensively used to predict the molecular mechanisms of ligand – ion channel interactions 636 (Yang et al., 2013; Yang et al., 2015; Nguyen et al., 2017; Nguyen et al., 2019; Shim et al., 2019; 637 Craig 2nd et al., 2022; Vu et al., 2020; Maly et al., 2022; Pumroy et al., 2022). While GALigandDock 638 has not yet been tested on ion channel structure-function relationships, it utilizes a new generalized 639 energy function tailored for small molecules while sampling ligand conformations with a genetic 640 algorithm, making it an attractive complement to RosettaLigand. 641 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 25 Indeed, this study demonstrates that high-resolution structures paired with RosettaLigand and 642 GALigandDock can be useful tools for formulating structural hypotheses and predicting binding 643 modes for drug discovery. Notably, both Rosetta methods can produce ligand poses near the native 644 ligand structure coordinates. Using a standard 2 Å RMSDMin cutoff (Maia et al., 2020; Park et al., 645 2021), the RosettaLigand docking set yielded an average RMSDMin of 1.2 ± 0.6 Å, and the 646 GALigandDock docking set yielded an average RMSDMin of 1.1 ± 0.5 Å (Table 2). However, the 647 ability to discriminate the RMSDMin pose remains challenging, and highlights important practical 648 considerations when applying Rosetta ligand docking methods to a chosen ion channel target. 649 When performing Rosetta ligand docking, the features of the ligand conformer library, the size of 650 the receptor site, and prior knowledge of critical functional residues, need to be considered to 651 determine the appropriate amount of pose sampling. For our study, we generated either 20,000 poses 652 or 100,000 poses per docking run. While we did notice a statistically significant difference in 653 RMSDMin between 20,000 and 100,000 total poses, the difference was within a few sub-Angstroms. 654 Furthermore, another RosettaLigand benchmarking study of the CASF-2016 dataset generated 1,000 655 total poses and was able to sufficiently sample a RMSDMin below 2 Å (Su et al., 2019; Smith and 656 Meiler, 2020). Thus, for ion channel docking, we suggest sampling in the ones to tens of thousands, 657 especially when docking in large receptor sites like the pore-forming domain. 658 Overall, our docking data was able to achieve a PNear greater than 0.5 for two test cases for each 659 docking method when using the top 10% of total-scoring poses. Both methods achieved PNear greater 660 than 0.5 for Bay K8644 to the rabbit CaV1.1 central pore cavity, while RosettaLigand achieved PNear 661 greater than 0.5 for Z944 to the human CaV3.1 central pore cavity, and GALigandDock achieved 662 PNear greater than 0.5 for GX-936 docking to hNaV1.7-NaVAb voltage sensing domain four (Table 3). 663 The reasons for a low PNear value in most test cases are possibly due to 1) improper scoring by 664 Rosetta score functions to discriminate native-like poses by energy, 2) multiple favorable ligand 665 binding states in the receptor site that have not been structurally resolved, and/or 3) insufficient pose 666 sampling. 667 Currently, we suggest that Rosetta score functions are unable to sufficiently score near-native 668 poses accurately in ion channel docking. For each docking run, comparing the RMSDMin pose and the 669 lowest interface energy pose indicates that the RMSDMin pose is under-scored (Tables S2.1-2.5). 670 Conversely, a previous RosettaLigand study using the CASF-2016 dataset (containing 285 crystal 671 structures of protein – ligand complexes with an overall resolution < 2.5 Å and an R-factor < 0.25) 672 and 1,000 total poses per docking run was able to frequently achieve PNear values between 0.8 and 1.0 673 (Su et al., 2019; Smith and Meiler, 2020), suggesting that the Rosetta score functions can be utilized 674 for other docking studies, but should be verified for accuracy with test cases. Furthermore, the 675 CASF-2016 dataset does not contain ion channels, while the sample size per docking run is 676 sufficiently lower in the RosettaLigand docking study using CASF-2016 than those in our study, 677 suggesting that sampling is not the predominant issue, but rather favorable scoring of near-native 678 poses is the primary issue. 679 For the voltage-sensing domains, GX-936 in complex with the VSD4 of hNaV1.7-NaVAb chimera 680 (PDB: 5EK0) was the only small-molecule docked yet that was consistently near or below 1 Å 681 RMSDMin for each docking set regardless of the method (Tables S1.1-1.2). This may be due to the 682 VSD4 receptor site being the narrowest binding pocket tested, thereby limiting the number of binding 683 configurations. This suggests that VSDs are generally well-suited for Rosetta ligand docking since 684 the receptor is constrained to allow shape-complementary between ligand and protein while reducing 685 the required sampling compared to pore-forming receptor sites. Notably, the Rosetta ligand docking 686

employed did not use any implicit membrane parameters, while GX-936 in a biologically 687 realistic context is partially exposed to a lipid head group (Ahuja et al., 2015). Thus, in this case, 688 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 26 docking without an implicit membrane energy function was still able to achieve the experimental 689 structurally resolved position, however, artifactual, low interface energy poses where GX-936 would 690 be overlapping the phospholipid space are present (Figure S13.1-13.2). Further, GALigandDock was 691 consistently able to achieve a PNear > 0.5 when using a padding of 2.0 Å, 4.0 Å, or 7.0 Å, suggesting 692 that ligands docked to VSDs could be screened by energy funnel convergence (Tables S4.5-4.10). 693 More small molecules structurally resolved and docked to the VSD are needed to validate this 694 observation. 695 For the pore-forming domain, TTX and STX were docked to hNaV1.7, as they are canonical 696 small-molecule channel blockers used when studying the NaV family. TTX and STX had the fewest 697 rotatable bonds in this study, suggesting that little conformer sampling would be needed to achieve a 698 near-native pose. Both TTX and STX docking achieved sub-Angstrom RMSDMin, except STX 699 docking with GALigandDock at a padding of 7 Å (RMSDMin =2.3 Å, padding 7 Å/100,000 poses vs 700 RMSDMin = 1.0 Å, padding 4 Å/100,000 poses; Table S1.2). However, all docking methods 701 employed were unable to achieve a PNear > 0.5, suggesting that 1) STX and TTX could have 702 alternative binding conformations, 2) the docking methods employed have wrongfully scored 703 alternative low-energy conformations, or 3) that inherent loop flexibility in the pore-forming domain 704 is a requirement for docking an energy-optimized, induced-fit conformation. It should be noted that 705 the hNaV1.7 selectivity filter is lined with polar residues that could contribute to a vast hydrogen 706 bonding network, in addition to water hydrating the selectivity filter and the filter opening. Thus, it is 707 unclear if the lowest energy poses, with potential salt bridge interactions, are possible alternate 708 binding modes. It is thus possible that the sum of hydrogen bonding interactions and water-bridging 709 effects could bias the energetic potential to a certain state, such as the one structurally resolved. Thus, 710 further experimental characterization is needed to test these structural hypotheses. 711 712 For the central pore cavity, while there should be little to no lipid interactions (excluding the 713 fenestration regions), docking of the central pore cavity produced the most pose variability, likely 714 due to it being the widest receptor site compared to the pore-forming domain and voltage sensing 715 domain. For example, the second orientation of verapamil positioned primarily in the central cavity 716 of rabbit CaV1.1 achieved a RMSDMin of 1.4 Å for RosettaLigand and 2.0 Å for GALigandDock 717 (Table 2). However, Z944 bound to hCaV3.1, with the wide aromatic end of Z944 in the narrow 718 fenestration and the narrow amide end in the wide cavity achieved a RMSDMin of 0.73 Å for 719 RosettaLigand and 0.84 Å for GALigandDock (Table 2). Further, RosettaLigand was able to achieve 720 a PNear value greater than 0.5 in some docking conditions for Z944 (Table 3; Table S4.1-4.4), 721 suggesting that the Rosetta docking methods could prove useful when docking similar ligands that 722 bridge between the fenestration and central pore, and target the narrow fenestration with an aromatic 723 moiety. 724 Small molecules docked in the central cavity were bound to the central pore (6JPA, orientation 2, 725 6JPB), the channel fenestration region (6JP5, 6JP8, 6U88), or both regions (6JPA orienation1, 6KZP, 726 6UZ0). It appears that those bound in the channel fenestration produced a lower RMSDMin, those 727 primarily bound in the pore center produced the largest RMSDMin, while those bound to both regions 728 produced intermediate RMSDMin values, with Z944 docked to human CaV3.1 central pore cavity 729 (6KZP) being an exception. The same trend is roughly observed for PNear, where fenestration-only 730 cases have PNear values orders of magnitude greater than pore center cases (Table S5). Due to the 731 limited number of cases for each sub-classification, further studies will need to be performed to 732 assess a correlation. 733 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 27 It appears that molecules with predominantly planar aromatic rings and space-filling, “bulky” 734 structures targeting the fenestration can be scored and ranked effectively; both Bay K 8644 (6JP8) 735 and Z944 (6KZP) possess these features in contrast to other small molecules, that possess aromatic or 736 aliphatic rings but are generally more linear and “floppy” targeting the fenestration, like flecainide 737 (6UZ0) and the first orientation of verapamil (6JPA). (Table S5). This general trend in increased 738 RMSDMin for non-aromatic containing small molecules when the receptor site is solely the central 739 pore could be due to 1) greater ligand flexibility within a larger area, which requires increased 740 sampling, 2) a bulkier ligand having fewer conformations and orientations relative to the channel, 3) 741 the absence of explicit water molecules and ions that would crowd the pore or directly bind to the 742 ligand (Tikhonov and Zhorov, 2008). Furthermore, if the ligand is expected to bind to a central cavity 743 motif present on multiple subunits, then post-hoc tools implementing symmetry to discriminate 744 unique binding modes will be necessary to calculate an appropriate adjusted RMSDMin. 745 746

747 In this study, we aimed to assess if 1) the RosettaLigand and GALigandDock molecular docking 748 methodologies can recapitulate structurally resolved ion channel – small-molecule binding 749 orientations with known pharmacological significance, 2) if their scoring functions can be used to 750 accurately rank small molecule binding orientation for the purpose of blindly screening small 751 molecules, and 3) if there are practical considerations when docking to specific domains of ion 752 channels. With 2.0 Å RMSDMin as a performance cutoff (Maia et al., 2020; Park et al., 2021), our 753

demonstrate that both RosettaLigand and GALigandDock can frequently sample the 754 experimentally resolved ligand binding mode with less than 2.0 Å RMSD. However, when using an 755 estimate of the Boltzmann probability for energy “funnel-likeness” (PNear) as a scoring function 756 assessment, we currently perceive Rosetta score functions as unable to sufficiently score near-native 757 poses accurately in ion channel docking; from this study, small molecules targeting voltage-sensing 758 domains and bulky small molecules primarily composed of aromatic rings targeting fenestration 759 regions appear to be most suited for score-based ranking. Thus, when performing an ion channel 760 virtual drug discovery campaign, special consideration should be given to sufficient pose sampling to 761 account for multiple rotameric and conformational states, to identify the size of the sampling required 762 for sufficient interface scoring of the receptor site, to identify the appropriate state of the ion channel, 763 to identify inherent channel flexibility that could influence ligand binding, and to potentially identify 764 specific chemical functional groups known experimentally to influence binding to the target when 765 selecting candidate conformations. 766 767

768 We thank members of the VY-Y, FD, and HW laboratories for helpful discussions. We also thank 769 Igor Vorobyov for helpful discussions and Jerome Manera for help with generating 2d structures of 770 ligands. The work in VY-Y lab was supported by NIH grants R61NS127285, R01HL128537, 771 R01HL159304, R01GM132110, R56NS9706, and R01NS128180. BJH was supported by NIH F31 772 Predoctoral Fellowship F31NS124337 and NIH T32 Predoctoral Fellowship GM007377. We are 773 grateful to OpenEye Scientific for OpenEye Academic License. 774 775 776 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 28

777 Ahuja, S., Mukund, S., Deng, L., Khakh, K., Chang, E., Ho, H., et al., and Payandeh, J. (2015) 778 Structural basis of NaV1.7 inhibition by an isoform-selective small-molecule antagonist. Science. 779 350(6267), aac5464. https://doi.org/10.1126/science.aac5464 780 781 Alexander, S. P., Kelly, E., Mathie, A., Peters, J. A., Veale, E. L., Armstrong, J. F., et al., and Verri, 782 T. (2021). The concise guide to pharmacology 2021/22: transporters. Br. J. Pharmacol. 178, S412–783 S513. https://doi.org/10.1111/bph.15543 784 785 Alford, R. F., Leaver-Fay, A., Jeliazkov, J. R., O'Meara, M. J., DiMaio, F. P., Park, H., Shapovalov, 786 M. V., Renfrew, P. D., Mulligan, V. K., Kappel, K., Labonte, J. W., Pacella, M. S., Bonneau, R., 787 Bradley, P., Dunbrack, R. L., Jr, Das, R., Baker, D., Kuhlman, B., Kortemme, T., and Gray, J. J. 788 (2017). The Rosetta all-atom energy function for macromolecular modeling and design. J. Chem. 789 Theory Comput., 13(6), 3031–3048. https://doi.org/10.1021/acs.jctc.7b00125 790 791 Bagal, S. K., Brown, A. D., Cox, P. J., Omoto, K., Owen, R. M., Pryde, D. C., Sidders, B., Skerratt, 792 S. E., Stevens, E. B., Storer, R. I., and Swain, N. A. (2013). Ion channels as therapeutic targets: a 793 drug discovery perspective. J. Med. Chem., 56(3), 593–624. https://doi.org/10.1021/jm3011433 794 795 Bhardwaj, G., Mulligan, V. K., Bahl, C. D., Gilmore, J. M., Harvey, P. J., Cheneval, O., et al., and 796 Baker, D. (2016). Accurate de novo design of hyperstable constrained peptides. Nature, 538(7625), 797 329–335. https://doi.org/10.1038/nature19791 798 799 Bennett, D. L., Clark, A. J., Huang, J., Waxman, S. G., and Dib-Hajj, S. D. (2019). The Role of 800 Voltage-Gated Sodium Channels in Pain Signaling. Physiological reviews, 99(2), 1079–1151. 801 https://doi.org/10.1152/physrev.00052.2017 802 803 Berman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., Shindyalov, I. N., 804 Bourne, P. E. (2000). The Protein Data Bank. Nucleic Acids Res. 28(1), 235-242. 805 https://doi.org/10.1093/nar/28.1.235 806 807 Case, D. A., Aktulga, H. M., Belfon, K., Ben-Shalom, I. Y., Brozell, S. R., Cerutti, D. S., et al., and 808 Kollman, P. A. (2021). Amber 2021. University of California, San Francisco. 809 810 Catterall, W. A. (1995) Structure and function of voltage-gated ion channels. Annu. Rev. Biochem 64, 811 493–531. https://doi.org/10.1146/annurev.bi.64.070195.002425 812 813 Catterall, W. A. (2011). Voltage-gated calcium channels. Cold Spring Harb. Perspect. Biol. 3(8), 814 a003947. https://doi.org/10.1101/cshperspect.a003947 815 816 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 29 Catterall, W. A., Lenaeus, M. J., and Gamal El-Din, T. M. (2020) Structure and pharmacology of 817 voltage-gated sodium and calcium channels. Annu. Rev. Pharmacol. Toxicol. 60, 133–154. 818 https://doi.org/10.1146/annurev-pharmtox-010818-021757 819 820 Combs, S., Kaufmann, K., Field, J. R., Blakely, R. D., and Meiler, J. (2011). Y95 and E444 821 interaction required for high-affinity S-citalopram binding in the human serotonin transporter. ACS 822 Chem Neurosci., 2(2), 75–81. https://doi.org/10.1021/cn100066p 823 824 Craig, R. A., 2nd, Garrison, C. E., Nguyen, P. T., Yarov-Yarovoy, V., & Du Bois, J. (2020). 825 Veratridine: A janus-faced modulator of voltage-gated sodium ion channels. ACS Chem Neurosci, 826 11(3), 418–426. https://doi.org/10.1021/acschemneuro.9b00621 827 828 Davis, I. W., Baker, D. (2009). RosettaLigand docking with full ligand and receptor flexibility. J. 829 Mol. Biol. 385(2), 381-392. https://doi.org/10.1016/j.jmb.2008.11.010 830 831 DeLuca, S., Khar, K., Meiler, J. (2015). Fully flexible docking of medium sized ligand libraries with 832 RosettaLigand. PLoS One. 10(7), e0132508. https://doi.org/10.1371/journal.pone.0132508 833 834 Gees, M., Owsianik, G., Nilius, B., and Voets, T. (2012). TRP channels. Compr. Physiol., 2(1), 563–835 608. https://doi.org/10.1002/cphy.c110026 836 837 Feske, S., Wulff, H., and Skolnik, E. Y. (2015). Ion channels in innate and adaptive immunity. Annu. 838 Rev. Immunol., 33, 291–353. https://doi.org/10.1146/annurev-immunol-032414-112212 839 840 Hagen, N. A., Cantin, L., Constant, J., Haller, T., Blaise, G., Ong-Lam, M., du Souich, P., Korz, W., 841 and Lapointe, B. (2017). Tetrodotoxin for moderate to severe cancer-related pain: a multicentre, 842 randomized, double-blind, placebo-controlled, parallel-design trial. Pain Res. Manag. 2017:7212713. 843 https://doi.org/10.1155/2017/7212713 844 845 Halgren, T. A. (1996). Merck molecular force field. I. Basis, form, scope, parameterization, and 846 performance of MMFF94. J. Comput. Chem. 17, 490-519. doi: 10.1002/(SICI)1096-847 987X(199604)17:5/63.0.CO;2-P 848 849 Halgren, T. A. (1996a). Merck molecular force field. II. MMFF94 van der Waals and electrostatic 850 parameters for intermolecular interactions. J. Comput. Chem.17520–552. doi: 10.1002/(SICI)1096-851 987X(199604)17:5/63.0.CO;2-W 852 853 Halgren, T. A. (1996b). Merck molecular force field. III. Molecular geometries and vibrational 854 frequencies for MMFF94. J. Comput. Chem. 17, 553–586. doi: 10.1002/(SICI)1096-855 987X(199604)17:5/63.0.CO;2-T 856 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 30 857 Halgren, T. A. (1996c). Merck molecular force field. V. Extension of MMFF94 using experimental 858 data, additional computational data, and empirical rules. J. Comput. Chem. 17:616–641. doi: 859 10.1002/(SICI)1096-987X(199604)17:5/63.0.CO;2-X 860 861 Halgren, T. A. and Nachbar, R. B. (1996). Merck molecular force field. IV. Conformational energies 862 and geometries for MMFF94. J. Comput. Chem. 17, 587–615. doi: 10.1002/(SICI)1096-863 987X(199604)17:5/6%3C587::AID-JCC4%3E3.0.CO;2-Q 864 865 Halgren, T. A. (1999a). MMFF VI. MMFF94s option for energy minimization studies. J. Comput. 866 Chem. 20, 720–729. doi: 10.1002/(SICI)1096-987X(199905)20:73.0.CO;2-X 867 868 Halgren, T. A. (1999b). MMFF VII. Characterization of MMFF94, MMFF94s, and other widely 869 available force fields for conformational energies and for intermolecular-interaction energies and 870 geometries. J. Comput. Chem. 20, 730–748. doi: 10.1002/(SICI)1096-987X(199905)20:73.0.CO;2-T 872 873 Hanwell, M. D., Curtis, D. E., Lonie, D. C., Vandermeersch, T., Zurek, E., Hutchison, G. R. (2012). 874 Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. J. 875 Cheminform. 4(1), 17. https://doi.org/10.1186/1758-2946-4-17. 876 877 Hawkins, P. C. D., Skillman, A. G., Warren, G. L., Ellingson, B. A., Stahl, M. T. (2010). Conformer 878 generation with OMEGA: Algorithm and validation using high quality structures from the Protein 879 Databank and Cambridge Structural Database. J. Chem. Inf. Model. 50(4), 572-584. 880 https://doi.org/10.1021/ci100031x 881 882 Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., et al., 883 and Oliphant, T. E. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. 884 https://doi.org/10.1038/s41586-020-2649-2 885 886 Hille, B. (1977). Local anesthetics: hydrophilic and hydrophobic pathways for the drug-receptor 887 reaction. J. Gen. Physiol. 69, 497–515. https://doi.org/10.1085/jgp.69.4.497 888 889 Hille, B. (2001). Ion channels of excitable membranes (3rd ed., pp. 503–537). Sunderland, MA, 890 USA. Sinauer Associates, Inc., Sunderland. 891 892 Hinman, A., and Du Bois, J. (2003). A stereoselective synthesis of (-)-tetrodotoxin. J. Am. Chem. 893 Soc., 125(38), 11510–11511. https://doi.org/10.1021/ja0368305 894 895 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 31 Hockerman, G. H., Peterson, B. Z., Johnson, B. D., and Catterall, W. A. (1997). Molecular 896 determinants of drug binding and action on L-type calcium channels. Annual review of 897 pharmacology and toxicology, 37, 361–396. https://doi.org/10.1146/annurev.pharmtox.37.1.361 898 899 Jiang, D., Shi, H., Tonggu, L., Gamal El-Din, T. M., Lenaeus, M. J., Zhao, Y., Yoshioka, C., Zheng, 900 N., Catterall, W. A. (2020) Structure of the cardiac sodium channel. Cell. 180, 122-134.e10. 901 https://doi.org/10.1016/j.cell.2019.11.041 902 903 Kaufmann, K. W., and Meiler, J. (2012). Using RosettaLigand for small molecule docking into 904 comparative models. PloS One, 7(12), e50769. https://doi.org/10.1371/journal.pone.0050769 905 906 Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., Li, Q., Shoemaker, B. A., Thiessen, P. A., 907 Yu, B., Zaslavsky, L., Zhang, J., Bolton, E. E. (2023). PubChem 2023 update. Nucleic Acids Res. 908 51(D1):D1373-D1380. https://doi.org/10.1093/nar/gkac956. 909 910 Kwong, K., and Carr, M. J. (2015). Voltage-gated sodium channels. Curr. Opin. Pharmacol. 22, 911 131–139. https://doi.org/10.1016/j.coph.2015.04.007 912 913 Lemmon, G., Kaufmann, K., and Meiler, J. (2012). Prediction of HIV-1 protease/inhibitor affinity 914 using RosettaLigand. Chem Biol Drug Des., 79(6), 888–896. https://doi.org/10.1111/j.1747-915 0285.2012.01356.x 916 917 Liu, H., Hofmann, J., Fish, I., Schaake, B., Eitel, K., Bartuschat, A., Kaindl, J., Rampp, H., Banerjee, 918 A., Hübner, H., Clark, M. J., Vincent, S. G., Fisher, J. T., Heinrich, M. R., Hirata, K., Liu, X., 919 Sunahara, R. K., Shoichet, B. K., Kobilka, B. K., and Gmeiner, P. (2018). Structure-guided 920 development of selective M3 muscarinic acetylcholine receptor antagonists. Proc Natl Acad Sci USA, 921 115(47), 12046–12050. https://doi.org/10.1073/pnas.1813988115 922 923 Maia EHB, Assis LC, de Oliveira TA, da Silva AM, Taranto AG. Structure-Based Virtual Screening: 924 From Classical to Artificial Intelligence. Front. Chem. 2020 Apr 28;8:343. 925 https://doi.org/10.3389/fchem.2020.00343. 926 927 Maly, J., Emigh, A.M., DeMarco, K.R., Furutani, K., Sack, J.T., Clancy, C.E., Vorobyov, I., Yarov-928 Yarovoy, V. (2022) Structural modeling of the hERG potassium channel and associated drug 929 interactions. Front Pharmacol. Sep 16;13:966463. https://doi.org/10.3389/fphar.2022.966463. 930 931 McCormack, K., Santos, S., Chapman, M. L., Krafte, D. S., Marron, B. E., West, C. W., Krambis, M. 932 J., Antonio, B. M., Zellmer, S. G., Printzenhoff, D., Padilla, K. M., Lin, Z., Wagoner, P. K., Swain, 933 N. A., Stupple, P. A., de Groot, M., Butt, R. P., & Castle, N. A. (2013). Voltage sensor interaction 934 site for selective small molecule inhibitors of voltage-gated sodium channels. Proc. Natl. Acad. Sci. 935 U S A. 110(29), E2724–E2732. https://doi.org/10.1073/pnas.1220844110 936 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 32 937 Meiler, J., and Baker, D. (2006). RosettaLigand: protein – small molecule docking with full side-938 chain flexibility. Proteins. 65(3), 538–548. https://doi.org/10.1002/prot.21086 939 940 Moretti, R., Bender, B. J., Allison, B., and Meiler, J. (2016). Rosetta and the design of ligand binding 941 sites. Methods Mol. Biol. (Clifton, N.J.), 1414, 47–62. https://doi.org/10.1007/978-1-4939-3569-7_4 942 943 Mulligan, V. K., Workman, S., Sun, T., Rettie, S., Li, X., Worrall, et al., Baker, D. (2021). 944 Computationally designed peptide macrocycle inhibitors of New Delhi metallo-β-lactamase 1. Proc. 945 Natl. Acad. of Sci. USA, 118(12), e2012800118. https://doi.org/10.1073/pnas.2012800118 946 947 Nanou, E., and Catterall, W. A. (2018). Calcium channels, synaptic plasticity, and neuropsychiatric 948 disease. Neuron 98(3), 466–481. https://doi.org/10.1016/j.neuron.2018.03.017 949 950 Nerbonne, J. M., and Kass, R. S. (2005). Molecular physiology of cardiac repolarization. Physiol. 951 Rev. 85(4), 1205–1253. https://doi.org/10.1152/physrev.00002.2005 952 953 Nguyen, H. M., Singh, V., Pressly, B., Jenkins, D. P., Wulff, H., & Yarov-Yarovoy, V. (2017). 954 Structural insights into the atomistic mechanisms of action of small molecule inhibitors targeting the 955 KCa3.1 channel pore. Mol. Pharmacol., 91(4), 392–402. https://doi.org/10.1124/mol.116.108068 956 957 Nguyen, P. T., DeMarco, K. R., Vorobyov, I., Clancy, C. E., & Yarov-Yarovoy, V. (2019). Structural 958 basis for antiarrhythmic drug interactions with the human cardiac sodium channel. Proc. Natl. Acad. 959 of Sci. USA, 116(8), 2945–2954. https://doi.org/10.1073/pnas.1817446116 960 961 Nguyen, P. T., and Yarov-Yarovoy, V. (2022). Towards Structure-Guided Development of Pain 962 Therapeutics Targeting Voltage-Gated Sodium Channels. Frontiers in pharmacology, 13:842032. 963 https://doi.org/10.3389/fphar.2022.842032 964 965 Nilius B. (2007). TRP channels in disease. Biochim. Biophys. Acta., 1772(8), 805–812. 966 https://doi.org/10.1016/j.bbadis.2007.02.002 967 968 Nivón, L.G., Moretti, R., Baker, D. (2013). A Pareto-optimal refinement method for protein design 969 scaffolds. PLoS One. 8(4), e59004. https://doi.org/10.1371/journal.pone.0059004 970 971 Park, H., Zhou, G., Baek, M., Baker, D., and DiMaio, F. (2021). Force field optimization guided by 972 small molecule crystal lattice data enables consistent sub-Angstrom protein-ligand docking. J. Chem. 973 Theory Comput 17(3), 2000-2010. https://doi.org/10.1021/acs.jctc.0c01184 974 975 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 33 Pajouhesh, H., Beckley, J. T., Delwig, A., Hajare, H. S., Luu, G., Monteleone, D., Zhou, X., Ligutti, 976 J., Amagasu, S., Moyer, B. D., Yeomans, D. C., Du Bois, J., and Mulcahy, J. V. (2020). Discovery of 977 a selective, state-independent inhibitor of NaV1.7 by modification of guanidinium toxins. Sci Rep. 978 10(1):14791. https://doi.org/10.1038/s41598-020-71135-2 979 980 Pumroy, R. A., Samanta, A., Liu, Y., Hughes, T. E., Zhao, S., Yudin, Y., Rohacs, T., Han, S., 981 Moiseenkova-Bell, V. Y. (2019). Molecular mechanism of TRPV2 channel modulation by 982 cannabidiol. Elife. 8e48792. https://doi.org/10.7554/eLife.48792 983 984 Pumroy, R.A., Protopopova, A.D., Fricke, T.C., Lange, I.U., Haug, F.M., Nguyen, P.T., Gallo, P.N., 985 Sousa, B.B., Bernardes, G.J.L., Yarov-Yarovoy, V., Leffler, A., Moiseenkova-Bell, V.Y. (2022). 986 Structural insights into TRPV2 activation by small molecules. Nat Commun. 13(1):2334. 987 https://doi.org/.1038/s41467-022-30083-3. 988 989 Salomon-Ferrer, R., Case, D. A., Walker, R. C. (2012). An overview of the Amber biomolecular 990 simulation package. WIREs Comput. Mol. Sci. 3, 198-210. https://doi.org/10.1002/wcms.1121 991 992 Shen, H., Liu, D., Wu, K., Lei, J., and Yan, N. (2019). Structures of human NaV1.7 channel in 993 complex with auxiliary subunits and animal toxins. Science. 363(6433), 1303–1308. 994 https://doi.org/10.1126/science.aaw2493 995 996 Shim, H., Brown, B. M., Singh, L., Singh, V., Fettinger, J. C., Yarov-Yarovoy, V., and Wulff, H. 997 (2019). The Trials and Tribulations of Structure Assisted Design of KCa Channel Activators. Front 998 Pharmacol., 10, 972. https://doi.org/10.3389/fphar.2019.00972 999 1000 Smith, S. T., Meiler, J. (2020). Assessing multiple score function in Rosetta for drug discovery. PLoS 1001 One. 15(10), e0240450. https://doi.org/10.1371/journal.pone.0240450 1002 1003 Stevens, M., Peigneur, S., and Tytgat, J. (2011). Neurotoxins and their binding areas on voltage-1004 gated sodium channels. Front Pharmacol., 2(71). https://doi.org/10.3389/fphar.2011.00071 1005 1006 Su, M., Yang, Q., Du, Y., Feng, G., Liu, Z., Li, Y., & Wang, R. (2019). Comparative Assessment of 1007 Scoring Functions: The CASF-2016 Update. J. Chem. Inf. Model., 59(2), 895–913. 1008 https://doi.org/10.1021/acs.jcim.8b00545 1009 1010 Tikhonov, D. B., and Zhorov, B. S. (2008). Molecular modeling of benzothiazepine binding in the L-1011 type calcium channel. J. Biol. Chem., 283(25), 17594–17604. 1012 https://doi.org/10.1074/jbc.M800141200 1013 1014 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 34 Thomas-Tran, R., Du Bois, J. (2016). Mutant cycle analysis with modified saxitoxins reveals specific 1015 interactions critical to attaining high-affinity inhibition of hNaV1.7. Proc. Natl. Acad. Sci. U S A. 1016 113(21), 5856-61. https://doi.org/10.1073/pnas.1603486113 1017 1018 Vallat, R. (2018). Pingouin: statistics in Python. Journal of Open Source Software, 3(31), 1026. 1019 https://doi.org/10.21105/joss.01026 1020 1021 Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al., and 1022 van Mulregt, P. (2020). SciPy 1.0: fundamental algorithms for scientific computing in 1023 Python. Nature methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2 1024 1025 Vu, S., Singh, V., Wulff, H., Yarov-Yarovoy, V., Zheng, J. (2020). New capsaicin analogs as 1026 molecular rulers to define the permissive conformation of the mouse TRPV1 ligand-binding pocket. 1027 Elife. Nov 9;9:e62039. https://doi.org/10.7554/eLife.62039. 1028 1029 Wang, S. Y., Tikhonov, D. B., Mitchell, J., Zhorov, B. S., and Wang, G. K. (2007). Irreversible block 1030 of cardiac mutant Na+ channels by batrachotoxin. Channels (Austin), 1(3), 179–188. 1031 https://doi.org/10.4161/chan.4437 1032 1033 Wulff, H., Christophersen, P., Colussi, P., Chandy, K. G., and Yarov-Yarovoy, V. (2019). Antibodies 1034 and venom peptides: new modalities for ion channels. Nat. Rev. Drug Discov. 18(5), 339-357. 1035 https://doi.org/10.1038/s41573-019-0013-8 1036 1037 Yang, T., Smith, J. A., Leake, B. F., Sanders, C. R., Meiler, J., and Roden, D. M. (2013). An 1038 allosteric mechanism for drug block of the human cardiac potassium channel KCNQ1. Mol 1039 Pharmacol., 83(2), 481–489. https://doi.org/10.1124/mol.112.081513 1040 1041 Yang, F., Xiao, X., Cheng, W., Yang, W., Yu, P., Song, Z., Yarov-Yarovoy, V., & Zheng, J. (2015). 1042 Structural mechanism underlying capsaicin binding and activation of the TRPV1 ion channel. Nat. 1043 Chem. Biol., 11(7), 518–524. https://doi.org/10.1038/nchembio.1835 1044 1045 Zhang, Y., Wang, K., and Yu, Z. (2020). Drug Development in Channelopathies: Allosteric 1046 Modulation of Ligand-Gated and Voltage-Gated Ion Channels. J. Med. Chem., 63(24), 15258–15278. 1047 https://doi.org/10.1021/acs.jmedchem.0c01304 1048 1049 Zheng, J., and Trudeau, M.C. (2015). Chapter 28. TRPV channels. In Handbook of Ion Channels (1st 1050 ed., pp. 427-431). Boca Raton, FL, USA. CRC Press. https://doi.org/10.1201/b18027 1051 1052 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint Running Title 35 Zhao, Y., Huang, G., Wu, J., Wu, Q., Gao, S., Yan, Z., Lei, J., and Yan, N. (2019a). Molecular basis 1053 for ligand modulation of a mammalian voltage-gated Ca2+ channel. Cell. 177(6), 1495-1506.e12. 1054 https://doi.org/10.1016/j.cell.2019.04.043 1055 1056 Zhao, Y., Huang, G., Wu, Q., Wu, K., Li, R., Lei, J., Pan, X., Yan, N. (2019b). Cryo-EM structures 1057 of apo and antagonist-bound human Cav3.1. Nature. 576(7787), 492–497. 1058 https://doi.org/10.1038/s41586-019-1801-3 1059 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 3, 2024. ; https://doi.org/10.1101/2024.04.02.587818doi: bioRxiv preprint